これがベルヌーイ数を求める最速のアルゴリズムだと思います。
直接代入。これに限ります。
n=0~100までサポートしています。
実用的にはこの範囲で足りるんじゃないかなと思います。
値はwolfram alphaより。
Table[bernoulli b(i), {i,0,100}] with 17 digits
program main
implicit none
integer::i
double precision::B
external::B
do i=0,30
write(6,'(i5,f30.17)')i,B(i)
enddo
stop
end program main
function B(n)
! 1st bernoulli number B(n), 0 <= n <= 100
! from wolfram alpha
! " Table[bernoulli b(i), {i,0,100}] with 17 digits "
implicit none
integer,intent(in)::n
double precision::B
double precision::A(0:100)
if(n.lt.0.or.n.gt.100)then
write(6,*)"Not expected n at bernoulli number function. "
write(6,*)"program stop"
stop
endif
A(0:100)=(/1d0 &
, -0.50000000000000000d0, 0.16666666666666667d0, 0d0, -0.033333333333333333d0, 0d0 &
, 0.023809523809523810d0, 0d0, -0.033333333333333333d0, 0d0, 0.075757575757575758d0, 0d0 &
, -0.25311355311355311d0, 0d0, 1.1666666666666667d0, 0d0, -7.0921568627450980d0, 0d0 &
, 54.971177944862155d0, 0d0, -529.12424242424242d0, 0d0, 6192.1231884057971d0, 0d0 &
, -86580.253113553114d0, 0d0, 1.4255171666666667d+6, 0d0, -2.7298231067816092d+7, 0d0 &
, 6.0158087390064237d+8, 0d0, -1.5116315767092157d+10, 0d0, 4.2961464306116667d+11, 0d0 &
, -1.3711655205088333d+13, 0d0, 4.8833231897359317d+14, 0d0, -1.9296579341940068d+16, 0d0 &
, 8.4169304757368262d+17, 0d0, -4.0338071854059455d+19, 0d0, 2.1150748638081992d+21, 0d0 &
, -1.2086626522296526d+23, 0d0, 7.5008667460769644d+24, 0d0, -5.0387781014810689d+26, 0d0 &
, 3.6528776484818123d+28, 0d0, -2.8498769302450882d+30, 0d0, 2.3865427499683628d+32, 0d0 &
, -2.1399949257225334d+34, 0d0, 2.0500975723478098d+36, 0d0, -2.0938005911346378d+38, 0d0 &
, 2.2752696488463516d+40, 0d0, -2.6257710286239576d+42, 0d0, 3.2125082102718033d+44, 0d0 &
, -4.1598278166794711d+46, 0d0, 5.6920695482035280d+48, 0d0, -8.2183629419784576d+50, 0d0 &
, 1.2502904327166993d+53, 0d0, -2.0015583233248370d+55, 0d0, 3.3674982915364374d+57, 0d0 &
, -5.9470970503135448d+59, 0d0, 1.1011910323627978d+62, 0d0, -2.1355259545253501d+64, 0d0 &
, 4.3328896986641192d+66, 0d0, -9.1885528241669328d+68, 0d0, 2.0346896776329074d+71, 0d0 &
, -4.7003833958035731d+73, 0d0, 1.1318043445484249d+76, 0d0, -2.8382249570693707d+78/)
B=A(n)
return
end function B
implicit none
integer::i
double precision::B
external::B
do i=0,30
write(6,'(i5,f30.17)')i,B(i)
enddo
stop
end program main
function B(n)
! 1st bernoulli number B(n), 0 <= n <= 100
! from wolfram alpha
! " Table[bernoulli b(i), {i,0,100}] with 17 digits "
implicit none
integer,intent(in)::n
double precision::B
double precision::A(0:100)
if(n.lt.0.or.n.gt.100)then
write(6,*)"Not expected n at bernoulli number function. "
write(6,*)"program stop"
stop
endif
A(0:100)=(/1d0 &
, -0.50000000000000000d0, 0.16666666666666667d0, 0d0, -0.033333333333333333d0, 0d0 &
, 0.023809523809523810d0, 0d0, -0.033333333333333333d0, 0d0, 0.075757575757575758d0, 0d0 &
, -0.25311355311355311d0, 0d0, 1.1666666666666667d0, 0d0, -7.0921568627450980d0, 0d0 &
, 54.971177944862155d0, 0d0, -529.12424242424242d0, 0d0, 6192.1231884057971d0, 0d0 &
, -86580.253113553114d0, 0d0, 1.4255171666666667d+6, 0d0, -2.7298231067816092d+7, 0d0 &
, 6.0158087390064237d+8, 0d0, -1.5116315767092157d+10, 0d0, 4.2961464306116667d+11, 0d0 &
, -1.3711655205088333d+13, 0d0, 4.8833231897359317d+14, 0d0, -1.9296579341940068d+16, 0d0 &
, 8.4169304757368262d+17, 0d0, -4.0338071854059455d+19, 0d0, 2.1150748638081992d+21, 0d0 &
, -1.2086626522296526d+23, 0d0, 7.5008667460769644d+24, 0d0, -5.0387781014810689d+26, 0d0 &
, 3.6528776484818123d+28, 0d0, -2.8498769302450882d+30, 0d0, 2.3865427499683628d+32, 0d0 &
, -2.1399949257225334d+34, 0d0, 2.0500975723478098d+36, 0d0, -2.0938005911346378d+38, 0d0 &
, 2.2752696488463516d+40, 0d0, -2.6257710286239576d+42, 0d0, 3.2125082102718033d+44, 0d0 &
, -4.1598278166794711d+46, 0d0, 5.6920695482035280d+48, 0d0, -8.2183629419784576d+50, 0d0 &
, 1.2502904327166993d+53, 0d0, -2.0015583233248370d+55, 0d0, 3.3674982915364374d+57, 0d0 &
, -5.9470970503135448d+59, 0d0, 1.1011910323627978d+62, 0d0, -2.1355259545253501d+64, 0d0 &
, 4.3328896986641192d+66, 0d0, -9.1885528241669328d+68, 0d0, 2.0346896776329074d+71, 0d0 &
, -4.7003833958035731d+73, 0d0, 1.1318043445484249d+76, 0d0, -2.8382249570693707d+78/)
B=A(n)
return
end function B
B(256)まで増やしたのはこちら。
一応逐次的に求める事が出来るアルゴリズムがありまして実装しましたが、全然良くないです[1]。
4倍精度演算を使って\(B_{30}\)を求めようとするとたったの8桁しか一致しません。
倍精度演算のみだと\(B_{14}\)を求めようとするとたったの6桁しか一致しません。
参考文献
[1]Algorithmic description -Bernoulli number, wikipedia
[adsense2]
