分割数 p(721)

ハーディーの「ラマヌジャンその生涯と業績に想起された主題による十二の講義」に 721 の分割数計算の級数和が出ていたので（講義VIII, p.193）、検算のため比較してみました。

ラマヌジャンその生涯と業績に想起された主題による十二の講義 (数学クラシックス)

作者:G.H. ハーディ
発売日: 2016/09/10
メディア: 単行本

Ramanujan: Twelve Lectures on Subjects Suggested by His Life and Work (Ams Chelsea Publishing)

作者:Hardy, G. H.
発売日: 1999/11/25
メディア: ハードカバー

f:id:fortran66:20190810183240p:plain — 分割数 p(721)

実行結果

多分、本の方は総和の小数点以下の０の数が１個少ないのではないかと思います。

     1       161061755750279601828302117.84821                                -0.00000
     2                     -124192062781.96844                                 0.00000
     3                           -706763.61926                                 0.00000
     4                              2169.16830                                -0.00000
     5                                 0.00000                                 0.00000
     6                                14.20704                                -0.00000
     7                                 6.07837                                -0.00000
     8                                 0.18926                                 0.00000
     9                                 0.04914                                -0.00000
    10                                 0.00000                                 0.00000
    11                                 0.08814                                 0.00000
    12                                -0.03525                                -0.00000
    13                                 0.03248                                 0.00000
    14                                -0.00687                                 0.00000
    15                                 0.00000                                -0.00000
    16                                -0.01133                                -0.00000
    17                                 0.00000                                -0.00000
    18                                -0.00554                                 0.00000
    19                                -0.00860                                -0.00000
    20                                -0.00000                                -0.00000
    21                                -0.00526                                 0.00000
 
p(721)       161061755750279477635534762.00040                                -0.00000

Hardy-Ramanujan-Rademacher の公式

無限和を 21 で止めています。

$\displaystyle{ p(n)=\sum_{k=1}^\infty{\sqrt{k}\over\pi\sqrt{2}}\sum_{0\lt h\leq k\\ (h,k)=1}e^{{-2\pi i n h\over k}+\pi i\sum_{j=1}^{k-1}{j\over k}({h j\over k}-\lfloor{h j\over k}\rfloor-{1\over2})} \left. {d\over dz}{\sinh\left({\pi\over k}\sqrt{{2\over3}(z-{1\over24})}\right)\over \sqrt{z-{1\over24}} }\right|_{z=n} }$

プログラム

プログラムは、以前 Fortran で書いたものを用いることにします。倍精度では p(250) 以降で桁が足りなくなるので、４倍精度で計算することにします。フォーマット文や、無限級数の打ち切りなどを調節しました。

fortran66.hatenablog.com

    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 = 21
        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, 2f40.5)', 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 = 721, 721 
          z = cpartition(i)
        print *
        print '(a, i3, a, 2f40.5)', 'p(', i, ')', z
      end do  
    end program Pn

fortran66のブログ

fortran について書きます。

【メモ帳】分割数 p(721) チェック

分割数 p(721)

実行結果

Hardy-Ramanujan-Rademacher の公式

プログラム