読者です 読者をやめる 読者になる 読者になる

Maxima で綴る数学の旅

紙と鉛筆の代わりに、数式処理システムMaxima / Macsyma を使って、数学を楽しみましょう

-数学- q級数の計算の高速化

にほんブログ村 科学ブログ 数学へ
にほんブログ村

詳しい(そして実に面白いのですが)ことは、tsujimotterさんのブログ記事:


上記の記事をぜひ読んでください。

で、その中で、素数 \(p \) が \( p=6\,x^2+x\,y+y^2 \) と書けるための必要十分条件が、

$$ \displaystyle q\prod_{n=1}^{\infty} (1-q^n)(1-q^{23n}) = \sum_{n=1}^{\infty} a_n q^n $$

の左辺を展開した右辺の級数の係数\( a_n \) と密接な関係があることが書かれております。で是非とも\( a_n \) を計算したいわけです。

 

tsujimotterさんは上記のブログ記事の中で、

a[n] := (1-q^n)*(1-q^(23*n));

expand(q*product(a[n],n,1,100);

で計算しようとして、計算が爆発してしまったようです。

 

以下のようにプログラムを書いて、途中で現れる不要な項をどんどん消しながら計算を進めることでqの1000乗の項までの展開係数を正しく求めることができます。

 

(%i1) dcproduct(expr,var,lo,hi,var2,maxdeg):=
if (lo=hi) then subst(lo,var,expr) else
degcut(expand(subst(hi,var,expr)*dcproduct(expr,var,lo,hi-1,var2,maxdeg)),var2,maxdeg)$

dcproduct(expr,var,lo,hi,var2,maxdeg)という関数はexprの中の変数varをloからhiまで変化させながら掛け合わせていく関数です。product()とほとんど同じ動作をするのですが、var2という変数の次数がmaxdegを超えた項が途中の計算で現れると(計算に関係なくなるので)degcut() を使って消してしまいます。


(%i2) degcut(expr,var,deg):=
if atom(expr) then if hipow(expr,var)>deg then 0 else expr else
map(lambda([exp],if hipow(exp,var)>deg then 0 else exp),expr)$

このdegcut(expr,var,deg)はexprがvarの多項式とみて、次数がdegを超える項を0に置き換えて消してしまいます。


(%i3) q*dcproduct(expand((1-q^n)*(1-q^(23*n))),n,1,1000,q,1000), expand;
$$ \tag{%o3} 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}\boldsymbol{+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}\boldsymbol{+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}\boldsymbol{+2}\,q^{167}\boldsymbol{+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}\boldsymbol{+2}\,q^{211}+q^{213}+q^{216}+q^{219}\boldsymbol{+2}\,q^{223}-q^{232}-q^{233}-q^{239}-q^{242}-q^{246}-q^{248}+q^{254}-q^{257}+q^{262}-q^{269}\boldsymbol{+2}\,q^{271}-q^{277}+q^{278}-q^{282}+q^{289}+q^{294}-q^{299}+q^{302}-2\,q^{303}\boldsymbol{+2}\,q^{307}-q^{311}+q^{312}\boldsymbol{+2}\,q^{317}-q^{325}+q^{326}-q^{328}-q^{331}-2\,q^{334}-2\,q^{346}\boldsymbol{+2}\,q^{347}-q^{349}-q^{351}-q^{353}\boldsymbol{+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}\boldsymbol{+2}\,q^{449}+q^{453}-q^{461}\boldsymbol{+2}\,q^{463}+q^{464}+q^{466}\boldsymbol{+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}\boldsymbol{+2}\,q^{593}+q^{598}\boldsymbol{+2}\,q^{599}-q^{600}-q^{601}\boldsymbol{+2}\,q^{606}\boldsymbol{+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}\boldsymbol{+2}\,q^{691}-2\,q^{694}+q^{696}+q^{698}+q^{699}+q^{702}+q^{706}-q^{713}+q^{717}\boldsymbol{+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}\boldsymbol{+2}\,q^{808}\boldsymbol{+2}\,q^{809}-q^{811}-2\,q^{813}+q^{818}\boldsymbol{+2}\,q^{821}-q^{823}\boldsymbol{+2}\,q^{829}+q^{831}-q^{832}-q^{834}-q^{837}\boldsymbol{+2}\,q^{853}-q^{857}-q^{859}-q^{863}-q^{867}\boldsymbol{+2}\,q^{877}+q^{878}\boldsymbol{+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}\boldsymbol{+2}\,q^{991}+q^{993}\boldsymbol{+2}\,q^{997}+q^{998} $$
(%i4) time(%o3);
$$ \tag{%o4} \left[ 856.36 \right] $$

計算に14分ほどかかりました。

This is a test.
天気