分割数 p(n) の直接計算 - fortran66のブログ

Hardy & Ramanujan と Rademacher によって、分割数 p(n) を漸化式を使わずに、直接計算できる式が与えられています。これを使って分割数を求めてみたいと思います。なお式の証明には、かなり高度な知識が要求されるようなので、結果だけ頂戴します。

分割数を直接計算できる式は、整数の分割という素朴なものとは思えないような複雑な形をしているので、数学科の連中が我等純真な素人を驚かすのによく持ち出してきます。ところが、計算して確かめてみることが可能な形で、はっきりと与えられることが少なく、しかも無限級数になっているので、ありがたみも少なめです。

G.E.Andrews の The Theory of Partitions に計算例が載っていますが、式が複雑で計算する気が起きませんでした。

J.Remmel & A.Mendes の Counting with Symmetric Functions に計算してみる気が起きる程度に簡単にまとめられた式が載っていたので、ここでは実際に数値を入れて計算してみます。

f:id:fortran66:20160709154256p:plain

この式は、G.E.Andrews の与える式と微妙に違っていますが、別途計算して見たところ、結果は同じになるようなので、24乗根の取り方の違いなのかなんなのか（追記：デデキント和の別表現のようです。）、どっちでもよさそうな感じです。（一度計算できると、 Andrews の式も微修正で計算できました。）

準備

とりあえず、 $\sinh({\pi\over k}\sqrt{{2\over3}(z-{1\over24})}) \over \sqrt{z-{1\over24}}$ の微分を求めます。素直に微分すればいいのですが、プログラムすることを考えて、適当に変数置換してやると、 $f={\pi\over k}\sqrt{2\over3}$ , $x_s=\sqrt{z-{1\over24}}$ とおいて、 ${f \cosh(f x_s) - {\sinh(f x_s)\over x_s}\over 2x_s^2}$
と表わせます。

なお式中の h に関する総和のところは、実際は $0\lt h\le k$ で等号が抜けていると思われます。k=1の時に問題になります。

ソース・プログラム

本来は、変数 k に関する和は無限和ですが、10 までで止めておきます。これでも大体十分です。k=100 まで和をとった場合と比較して、n ～ 1000 までで、結果に 0.3 程度の誤差が出ますが、四捨五入すると関係ありません。

    module m_part
      implicit none
      integer , parameter :: kd = kind(1.0q0)
      real(kd), parameter :: pi = 4 * atan(1.0_kd)
    contains  
      complex(kd) pure function cpartition(n)
        integer, intent(in) :: n
        integer, parameter :: kmax = 10
        real(kd) :: f, s, t
        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 
              do j = 1, k - 1
                t = real(m * j, kd) / k
                s = s + real(j, kd) / k * (t - floor(t) - 0.5_kd)                
              end do
              c = c + exp(cmplx(0.0_kd, pi * s, kind = kd)) 
            end if
          end do
          cpartition = cpartition + f * c * da(k, n) 
        end do  
      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 = 1, 1000
        print '(i7, f45.3, i30)', i, real( cpartition(i), kd ), nint(real( cpartition(i), kd ), kind = 8)
      end do  
    end program Pn

実行結果

      1                                        1.006                             1
      2                                        2.004                             2
      3                                        3.002                             3
      4                                        4.996                             5
      5                                        6.997                             7
      6                                       10.989                            11
      7                                       15.000                            15
      8                                       21.994                            22
      9                                       30.004                            30
     10                                       42.008                            42

400 近傍で 8 バイト整数の範囲を越えます。

    400                      6727090051741041925.961           6727090051741041926
    401                      7154640222653942320.999           7154640222653942321
    402                      7608802843339879269.005           7608802843339879269
    403                      8091200276484465580.978           8091200276484465581
    404                      8603551759348655060.036           8603551759348655060
    405                      9147679068859117602.063           9147679068859117602
    406                      9725512513742021729.066          -9223372036854775808
    407                     10339097267123947240.969          -9223372036854775808
    408                     10990600063775926994.010          -9223372036854775808
    409                     11682316277192317779.980          -9223372036854775808

    990         16183531619906475296861224625026.766          -9223372036854775808
    991         16839773100833956878604913215476.879          -9223372036854775808
    992         17522282609145324707635966077021.986          -9223372036854775808
    993         18232098083140097717852712346115.225          -9223372036854775808
    994         18970297947002453464660671155989.793          -9223372036854775808
    995         19738002669751617842096992232436.119          -9223372036854775808
    996         20536376383452971700767593594020.746          -9223372036854775808
    997         21366628562913781584556907794729.258          -9223372036854775808
    998         22230015769169865076825741905554.660          -9223372036854775808
    999         23127843459154899464880444632250.039          -9223372036854775808
   1000         24061467864032622473692149727991.395          -9223372036854775808
続行するには何かキーを押してください . . .

漸化式によって求めた p(1000) の値は p(1000)=24061467864032622473692149727991 となっていて、四捨五入後の結果と一致します。