お花見

 

q級数の展開計算をしたい、と思ったのです。以前に書いた記事 

ではMaximaのプログラムで書いてみたわけですが、1000次まで正確に計算するのに8分もかかるという遅さでした。これでは色々と試すには遅すぎですよね。

そこで今回全部Lisp言語で書き直してみました。

さすがにそのプログラムが見たい、という人は少ないと思うので、ここには掲載しません。100行弱のCommon Lispのプログラムです。CLOS (Common Lisp Object System)を使って書いてあります。

まずは読み込んで、

(%i1) load("/Users/yasube/Dropbox/qexpand.lisp")$
(%i2) powerdisp:true;
$$ ag{%o2} mathbf{true} $$

tsujimotterさんの記事で出てきた無限積をやってみます。
(%i3) q*product((1-q^n)*(1-q^(23*n)),n,1,inf);

$$ ag{%o3} q, prod_{n=1}^{infty }{left(1-q^{n} ight),left(1-q^{23,n} ight)} $$

 

今回作ったqsexpand()が高速展開関数です。第1引数には展開したいqの無限積を、第2引数には求めたい最高次数を指定します。

(%i4) qsexpand(%,1000);

$$ ag{%o4} q-q^2-q^3+q^6+q^8-q^{13}-q^{16}+q^{23}-q^{24}+q^{25}+q^{26}+q^{27}-q^{29}-q^{31}+q^{39}-q^{41}-q^{46}-q^{47}+q^{48}+q^{49}-q^{50}-q^{54}+q^{58}+2,q^{59}+q^{62}+q^{64}-q^{69}-q^{71}-q^{73}-q^{75}-q^{78}-q^{81}+q^{82}+q^{87}+q^{93}+q^{94}-q^{98}+2,q^{101}-q^{104}-2,q^{118}+q^{121}+q^{123}-q^{127}-q^{128}-q^{131}+q^{138}-q^{139}+q^{141}+q^{142}+q^{146}-q^{147}+q^{150}-q^{151}+q^{162}-q^{163}+2,q^{167}+2,q^{173}-q^{174}-2,q^{177}-q^{179}+q^{184}-q^{186}-q^{192}-q^{193}-q^{197}+q^{200}-2,q^{202}+q^{208}+2,q^{211}+q^{213}+q^{216}+q^{219}+2,q^{223}-q^{232}-q^{233}-q^{239}-q^{242}-q^{246}-q^{248}+q^{254}-q^{257}+q^{262}-q^{269}+2,q^{271}-q^{277}+q^{278}-q^{282}+q^{289}+q^{294}-q^{299}+q^{302}-2,q^{303}+2,q^{307}-q^{311}+q^{312}+2,q^{317}-q^{325}+q^{326}-q^{328}-q^{331}-2,q^{334}-2,q^{346}+2,q^{347}-q^{349}-q^{351}-q^{353}+2,q^{354}+q^{358}+q^{361}-q^{363}-q^{368}-q^{376}+q^{377}+q^{381}+q^{384}+q^{386}+q^{392}+q^{393}+q^{394}-q^{397}-q^{400}+q^{403}-q^{409}+q^{417}-2,q^{422}-q^{426}-q^{432}-q^{438}-q^{439}-q^{443}-2,q^{446}+2,q^{449}+q^{453}-q^{461}+2,q^{463}+q^{464}+q^{466}+2,q^{472}+q^{478}-q^{487}+q^{489}-q^{491}+q^{496}-q^{499}-2,q^{501}-q^{509}+q^{512}+q^{514}-2,q^{519}+q^{529}+q^{533}+q^{537}+q^{538}-q^{541}-2,q^{542}-q^{547}-q^{552}+q^{554}-q^{568}+q^{575}-q^{577}-q^{578}+q^{579}-q^{584}-q^{587}+q^{591}+2,q^{593}+q^{598}+2,q^{599}-q^{600}-q^{601}+2,q^{606}+2,q^{607}+q^{611}-2,q^{614}+q^{621}+q^{622}-q^{624}+q^{625}-2,q^{633}-2,q^{634}-q^{637}-q^{647}-q^{648}+q^{650}-q^{653}+q^{656}+q^{662}-q^{667}-2,q^{669}-q^{673}+q^{675}-q^{683}+2,q^{691}-2,q^{694}+q^{696}+q^{698}+q^{699}+q^{702}+q^{706}-q^{713}+q^{717}+2,q^{719}-q^{722}-q^{725}+q^{726}+q^{729}-q^{739}+q^{744}+q^{752}-q^{754}-q^{761}-q^{762}-2,q^{767}+q^{771}-q^{775}-q^{783}-q^{784}-q^{786}+q^{794}-q^{806}+q^{807}+2,q^{808}+2,q^{809}-q^{811}-2,q^{813}+q^{818}+2,q^{821}-q^{823}+2,q^{829}+q^{831}-q^{832}-q^{834}-q^{837}+2,q^{853}-q^{857}-q^{859}-q^{863}-q^{867}+2,q^{877}+q^{878}+2,q^{883}+q^{886}-q^{887}+q^{897}-2,q^{898}+q^{899}-q^{906}-2,q^{921}+q^{922}+q^{923}-2,q^{926}-q^{929}+q^{933}-q^{943}-2,q^{944}-q^{947}+q^{949}-2,q^{951}-q^{967}+q^{968}+q^{974}+q^{975}-q^{978}+q^{982}+q^{984}+2,q^{991}+q^{993}+2,q^{997}+q^{998} $$
(%i5) time(%);
$$ ag{%o5} left[ 0.18 ight] $$

たった0.18秒で展開できました。maximaで書いたプログラムでは850秒位かかっていましたから実に4700倍くらいの高速化が出来ました!

 

次にラマヌジャンの保型形式を100次まで展開してみます。
(%i6) q*product((1-q^n)^24,n,1,inf);
$$ ag{%o6} q, prod_{n=1}^{infty }{left(1-q^{n} ight)^{24}} $$
(%i7) qsexpand(%,100);
$$ ag{%o7} q-24,q^2+252,q^3-1472,q^4+4830,q^5-6048,q^6-16744,q^7+84480,q^8-113643,q^9-115920,q^{10}+534612,q^{11}-370944,q^{12}-577738,q^{13}+401856,q^{14}+1217160,q^{15}+987136,q^{16}-6905934,q^{17}+2727432,q^{18}+10661420,q^{19}-7109760,q^{20}-4219488,q^{21}-12830688,q^{22}+18643272,q^{23}+21288960,q^{24}-25499225,q^{25}+13865712,q^{26}-73279080,q^{27}+24647168,q^{28}+128406630,q^{29}-29211840,q^{30}-52843168,q^{31}-196706304,q^{32}+134722224,q^{33}+165742416,q^{34}-80873520,q^{35}+167282496,q^{36}-182213314,q^{37}-255874080,q^{38}-145589976,q^{39}+408038400,q^{40}+308120442,q^{41}+101267712,q^{42}-17125708,q^{43}-786948864,q^{44}-548895690,q^{45}-447438528,q^{46}+2687348496,q^{47}+248758272,q^{48}-1696965207,q^{49}+611981400,q^{50}-1740295368,q^{51}+850430336,q^{52}-1596055698,q^{53}+1758697920,q^{54}+2582175960,q^{55}-1414533120,q^{56}+2686677840,q^{57}-3081759120,q^{58}-5189203740,q^{59}-1791659520,q^{60}+6956478662,q^{61}+1268236032,q^{62}+1902838392,q^{63}+2699296768,q^{64}-2790474540,q^{65}-3233333376,q^{66}-15481826884,q^{67}+10165534848,q^{68}+4698104544,q^{69}+1940964480,q^{70}+9791485272,q^{71}-9600560640,q^{72}+1463791322,q^{73}+4373119536,q^{74}-6425804700,q^{75}-15693610240,q^{76}-8951543328,q^{77}+3494159424,q^{78}+38116845680,q^{79}+4767866880,q^{80}+1665188361,q^{81}-7394890608,q^{82}-29335099668,q^{83}+6211086336,q^{84}-33355661220,q^{85}+411016992,q^{86}+32358470760,q^{87}+45164021760,q^{88}-24992917110,q^{89}+13173496560,q^{90}+9673645072,q^{91}-27442896384,q^{92}-13316478336,q^{93}-64496363904,q^{94}+51494658600,q^{95}-49569988608,q^{96}+75013568546,q^{97}+40727164968,q^{98}-60754911516,q^{99} $$
(%i8) time(%);
$$ ag{%o8} left[ 0.01 ight] $$

このくらいの速度になれば色々と実験してみたくなります。