ベルヌーイ数の数値計算

これがベルヌーイ数を求める最速のアルゴリズムだと思います。
直接代入。これに限ります。

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

B(256)まで増やしたのはこちら

一応逐次的に求める事が出来るアルゴリズムがありまして実装しましたが、全然良くないです[1]。
4倍精度演算を使って\(B_{30}\)を求めようとするとたったの8桁しか一致しません。
倍精度演算のみだと\(B_{14}\)を求めようとするとたったの6桁しか一致しません。

参考文献

[1]Algorithmic description -Bernoulli number, wikipedia

スポンサーリンク


コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です