オイラーの五角数を用いた分割数の漸化式

$\displaystyle{ p(n) = \sum\{ (-1)^k p( k(3k-1) / 2 ) + (-1)^k p( k(3k+1) / 2 ) \} }$

但し $p(0) = 1$ とする。

Julia の BigInt の試用目的で昔の計算をやり直してみます。

fortran66.hatenablog.com fortran66.hatenablog.com fortran66.hatenablog.com

整数の分割

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

Fortran

まずは Fortran で書いてそれを移植する形で進めます。昔、書いたものを修正していくつか書いてみます。

FMLIB によるもの 1

https://dmsmith.lmu.build/FMLIB.html

main ルーチン中で途中経過を出力しながら 10⁶ の分割数を計算します。

    program test
        use, intrinsic :: iso_fortran_env
        use FMZM
        implicit none
        integer, parameter :: ki = int64, kd = real64
        integer, parameter :: np = 10**6  
        integer(ki) :: i, j, k  
        type(im) :: m, mp(0:np), isgn
        character (1200) :: str
        integer :: i0, i1, ic

        call system_clock(i0, ic)
        
        mp(0) = to_im(1)
        do i = 1, np
            mp(i) = to_fm(0)
        
            isgn = 1
            do j = 1, i
                k = i - j * ( 3 * j - 1) / 2
                if ( k < 0 ) exit
                isgn = -isgn
                mp(i) = mp(i) - isgn * mp(k) 
            end do
            
            isgn = 1
            do j = 1, i
                k = i - j * ( 3 * j + 1) / 2
                if ( k < 0 ) exit
                isgn = -isgn
                mp(i) = mp(i) - isgn * mp(k)
            end do
            
            call system_clock(i1)
           
            if (mod(i, 100000) == 0) then
                print '(f10.2)', real(i1 - i0) / ic
                call im_form('i1200', mp(i), str)
                write(*, '(a, i0, a, a)') 'p(', i, ')=', trim(adjustl(str))  
            end if
        end do
    end program test

百万の分割数を求めるのに、1,300秒ほどかかっています。

     21.19
p(100000)=274935105697756965126775163209863526881734293159800547582031259843021473281149641730550507416607366215901578447742962489404930
6307020046179276449303351011607934245719015571894350972531246610845200636955893446424871682878983218234500926285383140459702130713067451
0624419227311238999702284408609370935531629697851569569892196108480158600569421098519
     72.89
p(200000)=114214904051662974903712153040213809685392616088194853250500460143748696904706984823353353722854755809378873900261214264222037
9461285625596845952597486179818841905527749813150474498721747633064058840437916036245712660272843036962395856213719373419360640609193628
3942010549908977812630142747991846178699586488520943183046810987373649389305433467051476030145429989912526319600515962135313873380081652
17582822943454514251881668474896643416706191144360737225405767060795784657847417376966116043766
    149.90
p(300000)=707231191253211783269678797371648021355896898801778488819878481853999591121010342062777960135585544112696500275191068402829440
2958885133222880074150887104007189077222537716902899257410868565252976053292521902335261302640972607035941095443514511533394363682942044
7207620712300342910926227643727652063531315625542043533937391579914697595847692642602389304456413163890291950196347957519632703663009738
5442027675365492269983966956035293393526033624448380678592257208073583388792470728259656509767379603310147225241433214566875901141022409
3900667991718349432829214751081017558623916283527343825597038266403824
    248.04
p(400000)=131200193097519184547878054955101523969171484292611852019673890204125685362946566888203283508608026475734774436996428503963079
2644773245582092483624221442799304752281482750918438418010738139288648360586639566832051525393384597757564679076197280864701715813151614
7517669951797522809107918018308768874775535533269940818393263789280344720531854046250814604053450678695292843726918048226063314668912747
4162328855306490311315894888475657543504270307118284195854530558217484033107741325064923991690400741985013058560530156526463591901662133
0066245699646662742339077184995558642483191207042436892569973284574717154558133505876638000247857437268700170653948348813533331685742280
32077643798204985841890352489
    368.64
p(500000)=152472756837487393980836253932537739673358531198624515103707682257710225250432128794134234098309674791015637286141603020954224
0570198952049062704683020278422858152165197916558694852049725432428468932232809302094994003558655152182308393379265430858610663562220084
3507950426396914371150605581672698723453938785249080739534342453230387461061277399261986845473149706190286966253477258541270882378883038
9131560752520336948035125650561579102063523646545804292344473557497602227788254852943345757524511120797435403640616000889676211666179685
8975726034154016220298006489130140671909679401818357659146000508583759348862721560808802152150584981450693146410992733843661469688507084
0612645816218904303599258957453134962007342296918866491353916260421455682165182980786133075375269887558885137123
    514.65
p(600000)=194246109953425535538050608716833313396760167247956881113566053640960778746649992380574021065957325255482366279583079626334400
0688235360233840020986998982759848171576576503540053469470830309947794788694912557093672382060447268558338419878790760653954268178713601
4120934065332690529031757035114011819456723132689262150648243072153294195790294829034924369883206827035292484379465645001669708133297661
1249381724677646702226265374944993176151917956330653856014096433998657182921178334418204772834592691810950330470486832118960271990041991
5527764553402654256368889952751546017985199220130284902247676779249347407057107341201821068774049492494146193396491679668737771739509801
3756324832271756855640305349659695029010714590076546800982604426447833131100163345287291437351284584855036863627165183356275421878705480
893353311985442330043653016070261615365333691495637
    682.41
p(700000)=229371517575093341030303173260855451538391945689863147378306956440664883199111870795895016480927919733233856340398378553790157
4560239228471318567185321614031788172462639580069516320813704695416217435179095144772222273959484121207858450283350470476285629088838370
7875844704162051242064914616706238063874502023089918942926028965702496122381395801596808626736665073503189726899177978407931952257219118
1762286727797344114191762007444536371919411895173814674238879831814870705670874101887160393208817048256435351556165930105155926803911802
5778693286686204304867015706348155812737790256376101759006630203965349518888800636895096881276761358805029990838365621668009901244252923
0239726904289881227662033018415184665423208504925249680122383585988099324957423134700511513089640089866896499993433340470155425943627334
285533173550243382268370388958483310219360523760553623598753604832680167701602234619998630939112954232446171826011381608
    871.90
p(800000)=452566116823565258695585068696501113551717307775965307839936772817662723867298147782035052972999652449565566480121471400679992
3142642884230142240285419733922111198991121904425031982421132793802403522629882089802690359826076813154197189622894204942816158220102334
3198943500601994260722041662443491441851835767943782318542930810707851765147667131842900875918988697668635814760171712037643971847944243
9903097579131304595353855909986282573203048998607543555666472097464619015880537531306706851151504161575846706154447997071036510246659554
1925903764642221186746067108023567886278585476505064860119203131644821691871797736473738744178682688735719246878718774066983123016782910
9490129928995347846256963069724317948617603992557775724947494501795910056392079004567103501552485630591265550850246898497192739555641158
9913803601716794605147261863986698025440650254373686947999390947666426493916211310853508357868547934792444613705925296116184005379890208
164862104062700863595375993253089140776732271339
   1113.99
p(900000)=111253376808453434190074536845794182492615583102980174640471276160680294890818615369488895110709673026270698330209930093006965
7123399539425996938163849774130651871282184105332428951388806198214248447177010205227394456591481924133550557487870755593081999486524151
7289250056469136126791621951367425505582225548698464851577245461814633846498188124025578631145855279971019378403657609330528610674974394
6899328362662438340629571983750051496442134373915210853283593676547559491967144602985423116300931279889283440148547352699521404393618884
8110077795083826516203167823415669903459123587264733983662316564782929538051030790002733420943314117257370282117829042101432728056593519
2509733128123058326329972603557396851749255678834621605094639942385825609127826268937863582427108482079618143972966986648159265412893614
3932817315940032559163959321487618189133976859493444566715291415485057541162953192544661053959811358767630276192942531425842033406724123
2677926500662564812353438383656458699160401198447547473389224467173878325103717499772856329196190504992341801
   1362.47
p(1000000)=14716849863582233986310047606098959434840304844391421253346127473516661174189186182763301488739835975558420153741306002880959
2938734712823227032784957800193278439607206422865904871302017097184076102567647986084690814282935670692978599129051989944549067221999782
3452874982974022288229850136767566294781887494687879003824699988197729200632068668735996662273816798266213482417208446631027428001918132
1981771806465112345425950267284244525922967811934481399946647301057425643591547949891814852853513705513994767199816914590220155991019596
0141747407571543075002218489581520933901248173446944831932328015066538404299405417958775176129491624814247999880293650719525707448504757
1662771763903391442495113823298195263008336489826045837712202455304996382144601028531832004519046591968302787537418118486000612016852593
5427419802150462672454732373218458334275125242274653991301740769412808474008315422179992860711083363033162982891024446496968053954167918
7548001085263677402202312846764691977502234856252074774184334365780153413070476197553037516970799928704028567784161934747236817177215404
6664303121315630003467104673818

FMLIB によるもの 2

分割を求める関数をサブルーチン化しました。

    module sub_m
        use, intrinsic :: iso_fortran_env
        use FMZM
        implicit none
        integer, parameter :: ki = int64, kd = real64
    contains
        
        type(im) function p(n)
            integer, intent(in) :: n
            type(im) :: m, mp(0:n), isgn
            integer(ki) :: i, j, k  

            mp(0) = to_im(1)
            do i = 1, n
                mp(i) = to_fm(0)
          
                isgn = 1
                do j = 1, i
                    k = i - j * ( 3 * j - 1) / 2
                    if ( k < 0 ) exit
                    isgn = -isgn
                    mp(i) = mp(i) - isgn * mp(k) 
                end do
             
                isgn = 1
                do j = 1, i
                    k = i - j * ( 3 * j + 1) / 2
                    if ( k < 0 ) exit
                    isgn = -isgn
                    mp(i) = mp(i) - isgn * mp(k)
                end do
            end do
            p = mp(n)
        end function p
    end module sub_m
    
    
    program test
        use sub_m
        implicit none
        integer, parameter :: np = 10**6  
        character (1000) :: str
        type(im) :: mp
        
        do i = 0, np, np              
            mp = p(i)
            call im_form('i1000', mp, str)
            call system_clock(i1)
            write(*, '(f8.2, a, i0, a, a)') real(i1 - i0) / ic , ': p(', i, ')=', trim(adjustl(str)) 
        end do
    end program test

三角数経由

五角数を三角数から求めてやる方法。アルゴリズムが手計算向きの素朴なもの。8byte 整数の範囲での計算。(405 で overflow)

0  1  2  3  4  5  6  7  8  9...
0  1  3  6 10 15 21 28 36 45
0  x  1  2  x  5  7  x 12 15  
+     -  -     +  +     -  -
0=p(n)-p(n-1)-p(n-2)+p(n-5)+p(n-7)-p(n-12)-p(n-15)...

三角数を求め（累積和）それを 3 で割り、割り切れないものは捨てる、すると五角数が出てくる。二個づつ組で符号の因子が反転する。

    module sub_m
        use, intrinsic :: iso_fortran_env
        implicit none
        integer, parameter :: ki = int64
    contains
        integer(ki) pure function p4(n)
            integer(ki), intent(in) :: n
            integer(ki) :: i, j, isgn, k, ipenta, itrig, p(0:n)
      
            p(0) = 1 ! 1 partition, i sigma_1
            do i = 1, n
                p(i) = 0
                k = 0
                isgn = 1
                itrig = 0
                do j = 1, 2 * i
                    itrig = itrig + j
                    if (mod(itrig, 3) /= 0) cycle
                    ipenta = itrig / 3
                    if (ipenta > i) exit
                    k = k + 1
                    p(i) = p(i) + isgn * p(i - ipenta)
                    if (mod(k, 2) == 0) isgn = -isgn
                end do    
            end do  
            p4 = p(n)
        end function p4
    end module sub_m
    
    
    program partition_main
        use sub_m
        implicit none
        integer(ki) :: i

        do i = 1, 410
            print *, i, p4(i)
        end do
    end program partition_main

まず五角数の表を作る方法

プログラムはスッキリしますが、ワーク配列が必要になります。（大した量ではないですが）

    module sub_m
        use, intrinsic :: iso_fortran_env
        implicit none
        integer, parameter :: ki = int64
    contains
        integer(ki) pure function p4(n)
            integer(ki), intent(in) :: n
            integer(ki) :: i, j, k, penm(n), penp(n), p(0:n)
      
            do i = 1, n
                penm(i) = i * (3 * i - 1) / 2 ! pentagonal numbers
                penp(i) = i * (3 * i + 1) / 2 
            end do  
       
            do i = 1, n
                p(0) = 1 ! 1 partition, i sigma_1
                p(i) = 0 
                do j = 1, int( ( 1 + sqrt(1.0 + 24 * i)) / 6 )
                    p(i) = p(i) - (-1)**j * p(i - penm(j)) 
                end do  
                do j = 1, int( (-1 + sqrt(1.0 + 24 * i)) / 6 )
                    p(i) = p(i) - (-1)**j * p(i - penp(j)) 
                end do
            end do
            p4 = p(n)
        end function p4
    end module sub_m
    
    
    program partition_main
        use sub_m
        implicit none
        integer(ki) :: i
        
        do i = 1, 100
            print *, i, p4(i)
        end do
    end program partition_main

Modern Fortran Explained: Incorporating Fortran 2018 (Numerical Mathematics and Scientific Computation)

作者:Metcalf, Michael,Reid, John,Cohen, Malcolm
発売日: 2018/11/06
メディア: ハードカバー

Julia への移植

五角数を計算しながら求める方法

function p6(n)
    p = zeros(BigInt, n+1)
    p[1] = 1
    for i = 1:n

        isgn = 1
        for j = 1:i
            k = i - j*(3j-1)÷2
            k < 0 && break
            isgn = -isgn
            p[i+1] -= isgn * p[k+1]
        end
        
        isgn = 1
        for j = 1:i
            k = i - j*(3j+1)÷2
            k < 0 && break
            isgn = -isgn
            p[i+1] -= isgn * p[k+1]
        end
    end
    p[n+1]
end

1から始める Juliaプログラミング

作者:進藤裕之,佐藤建太
発売日: 2020/03/26
メディア: 単行本

BigInt の計算ですが Fortran の多倍長ルーチンで遅いと言われている FMLIB よりも時間がかかっているので、何か問題があるかもしれません。素直に Fortran 版を移植したつもりですが･･･配列の更新などでコピーを毎回取っているのかもしれません。よく分かりません。

@time p6(10^5)

 51.417361 seconds (274.62 M allocations: 10.206 GiB, 16.42% gc time)
27493510569775696512677516320986352688173429315980054758203125984302147328114964173055050741660736621590157844774296248940493063070200461792764493033510116079342457190155718943509725312466108452006369558934464248716828789832182345009262853831404597021307130674510624419227311238999702284408609370935531629697851569569892196108480158600569421098519

@time p6(10^6)
1975.618013 seconds (8.70 G allocations: 751.678 GiB, 18.74% gc time)
1471684986358223398631004760609895943484030484439142125334612747351666117418918618276330148873983597555842015374130600288095929387347128232270327849578001932784396072064228659048713020170971840761025676479860846908142829356706929785991290519899445490672219997823452874982974022288229850136767566294781887494687879003824699988197729200632068668735996662273816798266213482417208446631027428001918132198177180646511234542595026728424452592296781193448139994664730105742564359154794989181485285351370551399476719981691459022015599101959601417474075715430750022184895815209339012481734469448319323280150665384042994054179587751761294916248142479998802936507195257074485047571662771763903391442495113823298195263008336489826045837712202455304996382144601028531832004519046591968302787537418118486000612016852593542741980215046267245473237321845833427512524227465399130174076941280847400831542217999286071108336303316298289102444649696805395416791875480010852636774022023128467646919775022348562520747741843343657801534130704761975530375169707999287040285677841619347472368171772154046664303121315630003467104673818

まず五角数の表を作る方法

function p5(n)
    penta_m = zeros(BigInt, n)
    penta_p = zeros(BigInt, n)
    for i = 1:n
        penta_m[i] = i*(3i-1) ÷ 2
        penta_p[i] = i*(3i+1) ÷ 2
    end
    
    p = zeros(BigInt, n+1)
    p[1] = 1
    for i = 1:n
        p[i+1] = 0
        for j = 1:floor(BigInt, ( 1+sqrt(1+24i))÷6)
            p[i+1] -=(-1)^j * p[i+1-penta_m[j]]
        end
        for j = 1:floor(BigInt, (-1+sqrt(1+24i))÷6)
            p[i+1] -=(-1)^j * p[i+1-penta_p[j]]
        end
    end
    p[n+1]
end

@time p5(10^5)

170.563968 seconds (1.31 G allocations: 26.976 GiB, 18.67% gc time)
27493510569775696512677516320986352688173429315980054758203125984302147328114964173055050741660736621590157844774296248940493063070200461792764493033510116079342457190155718943509725312466108452006369558934464248716828789832182345009262853831404597021307130674510624419227311238999702284408609370935531629697851569569892196108480158600569421098519

三角数を経由する方法

function p4(n)
    p = zeros(BigInt, n+1)
    p[1] = 1
    for i = 1:n
        p[i+1] = 0
        k = 0
        isign = 1
        itrig = 0
        for j = 1:2i
            itrig += j
            mod(itrig, 3) != 0 && continue 
            ipenta = itrig ÷ 3
            ipenta > i && break 
            k += 1
            p[i+1] += isign * p[i+1 - ipenta]
            if mod(k, 2) == 0 
                isign = -isign
            end
        end
    end
    p[n+1]
end

分割数の直接計算

分割数を求める第五の方法として、 Hardy-Ramanujan-Rademacher の公式でも計算します。以前のプログラムを少し書き直してみました。

fortran66.hatenablog.com

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

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+0im
    for h = 1:k 
        if (gcd(h,k)==1) 
           z = -2n*h/k + fb(k, h)
           s += exp(im*π*z)
        end    
    end 
    s
end

function p(n)
    s = 0+0im
    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 %20.7f \n", i, round(Int64, real(pn)), real(pn), imag(pn)) 
end

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