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

!ベルヌーイ多項式からΣ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

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