オイラーの五角数を用いた分割数の漸化式
但し とする。
Julia の BigInt の試用目的で昔の計算をやり直してみます。
fortran66.hatenablog.com
fortran66.hatenablog.com
fortran66.hatenablog.com
まずは Fortran で書いてそれを移植する形で進めます。昔、書いたものを修正していくつか書いてみます。
FMLIB によるもの 1
https://dmsmith.lmu.build/FMLIB.html
main ルーチン中で途中経過を出力しながら 106 の分割数を計算します。
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
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
penp(i) = i * (3 * i + 1) / 2
end do
do i = 1, n
p(0) = 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
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
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
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