和の公式

ΣX^Nを求める

PUBLIC NUMERIC MAXLEVEL
LET MAXLEVEL=20
DIM Y(MAXLEVEL+1),X(MAXLEVEL+1)
FOR I=1 TO MAXLEVEL
PRINT "ΣX^";STR$(I);"=";
CALL SUM(X,I)
CALL DISPLAY(X,I+1)
NEXT I
END

EXTERNAL SUB SUM(X(),N)
!'Σ(X^M)=A*N+B*N^2+C*N^3+...+D*N^(M+1)
!' 1
!'1^1*A+1^2*B+1^3*C=Σ(N^2)=1^2
!' N=1
!' 2
!'2^1*A+2^2*B+2^3*C=Σ(N^2)=1^2+2^2
!' N=1
!' 3
!'3^1*A+3^2*B+3^3*C=Σ(N^2)=1^2+2^2+3^2
!' N=1
DIM A(N+1,N+1),B(N+1)
FOR I=1 TO N+1
FOR J=1 TO N+1
LET A(I,J)=I^J
NEXT J
LET B(I)=SUMN(N,I)
NEXT I
MAT A=INV(A)
MAT X=A*B
END SUB

EXTERNAL FUNCTION SUMN(M,X)
FOR I=1 TO X
LET S=S+I^M
NEXT I
LET SUMN=S
END FUNCTION

EXTERNAL SUB DISPLAY(A(),N)
IF N>1 THEN
IF A(N)<0 THEN PRINT "-";
IF ABS(A(N))<>1 THEN
PRINT STR$(ABS(A(N)));"*X^";STR$(N);
ELSE
PRINT "X^";STR$(N);
END IF
END IF
FOR I=N-1 TO 2 STEP -1
IF A(I)<>0 THEN
IF A(I)<0 THEN PRINT "-"; ELSE PRINT "+";
IF ABS(A(I))<>1 THEN
PRINT STR$(ABS(A(I)));"*X^";STR$(I);
ELSEIF ABS(A(I))=1 THEN
PRINT "X^";STR$(I);
END IF
END IF
NEXT I
IF A(1)<>0 THEN
IF N>1 THEN
IF A(1)<0 THEN PRINT "-"; ELSE PRINT "+";
END IF
IF ABS(A(1))<>1 THEN
PRINT STR$(ABS(A(1)));"*X";
ELSEIF ABS(A(1))=1 THEN
PRINT "X";
END IF
END IF
PRINT
END SUB

!ベルヌーイ多項式からΣi^kを求める

OPTION ARITHMETIC rational

OPTION BASE 0

!ベルヌーイ数 Bk　※２項係数の２重級数で表現する
!　 k　　　　m
!Bk=Σ1/(m+1)Σ(-1)^j*mCj*j^k
!　 m=0　　　j=0

!　 k　　　　　 k
!Bk=Σ(-1)^j*j^kΣmCj/(m+1)
!　 j=0　　　　 m=j

FUNCTION B(k) !ベルヌーイ数
LET Bk=0
FOR m=0 TO k
LET s=0
FOR j=0 TO m
LET s=s+(-1)^j*comb(m,j)*j^k
NEXT j
LET Bk=Bk+s/(m+1)
NEXT m
LET B=Bk
END FUNCTION

!ベルヌーイ多項式 Bn(x)
!　　　n
!Bn(x)=ΣnCk*Bk*x^(n-k)　Bkはベルヌーイ数
!　　　k=0

SUB BernoulliB(n,Bn())
FOR k=0 TO n
LET Bn(n-k)=comb(n,k)*B(k)
NEXT k
END SUB

!多項式の計算
!１変数Xの多項式（N次）のCi*X^iの係数を、配列C(i)=Ciで表す。

!演算関連
FUNCTION poly_val(n,a(),v) !多項式を計算する（数値代入）　a(α)
LET k=a(n)
FOR i=n-1 TO 0 STEP -1 !Hornerの方法
LET k=k*v+a(i) !(…(((AnX+An-1)X+An-2)X+An-3)X+…+A1)X+A0
NEXT i
LET poly_val=k
END FUNCTION
SUB poly_integral(n,a(), s()) !積分
FOR i=n+1 TO 1 STEP -1
LET s(i)=a(i-1)/i
NEXT i
LET s(0)=0
END SUB
!------------------------------ ここまでがサブルーチン

!1^k+2^k+…+n^k=Σi^k=Wk(n)とすると
!Wk(n)=k∫Wk-1(x)dx + Bk(1)*n　ただし、Bk(x)はベルヌーイ多項式

!W0(n)=0*∫W?(x)dx + B0(1)*n=n
!W1(n)=1*∫W0(x)dx + B1(1)*n=1/2*n^2 + 1/2*n
!W2(n)=2*∫W1(x)dx + B2(1)*n=2*(1/6*n^3+1/4*n^2) + 1/6*n
!W3(n)=3*∫W2(x)dx + B3(1)*n=3*(1/12*n^4+1/6*n^3+1/12*n^2) + 0*n
!　:

LET maxlevel=20

DIM n(maxlevel+1) !nの定義
LET n(1)=1

DIM w(maxlevel+1) !wの定義

FOR p=0 TO maxlevel
CALL poly_integral(maxlevel,w, w) !p∫W0(x)dx
MAT w=p*w

DIM Bn(maxlevel) !ベルヌーイ多項式の係数
CALL BernoulliB(p,Bn)
DIM T1(maxlevel+1)
MAT T1=poly_val(maxlevel,Bn,1)*n !Bk(1)*n

MAT w=w+T1

PRINT "Σi^";p
MAT PRINT w !係数のみ表示する　c0+c1*n+c2*n^2+ … + cn*n^n
NEXT p

END

!ベルヌーイ数からΣi^kを求める

OPTION ARITHMETIC rational

OPTION BASE 0

!ベルヌーイ数 Bk　※２項係数の２重級数で表現する
!　 k　　　　m
!Bk=Σ1/(m+1)Σ(-1)^j*mCj*j^k
!　 m=0　　　j=0

!　 k　　　　　 k
!Bk=Σ(-1)^j*j^kΣmCj/(m+1)
!　 j=0　　　　 m=j

FUNCTION B(k) !ベルヌーイ数
LET Bk=0
FOR m=0 TO k
LET s=0
FOR j=0 TO m
LET s=s+(-1)^j*comb(m,j)*j^k
NEXT j
LET Bk=Bk+s/(m+1)
NEXT m
LET B=Bk
END FUNCTION
!------------------------------ ここまでがサブルーチン

LET maxlevel=20

DIM w(maxlevel+1) !wの定義

!n-1　　　　　m
!Σj^m=1/(m+1)Σm+1Ck*Bk*n^(m+1-k)　ただし、Bkはベルヌーイ数
!j=0　　　　　k=0

FOR p=0 TO maxlevel
MAT w=ZER
PRINT "Σ[i=0,n] i^";p;"="
FOR k=0 TO p
LET w(p+1-k)=comb(p+1,k)*B(k)/(p+1) !Σ[j=0,n-1]
NEXT k

LET w(p)=w(p)+1 !+n^mで補正して、Σ[j=0,n]

MAT PRINT w !係数のみ表示する　c0+c1*n+c2*n^2+ … + cn*n^n
NEXT p

END

!ライプニッツの方法からΣi^kを求める

!参考 http://www004.upp.so-net.ne.jp/s_honma/sum/sum.htm

!an=nC0*a0 + nC1*b1 + nC2*c2 + nC3*d3 + …

!Σi^2の場合
!(0) 1 5 14 30 55 91 140 204 285 385 …　a0,a1,a2,a3,…　ak=Σi^k
! 0 (1) 4 9 16 25 36 49 64 81 100 …　b0,b1,b2,b3,…　階差数列
! 0 0 (3) 5 7 9 11 13 15 17 19 …　c0,c1,c2,c3,…　階差数列
! 0 0 0 (2) 2 2 2 2 2 2 2 …　d0,d1,d2,d3,…　階差数列
! 0 0 0 0 (0) 0 0 0 0 0 0 …　e0,e1,e2,e3,…　階差数列
!　 :
!　 :
!Sn=nC0*0+nC1*1+nC2*3+nC3*2　※nC4*e4以降は０
!　=1*0+n*1+n(n-1)/2!*3+n(n-1)(n-2)/3!*2

OPTION ARITHMETIC rational !有理数で計算する

OPTION BASE 0

LET maxlevel=30 !最高次数

!※多項式の計算
!１変数Xの多項式（N次）のCi*X^iの係数を、配列C(i)=Ciで表す。
SUB poly_mul(n,a(),b(), x()) !乗算x=a*b
DIM xx(100)
MAT xx=ZER(2*n)
FOR i=n TO 0 STEP -1
FOR j=n TO 0 STEP -1
LET xx(i+j)=xx(i+j)+a(i)*b(j) !分配する
NEXT j
NEXT i
MAT x=xx !copy it
END SUB

!その他
DIM c1(100)
MAT c1=ZER(maxlevel)
LET c1(0)=1 !定数１

SUB xxx(n,k, w()) !変数nの多項式 n*(n-1)*(n-2)*…*(n-k+1)/k!　※nCk
DIM T1(100)
MAT T1=ZER(n)
LET T1(1)=1 !n-0

MAT w=ZER
LET w(0)=1 !1

IF k>0 THEN
FOR p=1 TO k
CALL poly_mul(n,w,T1, w) !1*(n-0)*(n-1)*…*(n-k+1)
MAT T1=T1-c1 !n-p
NEXT p

MAT w=(1/fact(k))*w !/k!
ELSEIF k=0 THEN
EXIT SUB
ELSE
PRINT "未定義です。";k
STOP
END IF
END SUB
!------------------------------ ここまでがサブルーチン

LET k=2 !Σ[i=1,n]i^k

DIM a(maxlevel,maxlevel)
LET s=0
FOR i=1 TO maxlevel !a0=0
LET s=s+i^k !ai
LET a(0,i)=s
NEXT i
FOR i=1 TO maxlevel !階差数列 bn,cn,dn,…
FOR J=i TO maxlevel
LET a(i,J)=a(i-1,J)-a(i-1,J-1)
NEXT J
NEXT i

MAT PRINT a

DIM w(maxlevel*2)
CALL xxx(maxlevel,q, w) !nC0
MAT w=a(0,0)*w !*a0

LET q=1
DO WHILE a(q,q)>0 !対角線の成分
DIM T2(maxlevel*2)
CALL xxx(maxlevel,q, T2) !nCk

MAT T2=a(q,q)*T2 !*b1,*c2,*d3,…
MAT w=w+T2

LET q=q+1 !次へ
LOOP

MAT PRINT w !係数のみ表示する　c0+c1*n+c2*n^2+ … + cn*n^n

END

和の公式	返事を書くノートメニュー
しばっち <dihjvcfsyu> 2008/02/10 17:36:38
ΣX^Nを求める PUBLIC NUMERIC MAXLEVEL LET MAXLEVEL=20 DIM Y(MAXLEVEL+1),X(MAXLEVEL+1) FOR I=1 TO MAXLEVEL PRINT "ΣX^";STR$(I);"="; CALL SUM(X,I) CALL DISPLAY(X,I+1) NEXT I END EXTERNAL SUB SUM(X(),N) !'Σ(X^M)=AN+BN^2+CN^3+...+DN^(M+1) !' 1 !'1^1A+1^2B+1^3C=Σ(N^2)=1^2 !' N=1 !' 2 !'2^1A+2^2B+2^3C=Σ(N^2)=1^2+2^2 !' N=1 !' 3 !'3^1A+3^2B+3^3C=Σ(N^2)=1^2+2^2+3^2 !' N=1 DIM A(N+1,N+1),B(N+1) FOR I=1 TO N+1 FOR J=1 TO N+1 LET A(I,J)=I^J NEXT J LET B(I)=SUMN(N,I) NEXT I MAT A=INV(A) MAT X=AB END SUB EXTERNAL FUNCTION SUMN(M,X) FOR I=1 TO X LET S=S+I^M NEXT I LET SUMN=S END FUNCTION EXTERNAL SUB DISPLAY(A(),N) IF N>1 THEN IF A(N)<0 THEN PRINT "-"; IF ABS(A(N))<>1 THEN PRINT STR$(ABS(A(N)));"X^";STR$(N); ELSE PRINT "X^";STR$(N); END IF END IF FOR I=N-1 TO 2 STEP -1 IF A(I)<>0 THEN IF A(I)<0 THEN PRINT "-"; ELSE PRINT "+"; IF ABS(A(I))<>1 THEN PRINT STR$(ABS(A(I)));"X^";STR$(I); ELSEIF ABS(A(I))=1 THEN PRINT "X^";STR$(I); END IF END IF NEXT I IF A(1)<>0 THEN IF N>1 THEN IF A(1)<0 THEN PRINT "-"; ELSE PRINT "+"; END IF IF ABS(A(1))<>1 THEN PRINT STR$(ABS(A(1)));"*X"; ELSEIF ABS(A(1))=1 THEN PRINT "X"; END IF END IF PRINT END SUB

Re: 和の公式	返事を書くノートメニュー
山中和義 <drdlxujciw> 2008/02/12 15:00:46
!ベルヌーイ多項式からΣi^kを求める OPTION ARITHMETIC rational OPTION BASE 0 !ベルヌーイ数 Bk　※２項係数の２重級数で表現する !　 k　　　　m !Bk=Σ1/(m+1)Σ(-1)^jmCjj^k !　 m=0　　　j=0 !　 k　　　　　 k !Bk=Σ(-1)^jj^kΣmCj/(m+1) !　 j=0　　　　 m=j FUNCTION B(k) !ベルヌーイ数 LET Bk=0 FOR m=0 TO k LET s=0 FOR j=0 TO m LET s=s+(-1)^jcomb(m,j)j^k NEXT j LET Bk=Bk+s/(m+1) NEXT m LET B=Bk END FUNCTION !ベルヌーイ多項式 Bn(x) !　　　n !Bn(x)=ΣnCkBkx^(n-k)　Bkはベルヌーイ数 !　　　k=0 SUB BernoulliB(n,Bn()) FOR k=0 TO n LET Bn(n-k)=comb(n,k)B(k) NEXT k END SUB !多項式の計算 !１変数Xの多項式（N次）のCiX^iの係数を、配列C(i)=Ciで表す。 !演算関連 FUNCTION poly_val(n,a(),v) !多項式を計算する（数値代入）　a(α) LET k=a(n) FOR i=n-1 TO 0 STEP -1 !Hornerの方法 LET k=kv+a(i) !(…(((AnX+An-1)X+An-2)X+An-3)X+…+A1)X+A0 NEXT i LET poly_val=k END FUNCTION SUB poly_integral(n,a(), s()) !積分 FOR i=n+1 TO 1 STEP -1 LET s(i)=a(i-1)/i NEXT i LET s(0)=0 END SUB !------------------------------ ここまでがサブルーチン !1^k+2^k+…+n^k=Σi^k=Wk(n)とすると !Wk(n)=k∫Wk-1(x)dx + Bk(1)n　ただし、Bk(x)はベルヌーイ多項式 !W0(n)=0∫W?(x)dx + B0(1)n=n !W1(n)=1∫W0(x)dx + B1(1)n=1/2n^2 + 1/2n !W2(n)=2∫W1(x)dx + B2(1)n=2(1/6n^3+1/4n^2) + 1/6n !W3(n)=3∫W2(x)dx + B3(1)n=3(1/12n^4+1/6n^3+1/12n^2) + 0n !　: LET maxlevel=20 DIM n(maxlevel+1) !nの定義 LET n(1)=1 DIM w(maxlevel+1) !wの定義 FOR p=0 TO maxlevel CALL poly_integral(maxlevel,w, w) !p∫W0(x)dx MAT w=pw DIM Bn(maxlevel) !ベルヌーイ多項式の係数 CALL BernoulliB(p,Bn) DIM T1(maxlevel+1) MAT T1=poly_val(maxlevel,Bn,1)n !Bk(1)n MAT w=w+T1 PRINT "Σi^";p MAT PRINT w !係数のみ表示する　c0+c1n+c2n^2+ … + cnn^n NEXT p END

Re: !ベルヌーイ多項式からΣi^kを求める	返事を書くノートメニュー
山中和義 <drdlxujciw> 2008/02/12 16:26:21
!ベルヌーイ数からΣi^kを求める OPTION ARITHMETIC rational OPTION BASE 0 !ベルヌーイ数 Bk　※２項係数の２重級数で表現する !　 k　　　　m !Bk=Σ1/(m+1)Σ(-1)^jmCjj^k !　 m=0　　　j=0 !　 k　　　　　 k !Bk=Σ(-1)^jj^kΣmCj/(m+1) !　 j=0　　　　 m=j FUNCTION B(k) !ベルヌーイ数 LET Bk=0 FOR m=0 TO k LET s=0 FOR j=0 TO m LET s=s+(-1)^jcomb(m,j)j^k NEXT j LET Bk=Bk+s/(m+1) NEXT m LET B=Bk END FUNCTION !------------------------------ ここまでがサブルーチン LET maxlevel=20 DIM w(maxlevel+1) !wの定義 !n-1　　　　　m !Σj^m=1/(m+1)Σm+1CkBkn^(m+1-k)　ただし、Bkはベルヌーイ数 !j=0　　　　　k=0 FOR p=0 TO maxlevel MAT w=ZER PRINT "Σ[i=0,n] i^";p;"=" FOR k=0 TO p LET w(p+1-k)=comb(p+1,k)B(k)/(p+1) !Σ[j=0,n-1] NEXT k LET w(p)=w(p)+1 !+n^mで補正して、Σ[j=0,n] MAT PRINT w !係数のみ表示する　c0+c1n+c2n^2+ … + cn*n^n NEXT p END

Re: !ベルヌーイ数からΣi^kを求める	返事を書くノートメニュー
山中和義 <drdlxujciw> 2008/02/12 22:29:13 ** この記事は1回修正されてます
!ライプニッツの方法からΣi^kを求める !参考 http://www004.upp.so-net.ne.jp/s_honma/sum/sum.htm !an=nC0a0 + nC1b1 + nC2c2 + nC3d3 + … !Σi^2の場合 !(0) 1 5 14 30 55 91 140 204 285 385 …　a0,a1,a2,a3,…　ak=Σi^k ! 0 (1) 4 9 16 25 36 49 64 81 100 …　b0,b1,b2,b3,…　階差数列 ! 0 0 (3) 5 7 9 11 13 15 17 19 …　c0,c1,c2,c3,…　階差数列 ! 0 0 0 (2) 2 2 2 2 2 2 2 …　d0,d1,d2,d3,…　階差数列 ! 0 0 0 0 (0) 0 0 0 0 0 0 …　e0,e1,e2,e3,…　階差数列 !　 : !　 : !Sn=nC00+nC11+nC23+nC32　※nC4e4以降は０ !　=10+n1+n(n-1)/2!3+n(n-1)(n-2)/3!2 OPTION ARITHMETIC rational !有理数で計算する OPTION BASE 0 LET maxlevel=30 !最高次数 !※多項式の計算 !１変数Xの多項式（N次）のCiX^iの係数を、配列C(i)=Ciで表す。 SUB poly_mul(n,a(),b(), x()) !乗算x=ab DIM xx(100) MAT xx=ZER(2n) FOR i=n TO 0 STEP -1 FOR j=n TO 0 STEP -1 LET xx(i+j)=xx(i+j)+a(i)b(j) !分配する NEXT j NEXT i MAT x=xx !copy it END SUB !その他 DIM c1(100) MAT c1=ZER(maxlevel) LET c1(0)=1 !定数１ SUB xxx(n,k, w()) !変数nの多項式 n(n-1)(n-2)…(n-k+1)/k!　※nCk DIM T1(100) MAT T1=ZER(n) LET T1(1)=1 !n-0 MAT w=ZER LET w(0)=1 !1 IF k>0 THEN FOR p=1 TO k CALL poly_mul(n,w,T1, w) !1(n-0)(n-1)…(n-k+1) MAT T1=T1-c1 !n-p NEXT p MAT w=(1/fact(k))w !/k! ELSEIF k=0 THEN EXIT SUB ELSE PRINT "未定義です。";k STOP END IF END SUB !------------------------------ ここまでがサブルーチン LET k=2 !Σ[i=1,n]i^k DIM a(maxlevel,maxlevel) LET s=0 FOR i=1 TO maxlevel !a0=0 LET s=s+i^k !ai LET a(0,i)=s NEXT i FOR i=1 TO maxlevel !階差数列 bn,cn,dn,… FOR J=i TO maxlevel LET a(i,J)=a(i-1,J)-a(i-1,J-1) NEXT J NEXT i MAT PRINT a DIM w(maxlevel2) CALL xxx(maxlevel,q, w) !nC0 MAT w=a(0,0)w !a0 LET q=1 DO WHILE a(q,q)>0 !対角線の成分 DIM T2(maxlevel2) CALL xxx(maxlevel,q, T2) !nCk MAT T2=a(q,q)T2 !b1,c2,d3,… MAT w=w+T2 LET q=q+1 !次へ LOOP MAT PRINT w !係数のみ表示する　c0+c1n+c2n^2+ … + cn*n^n END