C言語で書かれたプログラム(4146.)の計算部分の説明に移ります。
関数init1()はmain()から呼び出され、ビットマップの1024×1024点の個々の計算に移る前に、できるだけの前計算をする部分です。
まず、最初に設定されたk値(母数)からl値を経てq値を算出します。この部分は公式集のそのままです。q値から出発することも可能で、その場合は上記のk→l→qの計算をコメントアウトします。k値は後の計算で必要なので、すぐ下で楕円テータ関数経由で再計算されています。
変数、te2a, te3a, te0aはそれぞれ、θ2(0), θ3(0), θ0(0)の事です。元の公式ではqの級数になっていて、プログラムでは展開して、最初の10項を計算しています。k^2 = 0.5を中心とし、0～1の広範囲で最初の数項しか役立ちません。律儀にやるなら項が0に十分近くなったらそこで打ち切るべきですが、それぞれ1回しか計算しないのでそのままにしています。逆に、64bit IEEE倍精度浮動小数点演算の精度では10項より多く計算しても誤差が累積するので無駄のようです。
次に、k'に対応するq値、q'を算出して、k'に対応する楕円テータ関数の定数`θ2(0), `θ3(0), `θ0(0)を同様に計算します(tf2a, tf3a, tf0a)。
定数、θ20, θ30, θ00が得られたので、k (elk)、k'(elkc)、K(ellk)、K'(ellkc)を算出します。
K ≠ K'の場合はuとvの間隔を1:1にするために、wとzの縦横比を修正しないといけません。Kが大きい場合は縦(VD)を縮小し、逆にK'が大きい場合は横(UD)を縮小します。
次がメインの、θn(w)と`θn(z)の計算部分です。実は最初の段階では全ピクセルでこの級数の計算が必要と思い込んでいたので、やや冗長な記述となりました。元来は打ち切り項数を4段階くらいに分けようとしていました。
しかし、加法定理のおかげで、縦横4ラインずつの計算で済みます。横が実部のte1[i], te2[i], te3[i], te0[i]で横のピクセル数分、縦が虚部のtf1[i], tf2[i], tf3[i], tf0[i]で縦のピクセル数分、の計算量です。上述の理由によって、楕円テータ関数の級数、10項分の計算に固定しています。
分母は共通ですが、θn(w)と`θn(z)が混じっているのであらかじめ計算出来るのは自乗だけです。乗算は高速で、キャッシュなどを考えると節約になっているのかどうか怪しいのですが、あらかじめ計算してメモリに置く方を選択しました。
関数draw1()はビットマップの1点ずつの計算部分です。中身は単なる二重ループで、点を打つのは関数drawsn1(), drawcn1(), drawdn1()です。動作は同様なのでdrawsn1()で説明します。
まず、(a + bi) / dのそれぞれの値に相当する、sn1, sn2, sn0を計算します。後は順番に、
分母が極端に小さい値の場合は極とみなす。
分子が極端に小さい値の場合は零点とみなす。
分子/分母が南極域に入っている(値が小さい)場合は零点の色にする。
分母/分子が北極域に入っている(値が小さい)場合は極の色にする。
その後は普通に実部と虚部が計算出来ますから、
縦横の隣接するドットと実部の符号が異なる場合は純虚数がある。
縦横の隣接するドットと虚部の符号が異なる場合は純実数がある。
北半球と南半球に分けて、絶対値が等高線の範囲に入っている場合はその色づけをする。
残りのドットは背景扱いで、値の複素数の偏角に応じた色づけをする。
後はmain()の後半で絵をファイル(sn.bmp, cn.bmp, dn.bmp)に書き出し、定数などの必要な情報をテキストファイル(el.txt)に出力します。
