George E. Andrews, The Theory of Partitions (Encyclopedia of Mathematics and its Applications) の中に、Hardy-Ramanujan-Rademacher の p(n) の式の p(200) に対する具体的計算例が載っています。
追記H29/12/29 結果のズレに関して続編あり
分割数 p(200) 直接計算 (その２） - fortran66のブログ
amazon AA
Chap.5.1 p.70
以前、計算した時のプログラムで再現してみたのですが、第三項が少し違うという結果が･･･
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
続行するには何かキーを押してください . . .
プログラム
以前のとちょっとだけ変えました（一部 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
