分割数

Hardy-Ramanujan-Rademacher の公式の exp の複素数和の所を、どうせ虚部は x 軸に対して対称の位置に出て打ち消しあうので、cos に置き換えて実数のみで計算します。

Fortran と Julia の以前のプログラムを改訂します。

Hardy-Ramanujan-Rademacher の公式

$\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} }$

無限和の内必要な項の見積もりの式があったのですが、失念しました。仕方ないので、適当に大きめに入れています。

Introduction to Analytic Number Theory (Undergraduate Texts in Mathematics)

作者:Apostol, Tom M.
発売日: 2010/12/04
メディア: ペーパーバック

Modular Functions and Dirichlet Series in Number Theory (Graduate Texts in Mathematics) (Graduate Texts in Mathematics, 41)

作者:Apostol, Tom M.
発売日: 2012/12/21
メディア: ペーパーバック

Integer Partitions

作者:Andrews, George
発売日: 2010/07/12
メディア: ペーパーバック

The Theory of Partitions (Encyclopedia of Mathematics and its Applications, Series Number 2)

作者:Andrews, George E.
発売日: 2008/01/12
メディア: ペーパーバック

整数の分割

作者:ジョージ・アンドリュース,キムモ・エリクソン
発売日: 2006/05/01
メディア: 単行本

Fortran 版

P(721) を計算し、無限和の各項の大きさの変化を見ます。

fortran66.hatenablog.com

    module partition_m
        use, intrinsic :: iso_fortran_env
        implicit none
        integer , parameter :: kd = real128
        real(kd), parameter :: pi = 4 * atan(1.0_kd)
    contains  

        impure elemental real(kd) function partition(n)
            integer, intent(in) :: n
            integer, parameter :: kmax = 21
            integer :: k
            partition = 0.0_kd
            do k = 1, kmax
                partition = partition + sqrt(k / 2.0_kd) / pi * fa(n, k) * fs(k, n) 
                print '(i6, 2f40.5)', k, sqrt(k / 2.0_kd) / pi * fa(n, k) * fs(k, n) 
            end do         
        end function partition
    
        real(kd) pure function fs(k, n) ! d/dx sinh( f * xs ) / xs ; xs = sqrt(x - 1/24)
            integer, intent(in) :: k, n 
            real(kd) :: xs, f
            xs = sqrt( n - 1 / 24.0_kd )
            f  = pi / k * sqrt(2 / 3.0_kd)
            fs = ( f * cosh(f * xs) - sinh(f * xs) / xs ) / (2 * xs * xs) 
        end function fs
        
        pure real(kd) function fa(n, k)
            integer, intent(in) :: n, k
            real(kd) :: s
            integer :: m
            fa = 0.0_kd
            do m = 1, k
                if (igcd(m, k) == 1) then
                    s = -real(2 * n * m, kd) / k + fb(m, k)
                    fa = fa + cos(pi * s) 
                end if
            end do
        end function fa
        
        pure real(kd) function fb(h, k)
            integer, intent(in) :: h, k
            integer :: j
            fb = 0.0_kd
            do j = 1, k - 1
                fb = fb + j * (h * j / real(k, kd) - floor(h * j / real(k, kd)) - 0.5_kd) 
            end do
            fb = fb / k
        end function fb
  
        pure integer function igcd(ia, ib)
            integer, value :: ia, ib
            do 
                igcd = ib
                ib = mod(ia, ib)
                if (ib == 0) exit
                ia = igcd
            end do
        end function igcd  
    end module partition_m
  
    program Pn
        use partition_m
        implicit none
        integer :: i
        real(kd) :: r       
!        print *, nint(partition([1:10]))
        
        do i = 721, 721 
            r = partition(i)
            print *
            print '(a, i3, a, f40.5)', 'p(', i, ')', r
        end do  
    end program Pn

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

Julia 版

百までの分割数を計算します。

fortran66.hatenablog.com

function fs(k, z)
    xs = sqrt(z-1/24)
    fs = √(2/3) * π / k  
    (fs * cosh(fs * xs) - sinh(fs * xs) / xs) / 2xs^2
end

function fb(k, h)
    s = 0
    for j = 1:k-1
        s += j * (h * j / k - floor(h * j / k) - 1/2)
    end 
    s / k
end

function fa(n, k)
    s = 0
    for h = 1:k 
        if (gcd(h,k)==1) 
           z = -2n*h/k + fb(k, h)
           s += cos(π*z)
        end    
    end 
    s
end

function p(n)
    s = 0
    for k=1:10
        s += sqrt(k) * fa(n, k) * fs(k, n)
    end 
    s / √2pi  
end

using Printf
for i = 1:100
    pn = p(i)
    @printf( "%10d %20d %20.7f \n", i, round(Int64, real(pn)), real(pn)) 
end

         1                    1            1.0063342 
         2                    2            2.0044611 
         3                    3            3.0023039 
         4                    5            4.9960379 
         5                    7            6.9972722 
         6                   11           10.9886684 
         7                   15           14.9997219 
         8                   22           21.9944150 
         9                   30           30.0038956 
        10                   42           42.0078618 
        11                   56           56.0028080 
        12                   77           76.9987022 
        13                  101          101.0035727 
        14                  135          135.0035003 
        15                  176          175.9936875 
        16                  231          231.0001458 
        17                  297          296.9912596 
        18                  385          385.0031500 
        19                  490          489.9937920 
        20                  627          627.0072493 
        21                  792          792.0080395 
        22                 1002         1001.9931672 
        23                 1255         1254.9993226 
        24                 1575         1575.0069153 
        25                 1958         1958.0040403 
        26                 2436         2435.9892205 
        27                 3010         3010.0028149 
        28                 3718         3717.9938783 
        29                 4565         4564.9961219 
        30                 5604         5604.0040496 
        31                 6842         6842.0039095 
        32                 8349         8349.0100898 
        33                10143        10142.9938646 
        34                12310        12309.9901594 
        35                14883        14883.0072610 
        36                17977        17977.0032812 
        37                21637        21637.0036938 
        38                26015        26014.9961219 
        39                31185        31184.9842031 
        40                37338        37338.0070623 
        41                44583        44582.9966134 
        42                53174        53174.0087179 
        43                63261        63261.0087796 
        44                75175        75174.9936763 
        45                89134        89134.0016629 
        46               105558       105557.9980458 
        47               124754       124753.9977149 
        48               147273       147272.9904516 
        49               173525       173525.0070202 
        50               204226       204226.0036898 
        51               239943       239943.0021134 
        52               281589       281588.9910221 
        53               329931       329931.0021026 
        54               386155       386155.0101807 
        55               451276       451275.9887966 
        56               526823       526823.0079577 
        57               614154       614153.9990722 
        58               715220       715220.0009011 
        59               831820       831820.0001180 
        60               966467       966466.9958295 
        61              1121505      1121504.9930073 
        62              1300156      1300155.9942735 
        63              1505499      1505499.0088466 
        64              1741630      1741630.0099512 
        65              2012558      2012558.0033675 
        66              2323520      2323519.9945812 
        67              2679689      2679689.0022615 
        68              3087735      3087734.9963524 
        69              3554345      3554344.9974198 
        70              4087968      4087967.9953353 
        71              4697205      4697204.9976341 
        72              5392783      5392782.9993796 
        73              6185689      6185689.0052464 
        74              7089500      7089500.0034465 
        75              8118264      8118264.0027078 
        76              9289091      9289091.0111846 
        77             10619863     10619862.9921330 
        78             12132164     12132163.9936887 
        79             13848650     13848649.9941123 
        80             15796476     15796476.0006777 
        81             18004327     18004326.9991289 
        82             20506255     20506255.0020917 
        83             23338469     23338469.0030397 
        84             26543660     26543659.9989915 
        85             30167357     30167357.0013871 
        86             34262962     34262962.0040397 
        87             38887673     38887673.0087649 
        88             44108109     44108108.9878901 
        89             49995925     49995924.9983584 
        90             56634173     56634173.0039072 
        91             64112359     64112358.9928125 
        92             72533807     72533806.9986574 
        93             82010177     82010177.0021016 
        94             92669720     92669719.9965452 
        95            104651419    104651419.0057203 
        96            118114304    118114304.0000427 
        97            133230930    133230930.0105914 
        98            150198136    150198136.0060634 
        99            169229875    169229874.9836169 
       100            190569292    190569292.0010892

fortran66のブログ

fortran について書きます。

【メモ帳】分割数計算　改訂

分割数

Hardy-Ramanujan-Rademacher の公式

Fortran 版

Julia 版