fortran66のブログ

fortran について書きます。

分割数 p(200) 直接計算

George E. Andrews, The Theory of Partitions (Encyclopedia of Mathematics and its Applications) の中に、Hardy-Ramanujan-Rademacher の p(n) の式の p(200) に対する具体的計算例が載っています。

amazon AA

Chap.5.1 p.77
f:id:fortran66:20170819014720p:plain


以前、計算した時のプログラムで再現してみたのですが、第三項が少し違うという結果が・・・

     1             3972998993185.896                        -0.000
     2                     36282.978                        -0.000
     3                       -87.584                         0.000
     4                         5.147                         0.000
     5                         1.424                         0.000
     6                         0.071                         0.000
     7                        -0.000                        -0.000
     8                         0.044                         0.000

p(200)             3972999029387.976                        -0.000
続行するには何かキーを押してください . . .

原因は今のところ不明です。
無限級数を第 50 項まで計算すると、整数化するので間違っているわけでもなさそうな・・・?
はぁ、やな感じ~。

     1             3972998993185.896                        -0.000
     2                     36282.978                        -0.000
     3                       -87.584                         0.000
     4                         5.147                         0.000
     5                         1.424                         0.000
     6                         0.071                         0.000
     7                        -0.000                        -0.000
     8                         0.044                         0.000
     9                         0.026                         0.000
    10                         0.009                         0.000
    11                         0.000                         0.000
    12                        -0.004                        -0.000
    13                        -0.000                        -0.000
    14                         0.000                         0.000
    15                        -0.003                        -0.000
    16                        -0.004                        -0.000
    17                         0.000                        -0.000
    18                        -0.002                         0.000
    19                        -0.000                        -0.000
    20                         0.001                        -0.000
    21                        -0.000                         0.000
    22                        -0.000                        -0.000
    23                         0.002                         0.000
    24                        -0.001                         0.000
    25                         0.002                        -0.000
    26                        -0.000                         0.000
    27                        -0.001                        -0.000
    28                         0.000                         0.000
    29                        -0.000                        -0.000
    30                        -0.001                        -0.000
    31                        -0.000                         0.000
    32                         0.000                        -0.000
    33                        -0.000                        -0.000
    34                         0.000                        -0.000
    35                         0.000                        -0.000
    36                         0.000                         0.000
    37                         0.001                        -0.000
    38                        -0.000                         0.000
    39                        -0.000                        -0.000
    40                        -0.000                        -0.000
    41                         0.001                         0.000
    42                        -0.000                         0.000
    43                        -0.001                        -0.000
    44                         0.000                        -0.000
    45                        -0.000                         0.000
    46                        -0.000                        -0.000
    47                        -0.001                        -0.000
    48                        -0.000                        -0.000
    49                        -0.000                        -0.000
    50                         0.000                         0.000

p(200)             3972999029388.000                        -0.000
続行するには何かキーを押してください . . .

fortran66.hatenablog.com

プログラム

以前のとちょっとだけ変えました(一部 Andrews の定義に書き直し)が結果は変わってません。

無限級数を第 8 項で打ち切っています。

    module m_part
      implicit none
      integer , parameter :: kd = kind(1.0q0)
      real(kd), parameter :: pi = 4 * atan(1.0_kd)
    contains  
      complex(kd) function cpartition(n)
        integer, intent(in) :: n
        integer, parameter :: kmax = 8
        real(kd) :: f, s, t, u
        complex(kd) :: c
        integer :: j, m, k
        cpartition = (0.0_kd, 0.0_kd)
        do k = 1, kmax
          f = sqrt(k / 2.0_kd) / pi 
          c = (0.0_kd, 0.0_kd)
          do m = 1, k
            if (igcd(m, k) == 1) then
              s = -real(2 * n * m, kd) / k + s_hk(m, k)
              c = c + exp(cmplx(0.0_kd, pi * s, kind = kd)) 
            end if
          end do
          cpartition = cpartition + f * c * da(k, n) 
          print '(i6, 2f30.3)', k, f * c * da(k, n) 
        end do  
      contains
        real(kd) function s_hk(h, k)
          integer, intent(in) :: h, k
          real(kd) :: t, u
          integer :: j
          s_hk = 0.0_kd
          do j = 1, k - 1
             t = real(m * j, kd) / k
             u = real(    j, kd) / k
             s_hk = s_hk + (u - floor(u) - 0.5_kd) * (t - floor(t) - 0.5_kd)                
          end do
        end function s_hk
      end function cpartition
  
      integer pure function igcd(ia, ib)
        integer, intent(in) :: ia, ib
        integer :: m(2)
        m = [ia, ib]
        do 
          m = [m(2), mod(m(1), m(2))]
          if (m(2) == 0) exit
        end do
        igcd = m(1)
      end function igcd  
      
      real(kd) pure function da(k, n) ! d/dx sinh( f * xs ) / xs ; xs = sqrt(x - 1/24)
        integer, intent(in) :: k, n 
        real(kd) :: xs, f
        xs = sqrt( real(n, kd) - 1.0_kd / 24.0_kd )
        f  = pi / k * sqrt(2.0_kd / 3.0_kd)
        da = ( f * cosh(f * xs) - sinh(f * xs) / xs ) / (2 * xs * xs) 
      end function da
    end module m_part
  
    program Pn
      use m_part
      implicit none
      integer :: i
      complex(kind = kd) :: z
      do i = 200, 200 
          z = cpartition(i)
        print *
        print '(a, i3, a, 2f30.3)', 'p(', i, ')', z
      end do  
    end program Pn