4149. 楕円関数、複素平面、上級編、続き^3

　プログラムでの計算のやり方の解説なので、4139.の続きです。

　ヤコビの楕円関数、sn(u, k)などと書く時、グラフの形を決めるkは母数(modulus)と呼ばれ、実数(0 ≦ k ≦ 1)です。uが関数の入力値で、実数のこともあれば複素数のこともあり、少なくとも本ブログの実数編では実数→実数の関数でした。

　加法定理の公式は複素数でも成り立つので、sn(u + iv)で実数部分と虚数部分に分離が出来て、虚数部分はヤコビの虚変換の公式で補母数k'を用いた実数→実数の関数で計算出来るようになります。
　具体的計算には楕円テータ関数を仲介とすると効率が良く、こちらではkと同じ働きをする係数がτ(タウ)で表され、これは二重周期の周期の比率です。τは一般には普通の複素数ですが、ここでの議論ではτ = iK'/Kなので純虚数となります。Kはヤコビの楕円関数の第一周期で、K'は第二周期で、どちらも実数(kとk'(kの補母数、√(1-k^2) (^は累乗))を用いた完全楕円積分の値)です。
　計算時にはτから計算されるq値が用いられて、楕円テータ関数はqの級数の形になっています。τが純虚数の場合はqは実数になります。

　楕円テータ関数は4種あり、θ1(v, τ)、θ2(v, τ)、θ3(v, τ)、θ0(v, τ)がフルの表記となります。文脈に応じて楕円関数のkが省略されるのと同様にτも普通は省略され、θ1(v)、θ2(v)、θ3(v)、θ0(v)の表記になります。テータは異体字を使うこともあり、普通のθを使うこともあります。
　このvがくせ者で、ヤコビの楕円関数のuとは、v = u/(2K)の関係があります。ところが上述のsn(u + iv)のvと記号が重複してしまいます。さらに、楕円テータ関数は本来は複素関数なので、変数用の文字が足りません。

　混乱を防ぐには慣例を無視せざるを得ず、以下とりあえず、

　w = u/(2K), z = v/(2K')

　と表記することにします。u, v, w, zはいずれも実数です。

　θ2(0)、θ3(0)、θ0(0)は定数としてよく出てくるので、括弧を省略してθ20、θ30、θ00と書くことにします。公式集では数字の添え字は足文字と肩文字が使われます。
　さらに、k'の時のq値をq'と表記します。q'に対応するテータ関数は、θ'では導関数になってしまうので、ここでの独自表記、`θを使用することにします。

　やっと本題に入れます。結論だけ書くと、

　sn(u + iv) = [(k'/k) θ1(w) θ0(w) `θ3(z) `θ0(z) + (ik'/k) θ2(w) θ3(w) `θ1(z) `θ2(z)] / [(θ1(w))^2 (`θ1(z))^2 + (θ0(w))^2 (`θ2(z))^2]

　cn(u + iv) = [(k'/k) θ2(w) θ0(w) `θ2(z) `θ0(z) + (-ik'/k) θ1(w) θ3(w) `θ1(z) `θ3(z)] / [(θ1(w))^2 (`θ1(z))^2 + (θ0(w))^2 (`θ2(z))^2]

　dn(u + iv) = [(k') θ3(w) θ0(w) `θ2(z) `θ3(z) + (-ik') θ1(w) θ2(w) `θ1(z) `θ0(z)] / [(θ1(w))^2 (`θ1(z))^2 + (θ0(w))^2 (`θ2(z))^2]

　です。分母はどれも同じです。式の変換は難しくありません。ただし、途中で、

　k = (θ20/θ30)^2,  k' = (θ00/θ30)^2,  k'/k = (θ00/θ20)^2
　k' = (`θ20/`θ30)^2,  k = (`θ00/`θ30)^2,  k/k' = (`θ00/`θ20)^2

　を用います。
　これでめでたく、

　f(x + iy) = (a + bi) / d

　の形になりました。

4148. 楕円関数、複素平面、上級編、修正済み

　よく見たら、円分方程式の時と赤・緑が逆になっていました。明日にも修正します。実数の＋が赤系なのが主流のようなので、合わせます。

　追記: 修正が完了しました。修正したのは、4136.の3枚と4145.の6枚の図の入れ替えと偏角の色の説明の場所(草色と紅色の交換)、4137.と4146.のプログラムの3箇所の負号「-」の削除、です。

4147. 楕円関数、複素平面、上級編、続き^2

　プログラムはC言語で書かれています。コンパイラは私はMicrosoft Visual Studioを使っていて、無料版がオンラインで手に入るはずです。C++のコンソールプロクラムが開発できるようにダウンロード、インストールします。C++のコンソールプロクラムでHello, worldが出力されるプロジェクトを新規作成し、ソースプログラムを前回掲載したプログラムに置き換え、ビルド、実行します。指定のフォルダにグラフの3ファイル(bmp)と情報の1テキストファイル(txt)が作られているので、適当なソフト(最初からOSに含まれているものでOK)で表示させます。

　k値などのパラメータを変更すると再コンパイルが必要です。別に用意したファイルから読み込むとか、ウィンドウ・プログラミングすれば操作性は向上しますが、本質とは関係の無いUI部のプログラムが増えるので、とりあえずこのままとさせていただきます。

　C言語として難しい箇所は無いはずです。しかし、配列とポインタとポインタのポインタは出てきますので、慣れていないと読みにくいと思います。
　プログラムの最初の部分は、Windowsでサポートされているビットマップ形式のファイルを作るための定数とメモリ領域の確保です。円周率/4の定数pi4を定義している部分までがそう。最終部のmain()関数の大部分もビットマップファイルとテキストファイルの作成に費やされています。
　グラフの作成部分は、関数init1()と関数draw1()です。

　init1()の前に、プロクラムで使用する大域変数の定義があります。ここの指定部分を書き換えて再コンパイル、実行すると新たな模様のグラフが得られます。

　double elk = sqrt(0.9755);

　が楕円関数の母数(modulus) kの指定です。ここでは自乗値の0.9755を与えています。もちろん、k値を直接書いてもOKです。

　UCとVCはグラフの中心の値(関数の入力値)で、その実部と虚部です。単位は周期を表すKとK'です。もっと左右や上下の部分を見たい時に変更します。0でも良いのですが、多少の理由があってわずかにずらせた値を指定しています。

　UDとVDは隣接するピクセルとの値の差分値です。通常は同じ値を指定します。図を全体として拡張表示するにはこのピッチ値を(両者とも)大きく(遠くまで見えるので縮尺が大きくなる)し、細部を見たい時は(両者とも)小さい値にします。ただし、今回のプログラムでは、k値が√(0.5)でない場合は、K/K'が1でなくなるので周期が延びた方を視野に納めるために、それと垂直方向をプログラム内(行104付近)で縮めています。
　縮めるのはKとK'が単位となる楕円テータ関数への入力値なので、この補正で、元のヤコビの楕円関数(sn, cn, dn)の入力値に関しては実部と虚部は正しく1:1の比率になります。その結果として、グラフの等高線と偏角の線が垂直に交わる、本来の姿(等角写像)となっています(多分)。

　DDは関数値の絶対値の等高線の間隔です。DEは等高線の幅(間隔との比)です。等高線の間隔は、関数値のリーマン球面を地球儀に例えて、関数値1を赤道に、零点を南極に、極(無限遠点)を北極に対応させ、緯度で指定します。

　const double DD = pi4 * 15.0 / 90.0; // value hight pitch

　では15度間隔になっています。各等高線の緯度での具体的な関数値はテキストファイルに書き出されています。

　今回のプログラムでは、赤道での関数値を変更できるようにしました。これは、関数値1だけでなく、ヤコビの楕円関数では、snでは1/k、cnではk'/k、dnではk'が特別な値で、定義域の複素平面での場所とその周りの様子が知りたかったからです。少し下にある、

　const double alf = 1.0; // *** modify here

　の1.0を正の任意の実数に変更できます。任意の緯度の値との関係は線形なので、このalf値を単純に掛けると算出できます。計算後の値はテキストファイルに出ています。

　極と零点を表すドットの大きさをリーマン球面の北極と南極からの角度で指定できるようにしました。

　const double anp = tan(4.0 * pi4 / 90.0); // pole area (infinity) in degree
　const double asp = tan(4.0 * pi4 / 90.0); // zero point area in degree

　の部分です。これで4.0度の半径の円を指定したことになります。北緯86度より北は北極のエリア、南緯86度より南は南極のエリア、です。
　なお、本プログラムでは、計算途中で分母の値が極端に小さな場合は分子の多少に関わらず極扱いしています。　

　色の指定は前回のプログラムと同様です。極と零点のドットの色、純実数と純虚数の線の色、北半球と南半球の等高線の色が配列cl0[7][3]で青・緑、赤の順で0-255の数値で指定されています。また、偏角の虹模様の色は配列cl1[24][3]です。24エントリですが、12色が入っています。

　その下の、ヤコビの楕円関数の複素平面表示の計算本体の説明は次回回しとします。

4146. 楕円関数、複素平面、上級編、続き

　今回作成した、ヤコビの楕円関数の複素平面表示のプログラムです。適当なコンパイラでコンパイルし、実行するとグラフの3ファイル、sn.bmp, cn.bmp, dn.bmpがソースプログラムのあるフォルダに出力されます。同時に、グラフを読むのに必要な数値をまとめたテキストファイル、el.txtが出力されます。(230811改訂版)

// ConsoleApplication7.cpp
// ellptic3.cpp; Jacobi elliptic functions on a Gaussian plane
// based on bmp24.cpp
//  200828; making 24bit color bmp file
// 230711; 230722; 230729; 230811

#include <stdio.h>
#include <math.h>

const int WT = 1024; // width  (mod 8 = 0) *** modify here
const int HT = 1024; // height (mod 8 = 0) *** modify here
const int BZ = WT * HT * 3; // buffer size
const int FZ = BZ + 0x36; // file size

unsigned char b1[0x36];  // header*2
unsigned char b2[BZ];  // bitmap sn
unsigned char b3[BZ];  // bitmap cn
unsigned char b4[BZ];  // bitmap dn

const double pi4 = atan(1.0);

// elliptic k (0.0 <= k <= 1.0)
double elk = sqrt(0.9755); // k^2: ***modify here

// draw area definition *** modify here
// UD,VD will be automatically modified with K'/K
const double UC = 0.00125; // real base
const double VC = 0.00125; // imaginary base
double UD = 0.0025; // real pitch // hard
double VD = 0.0025; // imaginary pitch // hard

// contour information *** modify here
const double DD = pi4 * 15.0 / 90.0; // value hight pitch
const double DE = 0.125; // value edge width (0.01 - 0.99)

// pole area definition *** modify here
const double anp = tan(4.0 * pi4 / 90.0); // pole area (infinity) in degree
const double asp = tan(4.0 * pi4 / 90.0); // zero point area in degree

// Riemann sphere magnitude (absolute value at equator)
// it dose not adjusted by 1/k or k'/k or k'. please rewrite manually
// const double alf = sqrt(1.0 - elk * elk) / elk; // example
const double alf = 1.0; // *** modify here

// color index
const unsigned char cl0[7][3] = { // {b, g, r} 0-255
{ 0x2e, 0x2e, 0x2e }, // background
{ 0x66, 0x66, 0xbb }, // pole (infinity)
{ 0xbb, 0xbb, 0x22 }, // zero point
{ 0x11, 0x11, 0x11 }, // real part = 0; pure imaginary
{ 0xff, 0xff, 0xff }, // imaginary part = 0; pure real
{ 0x66, 0x66, 0xaa }, // contour edge nothern hemisphere
{ 0x99, 0x99, 0x66 } // contour edge southern hemisphere
};

const unsigned char cl1[24][3] = {
{ 0x99, 0xff, 0x33 },{ 0xff, 0xff, 0x33 },{ 0xff, 0xff, 0x33 },{ 0xff, 0x99, 0x33 },
{ 0xff, 0x99, 0x33 },{ 0xff, 0x33, 0x33 },{ 0xff, 0x33, 0x33 },{ 0xff, 0x33, 0x99 },
{ 0xff, 0x33, 0x99 },{ 0xff, 0x33, 0xff },{ 0xff, 0x33, 0xff },{ 0x99, 0x33, 0xff },
{ 0x99, 0x33, 0xff },{ 0x33, 0x33, 0xff },{ 0x33, 0x33, 0xff },{ 0x33, 0x99, 0xff },
{ 0x33, 0x99, 0xff },{ 0x33, 0xff, 0xff },{ 0x33, 0xff, 0xff },{ 0x33, 0xff, 0x99 },
{ 0x33, 0xff, 0x99 },{ 0x33, 0xff, 0x33 },{ 0x33, 0xff, 0x33 },{ 0x99, 0xff, 0x33 }
};

// elliptic function constants

double elq, elqc; // q, q'
double elkc, ellk, ellkc; // k', K, K'
double te2a, te3a, te0a, tf2a, tf3a, tf0a; // theta(0)s

// elliptic function lines
double te1[WT + 1], te2[WT + 1], te3[WT + 1], te0[WT + 1]; // using q
double tf1[HT + 1], tf2[HT + 1], tf3[HT + 1], tf0[HT + 1]; // using q'

// elliptic function buffers
double tea[WT + 1], teb[WT + 1], tfa[HT + 1], tfb[HT + 1]; // denominator: sn0 = cn0 = dn0

void init1(void) {
 int i;
 double u;

 // prepare constants
 double ell = (1.0 - pow(1.0 - elk * elk, 0.25)) / 2.0 / (1.0 + pow(1.0 - elk * elk, 0.25));
 elq = ell + 2.0 * pow(ell, 5.0) + 15.0 * pow(ell, 9.0) + 150.0 * pow(ell, 13.0) + 1707.0 * pow(ell, 17.0);
 // elq = 0.04321; // starting with q
 te2a = 2.0 * (pow(elq, 0.25) + pow(elq, 2.25) + pow(elq, 6.25) + pow(elq, 12.25) + pow(elq, 20.25)
  + pow(elq, 30.25) + pow(elq, 42.25) + pow(elq, 56.25) + pow(elq, 72.25) + pow(elq, 90.25));
 te3a = 1.0 + 2.0 * (elq + pow(elq, 4.0) + pow(elq, 9.0) + pow(elq, 16.0) + pow(elq, 25.0)
  + pow(elq, 36.0) + pow(elq, 49.0) + pow(elq, 64.0) + pow(elq, 81.0) + pow(elq, 100.0));
 te0a = 1.0 + 2.0 * (-elq + pow(elq, 4.0) - pow(elq, 9.0) + pow(elq, 16.0) - pow(elq, 25.0)
  + pow(elq, 36.0) - pow(elq, 49.0) + pow(elq, 64.0) - pow(elq, 81.0) + pow(elq, 100.0));
 elqc = exp(pi4 * pi4 * 16.0 / log(elq));
 tf2a = 2.0 * (pow(elqc, 0.25) + pow(elqc, 2.25) + pow(elqc, 6.25) + pow(elqc, 12.25) + pow(elqc, 20.25)
  + pow(elqc, 30.25) + pow(elqc, 42.25) + pow(elqc, 56.25) + pow(elqc, 72.25) + pow(elqc, 90.25));
 tf3a = 1.0 + 2.0 * (elqc + pow(elqc, 4.0) + pow(elqc, 9.0) + pow(elqc, 16.0) + pow(elqc, 25.0)
  + pow(elqc, 36.0) + pow(elqc, 49.0) + pow(elqc, 64.0) + pow(elqc, 81.0) + pow(elqc, 100.0));
 tf0a = 1.0 + 2.0 * (-elqc + pow(elqc, 4.0) - pow(elqc, 9.0) + pow(elqc, 16.0) - pow(elqc, 25.0)
  + pow(elqc, 36.0) - pow(elqc, 49.0) + pow(elqc, 64.0) - pow(elqc, 81.0) + pow(elqc, 100.0));
 elk = pow(te2a / te3a, 2.0); elkc = pow(te0a / te3a, 2.0);
 ellk = 2.0 * pi4 * pow(te3a, 2.0); // K
 ellkc = ellk / (4.0 * pi4) * log(1 / elq); // K'
 // automatic expansion
 if (ellk > ellkc) { VD *= ellk / ellkc; }
 else { UD *= ellkc / ellk; }

 // prepare lines
 for (i = 0; i <= WT; i++) {
  u = UC + UD * (double)(i - WT / 2);
  te1[i] = 2 * (pow(elq, 0.25) * sin(pi4 * 4.0 * u) - pow(elq, 2.25) * sin(pi4 * 4.0 * 3.0 * u)
   + pow(elq, 6.25) * sin(pi4 * 4.0 * 5.0 * u) - pow(elq, 12.25) * sin(pi4 * 4.0 * 7.0 * u)
   + pow(elq, 20.25) * sin(pi4 * 4.0 * 9.0 * u) - pow(elq, 30.25) * sin(pi4 * 4.0 * 11.0 * u)
   + pow(elq, 42.25) * sin(pi4 * 4.0 * 13.0 * u) - pow(elq, 56.25) * sin(pi4 * 4.0 * 15.0 * u)
   + pow(elq, 72.25) * sin(pi4 * 4.0 * 17.0 * u) - pow(elq, 90.25) * sin(pi4 * 4.0 * 19.0 * u));
  te2[i] = 2 * (pow(elq, 0.25) * cos(pi4 * 4.0 * u) + pow(elq, 2.25) * cos(pi4 * 4.0 * 3.0 * u)
   + pow(elq, 6.25) * cos(pi4 * 4.0 * 5.0 * u) + pow(elq, 12.25) * cos(pi4 * 4.0 * 7.0 * u)
   + pow(elq, 20.25) * cos(pi4 * 4.0 * 9.0 * u) + pow(elq, 30.25) * cos(pi4 * 4.0 * 11.0 * u)
   + pow(elq, 42.25) * cos(pi4 * 4.0 * 13.0 * u) + pow(elq, 56.25) * cos(pi4 * 4.0 * 15.0 * u)
   + pow(elq, 72.25) * cos(pi4 * 4.0 * 17.0 * u) + pow(elq, 90.25) * cos(pi4 * 4.0 * 19.0 * u));
  te3[i] = 1.0 + 2.0 * (elq * cos(pi4 * 8.0 * u) + pow(elq, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
   + pow(elq, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elq, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
   + pow(elq, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elq, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
   + pow(elq, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elq, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
   + pow(elq, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elq, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
  te0[i] = 1.0 + 2.0 * (-elq * cos(pi4 * 8.0 * u) + pow(elq, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
   - pow(elq, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elq, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
   - pow(elq, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elq, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
   - pow(elq, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elq, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
   - pow(elq, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elq, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
 }
 for (i = 0; i <= HT; i++) {
  u = VC + VD * (double)(i - HT / 2);
  tf1[i] = 2.0 * (pow(elqc, 0.25) * sin(pi4 * 4.0 * u) - pow(elqc, 2.25) * sin(pi4 * 4.0 * 3.0 * u)
   + pow(elqc, 6.25) * sin(pi4 * 4.0 * 5.0 * u) - pow(elqc, 12.25) * sin(pi4 * 4.0 * 7.0 * u)
   + pow(elqc, 20.25) * sin(pi4 * 4.0 * 9.0 * u) - pow(elqc, 30.25) * sin(pi4 * 4.0 * 11.0 * u)
   + pow(elqc, 42.25) * sin(pi4 * 4.0 * 13.0 * u) - pow(elqc, 56.25) * sin(pi4 * 4.0 * 15.0 * u)
   + pow(elqc, 72.25) * sin(pi4 * 4.0 * 17.0 * u) - pow(elqc, 90.25) * sin(pi4 * 4.0 * 19.0 * u));
  tf2[i] = 2.0 * (pow(elqc, 0.25) * cos(pi4 * 4.0 * u) + pow(elqc, 2.25) * cos(pi4 * 4.0 * 3.0 * u)
   + pow(elqc, 6.25) * cos(pi4 * 4.0 * 5.0 * u) + pow(elqc, 12.25) * cos(pi4 * 4.0 * 7.0 * u)
   + pow(elqc, 20.25) * cos(pi4 * 4.0 * 9.0 * u) + pow(elqc, 30.25) * cos(pi4 * 4.0 * 11.0 * u)
   + pow(elqc, 42.25) * cos(pi4 * 4.0 * 13.0 * u) + pow(elqc, 56.25) * cos(pi4 * 4.0 * 15.0 * u)
   + pow(elqc, 72.25) * cos(pi4 * 4.0 * 17.0 * u) + pow(elqc, 90.25) * cos(pi4 * 4.0 * 19.0 * u));
  tf3[i] = 1.0 + 2.0 * (elqc * cos(pi4 * 8.0 * u) + pow(elqc, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
   + pow(elqc, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elqc, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
   + pow(elqc, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elqc, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
   + pow(elqc, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elqc, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
   + pow(elqc, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elqc, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
  tf0[i] = 1.0 + 2.0 * (-elqc * cos(pi4 * 8.0 * u) + pow(elqc, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
   - pow(elqc, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elqc, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
   - pow(elqc, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elqc, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
   - pow(elqc, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elqc, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
   - pow(elqc, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elqc, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
 }

 // prepare denominator
 for (i = 0; i <= WT; i++) {
  tea[i] = te1[i] * te1[i];
  teb[i] = te0[i] * te0[i];
 }
 for (i = 0; i <= HT; i++) {
  tfa[i] = tf1[i] * tf1[i];
  tfb[i] = tf2[i] * tf2[i];
 }
}

void drawsn1(unsigned char** pp, int i, int j) {
 double sn1, sn2, sn0;
 double z;
 int k;

 sn0 = tea[i] * tfa[j] + teb[i] * tfb[j];
 sn1 = elkc / elk * te1[i] * te0[i] * tf3[j] * tf0[j];
 sn2 = elkc / elk * te2[i] * te3[i] * tf1[j] * tf2[j];
 z = sqrt(sn1 * sn1 + sn2 * sn2);
 if (sn0 < 1e-10) { // pole
  *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
  return;
 }
 if (z < 1e-10) { // zero point
  *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
  return;
 }
 if (z / sn0 / alf < asp) { // south pole area
  *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
  return;
 }
 if (sn0 / z * alf < anp) { // north pole area
  *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
  return;
 }
 if (sn1 * (te1[i + 1] * te0[i + 1] * tf3[j] * tf0[j]) <= 0) { // pure imaginary
  *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
  return;
 }
 if (sn1 * (te1[i] * te0[i] * tf3[j + 1] * tf0[j + 1]) <= 0) { // pure imaginary
  *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
  return;
 }
 if (sn2 * (te2[i + 1] * te3[i + 1] * tf1[j] * tf2[j]) <= 0) { // pure real
  *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
  return;
 }
 if (sn2 * (te2[i] * te3[i] * tf1[j + 1] * tf2[j + 1]) <= 0) { // pure real
  *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
  return;
 }
 z = atan(z / alf / sn0) - pi4; // -pi/4 <= z < pi/4
 if ((z > 0.0) && (fmod(z, DD) <= DE * DD)) { // north edge
  *(*pp)++ = cl0[5][0]; *(*pp)++ = cl0[5][1]; *(*pp)++ = cl0[5][2];
  return;
 }
 if ((z <= 0.0) && (fmod(-z, DD) <= DE * DD)) { // south edge
  *(*pp)++ = cl0[6][0]; *(*pp)++ = cl0[6][1]; *(*pp)++ = cl0[6][2];
  return;
 }
 k = (int)((atan2(sn2, sn1) / pi4) * 2.999 + 12.0); // hard
 *(*pp)++ = cl1[k][0]; *(*pp)++ = cl1[k][1]; *(*pp)++ = cl1[k][2];

}

void drawcn1(unsigned char** pp, int i, int j) {
 double cn1, cn2, cn0;
 double z;
 int k;

 cn0 = tea[i] * tfa[j] + teb[i] * tfb[j];
 cn1 = elkc / elk * te2[i] * te0[i] * tf2[j] * tf0[j];
 cn2 = -elkc / elk * te1[i] * te3[i] * tf1[j] * tf3[j];
 z = sqrt(cn1 * cn1 + cn2 * cn2);
 if (cn0 < 1e-10) { // pole
  *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
  return;
 }
 if (z < 1e-10) { // zero point
  *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
  return;
 }
 if (z / cn0 / alf < asp) { // south pole area
  *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
  return;
 }
 if (cn0 / z * alf < anp) { // north pole area
  *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
  return;
 }
 if (cn1 * (te2[i + 1] * te0[i + 1] * tf2[j] * tf0[j]) <= 0) { // pure imaginary
  *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
  return;
 }
 if (cn1 * (te2[i] * te0[i] * tf2[j + 1] * tf0[j + 1]) <= 0) { // pure imaginary
  *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
  return;
 }
 if (cn2 * (-te1[i + 1] * te3[i + 1] * tf1[j] * tf3[j]) <= 0) { // pure real
  *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
  return;
 }
 if (cn2 * (-te1[i] * te3[i] * tf1[j + 1] * tf3[j + 1]) <= 0) { // pure real
  *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
  return;
 }
 z = atan(z / alf / cn0) - pi4; // -pi/4 <= z < pi/4
 if ((z > 0.0) && (fmod(z, DD) <= DE * DD)) { // north edge
  *(*pp)++ = cl0[5][0]; *(*pp)++ = cl0[5][1]; *(*pp)++ = cl0[5][2];
  return;
 }
 if ((z <= 0.0) && (fmod(-z, DD) <= DE * DD)) { // south edge
  *(*pp)++ = cl0[6][0]; *(*pp)++ = cl0[6][1]; *(*pp)++ = cl0[6][2];
  return;
 }
 k = (int)((atan2(cn2, cn1) / pi4) * 2.999 + 12.0); // hard
 *(*pp)++ = cl1[k][0]; *(*pp)++ = cl1[k][1]; *(*pp)++ = cl1[k][2];
}

void drawdn1(unsigned char** pp, int i, int j) {
 double dn1, dn2, dn0;
 double z;
 int k;

 dn0 = tea[i] * tfa[j] + teb[i] * tfb[j];
 dn1 = elkc * te3[i] * te0[i] * tf2[j] * tf3[j];
 dn2 = -elkc * te1[i] * te2[i] * tf1[j] * tf0[j];
 z = sqrt(dn1 * dn1 + dn2 * dn2);
 if (dn0 < 1e-10) { // pole
  *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
  return;
 }
 if (z < 1e-10) { // zero point
  *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
  return;
 }
 if (z / dn0 / alf < asp) { // south pole area
  *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
  return;
 }
 if (dn0 / z * alf < anp) { // north pole area
  *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
  return;
 }
 if (dn1 * (te3[i + 1] * te0[i + 1] * tf2[j] * tf3[j]) <= 0) { // pure imaginary
  *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
  return;
 }
 if (dn1 * (te3[i] * te0[i] * tf2[j + 1] * tf3[j + 1]) <= 0) { // pure imaginary
  *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
  return;
 }
 if (dn2 * (-te1[i + 1] * te2[i + 1] * tf1[j] * tf0[j]) <= 0) { // pure real
  *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
  return;
 }
 if (dn2 * (-te1[i] * te2[i] * tf1[j + 1] * tf0[j + 1]) <= 0) { // pure real
  *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
  return;
 }
 z = atan(z / alf / dn0) - pi4; // -pi/4 <= z < pi/4
 if ((z > 0.0) && (fmod(z, DD) <= DE * DD)) { // north edge
  *(*pp)++ = cl0[5][0]; *(*pp)++ = cl0[5][1]; *(*pp)++ = cl0[5][2];
  return;
 }
 if ((z <= 0.0) && (fmod(-z, DD) <= DE * DD)) { // south edge
  *(*pp)++ = cl0[6][0]; *(*pp)++ = cl0[6][1]; *(*pp)++ = cl0[6][2];
  return;
 }
 k = (int)((atan2(dn2, dn1) / pi4) * 2.999 + 12.0); // hard
 *(*pp)++ = cl1[k][0]; *(*pp)++ = cl1[k][1]; *(*pp)++ = cl1[k][2];

}

void draw1(void) {
 int i, j;
 unsigned char* p2, *p3, *p4;

 p2 = b2; p3 = b3; p4 = b4;
 for (j = 0; j < HT; j++) {
  for (i = 0; i < WT; i++) {
   drawsn1(&p2, i, j);
   drawcn1(&p3, i, j);
   drawdn1(&p4, i, j);
  }
 }
}

int main()
{
 FILE* ftp;
 errno_t err;
 unsigned int ui;
 double z;

 // writing hedder
 b1[0x00] = 'B'; b1[0x01] = 'M';
 b1[0x02] = FZ & 0xff; b1[0x03] = FZ >> 8 & 0xff; b1[0x04] = FZ >> 16 & 0xff; b1[0x05] = FZ >> 24;
 b1[0x06] = 0; b1[0x07] = 0;  b1[0x08] = 0; b1[0x09] = 0;
 b1[0x0a] = 0x36; b1[0x0b] = 0; b1[0x0c] = 0; b1[0x0d] = 0;
 b1[0x0e] = 0x28; b1[0x0f] = 0; b1[0x10] = 0; b1[0x11] = 0;
 b1[0x12] = WT & 0xff; b1[0x13] = WT >> 8 & 0xff; b1[0x14] = WT >> 16 & 0xff; b1[0x15] = WT >> 24;
 b1[0x16] = HT & 0xff; b1[0x17] = HT >> 8 & 0xff; b1[0x18] = HT >> 16 & 0xff; b1[0x19] = HT >> 24;
 b1[0x1a] = 0x1; b1[0x1b] = 0;  b1[0x1c] = 0x18; b1[0x1d] = 0;
 b1[0x1e] = 0; b1[0x1f] = 0; b1[0x20] = 0; b1[0x21] = 0;
 b1[0x22] = BZ & 0xff; b1[0x23] = BZ >> 8 & 0xff; b1[0x24] = BZ >> 16 & 0xff; b1[0x25] = BZ >> 24;
 b1[0x26] = 0; b1[0x27] = 0; b1[0x28] = 0; b1[0x29] = 0;
 b1[0x2a] = 0; b1[0x2b] = 0; b1[0x2c] = 0; b1[0x2d] = 0;
 b1[0x2e] = 0; b1[0x2f] = 0; b1[0x30] = 0; b1[0x31] = 0;
 b1[0x32] = 0; b1[0x33] = 0; b1[0x34] = 0; b1[0x35] = 0;

 // initialize
 init1();

 // making bitmap data
 draw1();

 // output file
 err = fopen_s(&ftp, "sn.bmp", "wb");
 if (err != NULL) return err;
 for (ui = 0; ui < 0x36; ui++) {
  fputc(b1[ui], ftp);
 }
 for (ui = 0; ui < BZ; ui++) {
  fputc(b2[ui], ftp);
 }
 err = fclose(ftp);
 if (err != NULL) return err;
 err = fopen_s(&ftp, "cn.bmp", "wb");
 if (err != NULL) return err;
 for (ui = 0; ui < 0x36; ui++) {
  fputc(b1[ui], ftp);
 }
 for (ui = 0; ui < BZ; ui++) {
  fputc(b3[ui], ftp);
 }
 err = fclose(ftp);
 if (err != NULL) return err;
 err = fopen_s(&ftp, "dn.bmp", "wb");
 if (err != NULL) return err;
 for (ui = 0; ui < 0x36; ui++) {
  fputc(b1[ui], ftp);
 }
 for (ui = 0; ui < BZ; ui++) {
  fputc(b4[ui], ftp);
 }
 err = fclose(ftp);
 if (err != NULL) return err;
 err = fopen_s(&ftp, "el.txt", "w");
 if (err != NULL) return err;
 fprintf(ftp, "Jacobi elliptic function constants\n");
 fprintf(ftp, "k = %10.5g, k^2 = %10.5g, k' = %10.5g, k'^2 = %10.5g\n", elk, elk * elk, elkc, elkc * elkc);
 fprintf(ftp, "1/k = %10.5g, k'/k = %10.5g\n", 1 / elk, elkc / elk);
 fprintf(ftp, "K = %10.5g, K' = %10.5g, K/K' = %10.5g\n", ellk, ellkc, ellk / ellkc);
 fprintf(ftp, "K/(pi/2) = %10.5g, K'/(pi/2) = %10.5g\n", ellk / pi4 / 2.0, ellkc / pi4 / 2.0);
 fprintf(ftp, "q  = %10.5g, te2(0) = %10.5g, te3(0) = %10.5g, te0(0) = %10.5g\n", elq, te2a, te3a, te0a);
 fprintf(ftp, "q' = %10.5g, tf2(0) = %10.5g, tf3(0) = %10.5g, tf0(0) = %10.5g\n", elqc, tf2a, tf3a, tf0a);
 fprintf(ftp, "horizontal: WT = %6i, CT = %10.5g, D = %10.5g\n", WT, UC, UD);
 fprintf(ftp, "vertical:   HT = %6i, CT = %10.5g, D = %10.5g\n", WT, VC, VD);
 fprintf(ftp, "northern hemisphere parallels: (degree) value\n");
 for (z = 0; z < pi4 - 0.001; z += DD) {
  fprintf(ftp, "(%5.2f) %10.5g\n", z / pi4 * 90.0, tan(pi4 + z) * alf);
 }
 fprintf(ftp, "north polar (infinity) area = %10.5g(deg)\n", atan(anp) / pi4 * 90.0);
 fprintf(ftp, "southern hemisphere parallels: (degree) value\n");
 for (z = 0; z < pi4 - 0.001; z += DD) {
  fprintf(ftp, "(%5.2f) %10.5g\n", z / pi4 * 90.0, tan(pi4 - z) * alf);
 }
 fprintf(ftp, "south polar (zero point) area = %10.5g(deg)\n", atan(asp) / pi4 *90.0);
 err = fclose(ftp);
 if (err != NULL) return err;
 printf("complete\n");
 return err;
}

　上記のプログラムと本ブログに掲載の関連図の著作権は私にあります。しかし、改造などはご自由になさってください。明らかな改造後のプログラムの著作権はあなたのものとなります。係数などを必要に応じて変更した、プロクラムが出力するグラフもあなたのものです。私に連絡する必要はありませんし、許諾なども必要ありません。

4145. 楕円関数、複素平面、上級編

　ヤコビの楕円関数の複素平面表示の改良版が一応の完成を見たと思いますので公表します。
　商品や論文ではありませんので徹底的な検証はしておらず、ですからソースプログラムを開示するのです。質問などには答える時間が取れないと思います。今から必要十分な説明をするつもりですので、それでご容赦下さい。
　わずかな改良ですが、多少ややこしくて、上級者編になってしまいました。そのため、初期のプログラム(4137.)はそのままにしておきます。

　楕円関数の母数k = √(1/2)の場合の図の印象は前回(4136.)とほとんど変わりません。
　最初の3図はk = √(0.9755)のものです。この不思議な数値は、「岩波 数学公式III」の264ページの表に載っている数値です。k = 1に近い典型値として採用しました。
　まずは基本となるsn(u, √(0.9755))の図です。
(sn)

　画面中央がガウス平面(複素平面/複素数平面)の中心(0+0i)の場所です。青系の灰色の楕円が収束しているように見えるのは、関数値の複素数の絶対値が0に近づいていることを示し、中央には少し明るい空色のドットが見えるはずです。ここが関数値0、つまり零点です。
　虹色の帯は関数値の複素数の偏角を示します。中央から右方向の紅色が0度付近で、白線は純実数(虚数部=0)の線です。中央から上には黄色の領域があって、90度付近です。黒線は純虚数(実数部=0)の線です。ヤコビの楕円関数では純実数と純虚数の線は水平線と垂直線になるようです。
　中央から左の草色は偏角が180度付近で、白線がありますから実数が負に進んでいることが分かります。中央から下の青色は270度付近です。

　再び、中央から右に進んで行くと、白線の十字路があります。ここが周期(K)、つまり定義域の(3.2547+0i)の位置です。sn(u)の値はちょうど1で、背景の青系の灰色と赤系の灰色がクロスしています。この灰色の領域は本来は複素数の絶対値の等高線なので線なのですが、ぐぐっと延びて面状になってしまっています。これは、kが1に近いsn(u)の値が±1付近は頭打ちとなって平らに近い、全体としてみると方形波に近づいていることを示します。

　そのままさらに右に進むと、再び関数値が減って零点が現れます。ここが周期の(2K)の場所です。水平線の反対側に、同じ模様の所があって、こちらも(2K)です。その右に関数値-1の(3K)があって、原点に戻ります。つまりsn(u)の実数方向の周期は4Kです。

　今度は中央から上に進むと、赤系の灰色の楕円が集中している場所があります。中央には極(無限遠点)を示す少し明るい土色のドットが見えるはずです。ここは関数の値の分母が0になる点(の近傍)です。ここが第2周期(K')、つまり定義域の(0+1.5806i)の場所です。
　その直上に再び青系の灰色の楕円が集中している場所があって、中央付近と同じ図柄になっています。つまり、虚数方向の周期は2K'です。4Kと2K'、これが楕円関数の二重周期と呼ばれる現象です。前回はk = √(0.5)だったのでK = K'でしたが、今回は異なっていて、結果として格子が正方形では無く長方形になっています。

　cn(u, √(0.9755))の図に移ります。
(cn)

　前回のcn(u, √(0.5))とずいぶん印象が異なります。しかし、周期性は同じで、横方向、つまり実軸方向が4Kで、虚数方向は斜め上下の(2K+2iK')であるのは同じです。また、k=1に近いcnは実軸に沿って奇数Kの場所で0となり、そのあたりはなだらかで、偶数Kの所でスパイク状に立ち上がって±1になる様子が、かなり想像の補強が必要ですが、見えると思います。極のドットはsnと同じ位置に打たれていますが、かなり小さいです。

　dn(u, √(0.9755))の図です。
(dn)

　(K)の位置(3.2547+0i)の関数値は0にはなりませんが、0に近づいていることは分かると思います。極の位置はこれも同じです。周期は2Kと4iK'です。基本周期平行四辺形(基本領域)が正方形に近い感じに見えますが、これはK/K' = 2.0592と2に近いためで、偶然です。

　以下は、以上の図と同時に出力されるテキストファイルです。楕円関数に関わる定数と、グラフの中心とドット間のピッチ、関数値の等高線(外縁)の値です。詳しい説明はプログラムの解説時にする予定です。

Jacobi elliptic function constants
k =    0.98767, k^2 =     0.9755, k' =    0.15652, k'^2 =     0.0245
1/k =     1.0125, k'/k =    0.15848
K =     3.2547, K' =     1.5806, K/K' =     2.0592
K/(pi/2) =      2.072, K'/(pi/2) =     1.0062
q  =    0.21749, te2(0) =     1.4306, te3(0) =     1.4395, te0(0) =    0.56949
q' =  0.0015503, tf2(0) =    0.39686, tf3(0) =     1.0031, tf0(0) =     0.9969
horizontal: WT =   1024, CT =    0.00125, D =     0.0025
vertical:   HT =   1024, CT =    0.00125, D =  0.0051481
northern hemisphere parallels: (degree) value
( 0.00)          1
(15.00)     1.3032
(30.00)     1.7321
(45.00)     2.4142
(60.00)     3.7321
(75.00)     7.5958
north polar (infinity) area =          4(deg)
southern hemisphere parallels: (degree) value
( 0.00)          1
(15.00)    0.76733
(30.00)    0.57735
(45.00)    0.41421
(60.00)    0.26795
(75.00)    0.13165
south polar (zero point) area =          4(deg)

　参考に、母数を補母数と入れ替えた、k=√(0.0245)のsn, cn, dnのグラフを掲載します。上図との類似性にご注目下さい。
(sn)

(cn)

(dn)

Jacobi elliptic function constants
k =    0.15652, k^2 =     0.0245, k' =    0.98767, k'^2 =     0.9755
1/k =     6.3888, k'/k =       6.31
K =     1.5806, K' =     3.2547, K/K' =    0.48562
K/(pi/2) =     1.0062, K'/(pi/2) =      2.072
q  =  0.0015503, te2(0) =    0.39686, te3(0) =     1.0031, te0(0) =     0.9969
q' =    0.21749, tf2(0) =     1.4306, tf3(0) =     1.4395, tf0(0) =    0.56949
horizontal: WT =   1024, CT =    0.00125, D =  0.0051481
vertical:   HT =   1024, CT =    0.00125, D =     0.0025
northern hemisphere parallels: (degree) value
( 0.00)          1
(15.00)     1.3032
(30.00)     1.7321
(45.00)     2.4142
(60.00)     3.7321
(75.00)     7.5958
north polar (infinity) area =          4(deg)
southern hemisphere parallels: (degree) value
( 0.00)          1
(15.00)    0.76733
(30.00)    0.57735
(45.00)    0.41421
(60.00)    0.26795
(75.00)    0.13165
south polar (zero point) area =          4(deg)

4144. 土曜日

　複素平面の方のヤコビの楕円関数グラフの改良は途中まで進んでいて、多分、明日にも公開できると思います。
　来月は例年通り忙しくなりそうなので、何とかまとめないといけません。

　国内は相変わらずで、特に変化無し。先週の出張先はおそらく製造業に属する企業(複数)と思います。案内された部屋は小綺麗になっていて、緊張感はあるものの明るい雰囲気で、想像するに笑いが止まらない状態なのだと思います。こちらは他企業なので尋ねたりはしません。仕事の範囲で拝見するだけです。

　2週間ほど前からネットの論調が明らかに変化していて、中国の経済的な落ち込みが話題になっています。少し前までは中国に不利な論調は避けられていました。もっぱらの噂では、米中対立が新たな段階に入ったから、だそうです。とはいえ、東アジア情勢の情報不足は続いていて、外部からはよく分かりません。
　しばらくしたら落ち着くだろう、と思っていたら、ますます対立が先鋭化しているそうです。このままでは済まなそうで、一体どこまで行くのやら。ただ、我が政府はやたらと落ち着いているので、とっくに複数のシナリオは十分に検討されているのだと思います。

4143. キャラの誕生日、他

　明日、7月29日はPS4の最新アイマスゲーム、スターリットシーズンに出てくる仮想アイドル、小宮果穂(シャイニーカラーズ)の誕生日だそうです。
　この日は同仮想アイドル、田中摩美々(シャイニーカラーズ)々の声優、菅沼千紗さんの誕生日でもあります。

　翌、7月30日は同仮想アイドル、伊吹翼(ミリオンライブ)の誕生日だそうです。

　いつものようにスターリットシーズンのPV新着欄で有志・精鋭Pがお祝いのPVを上げるはずです。

4142. 木曜日

　本日は午前の半日出張でした。激務ではありませんが、作業量が多くて帰社した頃にはへとへと。内部の仕事をやっつけた以外は何も出来ず。
　今は趣味の数学も進まないので、かえって暴走気味の思考が止まって冷却するチャンスと思っておきます。

　ノートを繰ってやり残したことが無いかどうかをチェック、と。楕円関数に関してはまだ手も足も出ない領域がたくさん残っています。これを進めると時間が取られそうなので、どうしようかな。公式集を見て、すぐに分かりそうな数式などが出てきたら良いのですが。
　楕円関数自体は、三角関数や双曲線関数程度には身近なところで観察できるようです。なぜかというと、線形では無いものの比較的に簡単な(常)微分方程式の解だから。そして、これに近縁の有用な特殊関数はいくつかあって、出来ればそちらのグラフ表示も考えたいと思っていて。なあに、ぱっと見は高校の物理や化学で出てきた図表と似ています。

　こうした理由で、楕円関数は20世紀初頭までの数学らしい数学の時代と、今の高度な数学と機能的な数学の2方向に分かれた感じの数学を結びつける話題として、かなり良い線の話題と思います。私としても、なぜ現代数学がこのような感じなのか、少しは理解が進んだような気がしています。

4141. 楕円関数、複素平面編、続き^4

　現在、3点の改良を考えていて、どれもそれほど時間を取らないはずですが、一定の検証が必要ですので今度の週末の作業になると思います。出てくる図の印象は全く変わりません。

　1つは南極領域(零点)と北極領域(極)のキャップの青ドット・赤ドットの厳密化。今のは実は取って付けです。北極から半径100kmは北極領域とかそんな感じの定義にしたいです。

　2つめはK'/Kで縦横比の修正をすること。現在は1:1のままです。実数編では我慢できましたが、2次元では無いと誤解を招く恐れがあるからです。その内容は実装して、検証してから言います。

　3つめは等高線の南半球と北半球の境を移動可とすること。私はリーマン球面の大きさ調整と頭の中では理解していますが、こちらも誤解の無いように説明するのは難しそうです。効果は一目で分かるはずで、私の頭の中ではすっかり完成しているその効果が現実の絵として現れたら、多分、自慢げに説明すると思います。

　いずれも、一般の複素有理関数の表示に役立つはずです。つまり先行して、楕円関数を例題としての導入です。なんだか、読者の皆様に導かれいてるような気がします。ありがとうございます。

4140. 心白の誕生日

　明日、7月26日はPS4の最新アイマスゲーム、スターリットシーズンに出てくる仮想アイドル、奥空心白(新人)の誕生日だそうです。
　また、この日は最初のアイドルマスター、アーケード版(2005年)の稼働日なのでアイマス自体の誕生日とされています。
　いつものようにスターリットシーズンのPV新着欄で有志Pがお祝いのPVを上げるはずです。

4139. 楕円関数、複素平面編、続き^3

　プログラムでは複素数型を使っていません。極(無限大)の表現を探っている間にどこかに行ってしまいました。
　複素数の分数式の分母の虚数は通分の要領で分子に集めることが出来て、その結果、定義域が複素平面の複素関数の値は、

　f(x + iy) = (a + bi) / d

　の形になります。dがゼロの点が極です。現在、ここをもう少し洗練された形で表示する方法を探っています。つまり今回の楕円関数のプログラミングは私にとっては、とても役立ったということ。

　偶然なのか、ヤコビの楕円関数の加法定理の公式が上記の形になっています。uとvが実数として、

　sn(u + iv) = [sn(u) cn(iv) dn(iv) + cn(u) dn(u) sn(iv)] / [1 - k^2 (sn(u))^2 (sn(iv))^2]

　です。これにヤコビの虚変換の式、

　sn(iv, k) = i sn(v, k') / cn(v, k')
　cn(iv) = 1 / cn(v, k')
　dn(iv) = dn(v, k') / cn(v, k')

　を適用します。kは楕円関数の母数、k'は補母数で、　

　sn(v, k') = sn(v, √(1-k^2))

　の事です。上述のsn(u + iv)は母数を書くと、sn(u + iv, k)です。

　加算公式を見ると、分母が偶然なのか実数になっていて、

　1 + k^2 (sn(u, k))^2 (sn(v, k'))^2 / (cn(v, k'))^2

　となり、分子は前半が、

　sn(u, k) (1 / cn(v, k')) (dn(v, k') / cn(v, k'))

　と実数になり、後半が、

　i cn(u, k) dn(u, k) sn(v, k') / cn(v, k')

　と純虚数になります。つまり、複素数のsn(u + iv)は、

　sn(u, k), cn(u, k), dn(u, k), sn(v, k'), cn(v, k'), dn(v, k')

　の6関数から算出できます。これで複素変数のヤコビの楕円関数が計算できますが、楕円テータ関数を使う方針なので、さらに、

　sn(u) = [θ30 θ1(w)] / [θ20 θ0(w)]
　cn(u) = [θ00 θ2(w)] / [θ20 θ0(w)]
　dn(u) = [θ00 θ3(w)] / [θ30 θ0(w)]

　を適用します。ただし、w = u/(2K)で、楕円関数の周期を表すKで換算したものです。

　やや長くなったので、続きます。

4138. 楕円関数、複素平面編、続き^2

　もう少し何とか出来るだろう、の感じのプログラムですが、扱える関数がほぼヤコビの楕円関数に限られているので、これ以上の情熱が続きそうにありません。他の一般の複素関数にも役立つ形を考えた方が良いような気がするので、これはこれとして固定させていただきます。

　コンパイルにはコンパイラが必要です。
　私が使っているのはMicrosoft Visual Studioで、無料版がオンラインで手に入るはずです。C言語のコンソールプログラムが作れるようにダウンロードしてインストールします。コンソール画面でHello, worldが表示されるプログラムの新規作成を選び、表示されるソースプログラムを前回掲示した楕円関数のプログラムに置き換えます。
　メニュー直下の「ローカルWindowsデバッガー」のボタン(?)をクリックすると、ソースプログラムがセーブされ、コンパイルが実施され、実行されてコンソール画面が現れます。この時点で、グラフのファイル、sn.bmp, cn.bmp, dn.bmpがソースプログラムのフォルダーに書かれているはずです。適当な表示ソフトで表示します。通常はクリックするだけで絵が出てきます。
　さらに、その時の定数を書いたテキストファイル、el.txtが書かれます。
　もう一度実行するとこれらの4ファイルは上書きされますので、後で見たい場合は適当なフォルダーに移動またはコピーしておきます。

　k値を変更するには、行23の、

　double elk = sqrt(0.4); // k^2: ***modify here

　の0.4の所を書き換えます。0≦ k ≦1ですが、0と1に極端に近い場所は誤差が出ます。
　この0.4はk値の自乗です。自乗で指定する方が便利なのでこうしました。もちろん、k値をそのまま書いてもOKです。
　q値(0≦ q ≦1)から出発することも出来て、その場合はこの行をコメントアウトし、行78と79もコメントアウトし、その下の、

    // elq = 0.04321; // starting with q

　のコメントを外して実行できるようにします。

　本プログラムでは実際上、k値(またはq値)を変更することしか出来ることはありません。見所は沢山あるので、参考書を見ながらあれこれ試すと楽しめると思います。

　グラフの縦横のドット数を変更するには、最初の方にあるWTとHTの値を変更します。
　中心座標を変更するにはその下にあるUCとVCを変更します。拡大・縮小はドット間の差分値を指定するUDとVDの値を変更します。
　等高線の間隔を決めているのはその下のDDです。プロクラムの値では、地球儀の緯度の15°間隔になっていて、極地までに6本の等高線が見られます。DEは等高線の幅です。数値が小さいと細くなります。
　南極と北極付近は塗りつぶしています。その大きさを決めるのがanpとaspです。この値は微妙なので、グラフを見ながら適当に指定するしかありません。

　配列cl0[7][3]は、等高線などの色の指定です。青・緑・赤の順で0-255の値で指定します。最初が背景色ですが、普通は表示されません。cl0[1]は北極(極。無限遠点)のドットの色、cl0[2]は南極(零点)の色です。cl0[3]は純虚数の、cl0[4]は純実数の線の色です。cl0[5]は北半球(関数値1～∞)の等高線の、cl0[6]は南半球(関数値0～1)の等高線の色です。
　配列cl1[24][3]は関数値の偏角の色相です。24個ありますが、12色で指定しています。
　これらの色表現は私が適当に決めただけですので、見やすいように調整してみて下さい。

　最後に書いてあるmain()関数はbmpファイルの作成と出力が主な作業で、グラフの描画は関数init1()と関数draw1()がやっています。その解説は稿を改めます。

4137. 楕円関数、複素平面編、続き

　C言語で書かれたヤコビの楕円関数の複素平面表示のプログラムです。実行すると、sn(u, k)、cn(u, k)、dn(u, k)の3関数のグラフと、その時の各種パラメータの値のテキストファイル()が出力されます。グラフは標準的なbmpファイルで出力されるので、適当な表示ソフト(ペイントなど)で表示させます。

// ConsoleApplication6.cpp
// ellptic2.cpp; Jacobi elliptic functions on a Gaussian plane
// based on bmp24.cpp
//  200828; making 24bit color bmp file
// 230711; 230722

#include <stdio.h>
#include <math.h>

const int WT = 1024; // width  (mod 8 = 0) *** modify here
const int HT = 1024; // height (mod 8 = 0) *** modify here
const int BZ = WT * HT * 3; // buffer size
const int FZ = BZ + 0x36; // file size

unsigned char b1[0x36];  // header*2
unsigned char b2[BZ];  // bitmap sn
unsigned char b3[BZ];  // bitmap cn
unsigned char b4[BZ];  // bitmap dn

const double pi4 = atan(1.0);

// elliptic k
double elk = sqrt(0.4); // k^2: ***modify here

// draw area definition *** modify here
const double UC = 0.00125; // real base
const double UD = 0.0025; // real pitch // hard
const double VC = 0.00125; // imaginary base
const double VD = 0.0025; // imaginary pitch // hard

// contour information *** modify here
const double DD = pi4 * 15.0 /90.0 ; // value hight pitch
const double DE = 0.125; // value edge width (0.01 - 0.99)

// pole area definition
const double anp = 1e-3; // pole area (infinity)
const double asp = 1e-3; // zero point area

// color index
const unsigned char cl0[7][3] = { // {b, g, r} 0-255
    {0x2e, 0x2e, 0x2e}, // background
    {0x66, 0x66, 0xbb}, // pole (infinity)
    {0xbb, 0xbb, 0x22}, // zero point
    {0x11, 0x11, 0x11}, // real part = 0; pure imaginary
    {0xff, 0xff, 0xff}, // imaginary part = 0; pure real
    {0x66, 0x66, 0xaa}, // contour edge nothern hemisphere
    {0x99, 0x99, 0x66} // contour edge southern hemisphere
};

const unsigned char cl1[24][3] = {
    {0x99, 0xff, 0x33}, {0xff, 0xff, 0x33}, {0xff, 0xff, 0x33}, {0xff, 0x99, 0x33},
    {0xff, 0x99, 0x33}, {0xff, 0x33, 0x33}, {0xff, 0x33, 0x33}, {0xff, 0x33, 0x99},
    {0xff, 0x33, 0x99}, {0xff, 0x33, 0xff}, {0xff, 0x33, 0xff}, {0x99, 0x33, 0xff},
    {0x99, 0x33, 0xff}, {0x33, 0x33, 0xff}, {0x33, 0x33, 0xff}, {0x33, 0x99, 0xff},
    {0x33, 0x99, 0xff}, {0x33, 0xff, 0xff}, {0x33, 0xff, 0xff}, {0x33, 0xff, 0x99},
    {0x33, 0xff, 0x99}, {0x33, 0xff, 0x33}, {0x33, 0xff, 0x33}, {0x99, 0xff, 0x33}
};

// elliptic function constants

double elq, elqc; // k, q, q'
double elkc, ellk, ellkc; // k', K, K'
double te2a, te3a, te0a, tf2a, tf3a, tf0a; // theta(0)s

// elliptic function lines
double te1[WT + 1], te2[WT + 1], te3[WT + 1], te0[WT + 1]; // using q
double tf1[HT + 1], tf2[HT + 1], tf3[HT + 1], tf0[HT + 1]; // using q'

// elliptic function buffers
double tea[WT + 1], teb[WT + 1], tfa[HT + 1], tfb[HT + 1]; // denominator: sn0 = cn0 = dn0

void init1(void) {
    int i;
    double u;
    

    // prepare constants
    double ell = (1.0 - pow(1.0 - elk * elk, 0.25)) / 2.0 / (1.0 + pow(1.0 - elk * elk, 0.25));
    elq = ell + 2.0 * pow(ell, 5.0) + 15.0 * pow(ell, 9.0) + 150.0 * pow(ell, 13.0) + 1707.0 * pow(ell, 17.0);
    // elq = 0.04321; // starting with q
    te2a = 2.0 * (pow(elq, 0.25) + pow(elq, 2.25) + pow(elq, 6.25) + pow(elq, 12.25) + pow(elq, 20.25)
        + pow(elq, 30.25) + pow(elq, 42.25) + pow(elq, 56.25) + pow(elq, 72.25) + pow(elq, 90.25));
    te3a = 1.0 + 2.0 * (elq + pow(elq, 4.0) + pow(elq, 9.0) + pow(elq, 16.0) + pow(elq, 25.0)
        + pow(elq, 36.0) + pow(elq, 49.0) + pow(elq, 64.0) + pow(elq, 81.0) + pow(elq, 100.0));
    te0a = 1.0 + 2.0 * (-elq + pow(elq, 4.0) - pow(elq, 9.0) + pow(elq, 16.0) - pow(elq, 25.0)
        + pow(elq, 36.0) - pow(elq, 49.0) + pow(elq, 64.0) - pow(elq, 81.0) + pow(elq, 100.0));
    elqc = exp(pi4 * pi4 * 16.0 / log(elq));
    tf2a = 2.0 * (pow(elqc, 0.25) + pow(elqc, 2.25) + pow(elqc, 6.25) + pow(elqc, 12.25) + pow(elqc, 20.25)
        + pow(elqc, 30.25) + pow(elqc, 42.25) + pow(elqc, 56.25) + pow(elqc, 72.25) + pow(elqc, 90.25));
    tf3a = 1.0 + 2.0 * (elqc + pow(elqc, 4.0) + pow(elqc, 9.0) + pow(elqc, 16.0) + pow(elqc, 25.0)
        + pow(elqc, 36.0) + pow(elqc, 49.0) + pow(elqc, 64.0) + pow(elqc, 81.0) + pow(elqc, 100.0));
    tf0a = 1.0 + 2.0 * (-elqc + pow(elqc, 4.0) - pow(elqc, 9.0) + pow(elqc, 16.0) - pow(elqc, 25.0)
        + pow(elqc, 36.0) - pow(elqc, 49.0) + pow(elqc, 64.0) - pow(elqc, 81.0) + pow(elqc, 100.0));
    elk = pow(te2a / te3a, 2.0); elkc = pow(te0a / te3a, 2.0);
    ellk = 2.0 * pi4 * pow(te3a, 2.0); ellkc = ellk / (4.0 * pi4) * log(1 / elq);

    // prepare lines
    for (i = 0; i <= WT; i++) {
        u = UC + UD * (double)(i - WT / 2);
        te1[i] = 2 * (pow(elq, 0.25) * sin(pi4 * 4.0 * u) - pow(elq, 2.25) * sin(pi4 * 4.0 * 3.0 * u)
            + pow(elq, 6.25) * sin(pi4 * 4.0 * 5.0 * u) - pow(elq, 12.25) * sin(pi4 * 4.0 * 7.0 * u)
            + pow(elq, 20.25) * sin(pi4 * 4.0 * 9.0 * u) - pow(elq, 30.25) * sin(pi4 * 4.0 * 11.0 * u)
            + pow(elq, 42.25) * sin(pi4 * 4.0 * 13.0 * u) - pow(elq, 56.25) * sin(pi4 * 4.0 * 15.0 * u)
            + pow(elq, 72.25) * sin(pi4 * 4.0 * 17.0 * u) - pow(elq, 90.25) * sin(pi4 * 4.0 * 19.0 * u));
        te2[i] = 2 * (pow(elq, 0.25) * cos(pi4 * 4.0 * u) + pow(elq, 2.25) * cos(pi4 * 4.0 * 3.0 * u)
            + pow(elq, 6.25) * cos(pi4 * 4.0 * 5.0 * u) + pow(elq, 12.25) * cos(pi4 * 4.0 * 7.0 * u)
            + pow(elq, 20.25) * cos(pi4 * 4.0 * 9.0 * u) + pow(elq, 30.25) * cos(pi4 * 4.0 * 11.0 * u)
            + pow(elq, 42.25) * cos(pi4 * 4.0 * 13.0 * u) + pow(elq, 56.25) * cos(pi4 * 4.0 * 15.0 * u)
            + pow(elq, 72.25) * cos(pi4 * 4.0 * 17.0 * u) + pow(elq, 90.25) * cos(pi4 * 4.0 * 19.0 * u));
        te3[i] = 1.0 + 2.0 * (elq * cos(pi4 * 8.0 * u) + pow(elq, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
            + pow(elq, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elq, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
            + pow(elq, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elq, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
            + pow(elq, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elq, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
            + pow(elq, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elq, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
        te0[i] = 1.0 + 2.0 * (-elq * cos(pi4 * 8.0  * u) + pow(elq, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
            - pow(elq, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elq, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
            - pow(elq, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elq, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
            - pow(elq, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elq, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
            - pow(elq, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elq, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
    }
    for (i = 0; i <= HT; i++) {
        u = VC + VD * (double)(i - HT / 2);
        tf1[i] = 2.0 * (pow(elqc, 0.25) * sin(pi4 * 4.0 * u) - pow(elqc, 2.25) * sin(pi4 * 4.0 * 3.0 * u)
            + pow(elqc, 6.25) * sin(pi4 * 4.0 * 5.0 * u) - pow(elqc, 12.25) * sin(pi4 * 4.0 * 7.0 * u)
            + pow(elqc, 20.25) * sin(pi4 * 4.0 * 9.0 * u) - pow(elqc, 30.25) * sin(pi4 * 4.0 * 11.0 * u)
            + pow(elqc, 42.25) * sin(pi4 * 4.0 * 13.0 * u) - pow(elqc, 56.25) * sin(pi4 * 4.0 * 15.0 * u)
            + pow(elqc, 72.25) * sin(pi4 * 4.0 * 17.0 * u) - pow(elqc, 90.25) * sin(pi4 * 4.0 * 19.0 * u));
        tf2[i] = 2.0 * (pow(elqc, 0.25) * cos(pi4 * 4.0 * u) + pow(elqc, 2.25) * cos(pi4 * 4.0 * 3.0 * u)
            + pow(elqc, 6.25) * cos(pi4 * 4.0 * 5.0 * u) + pow(elqc, 12.25) * cos(pi4 * 4.0 * 7.0 * u)
            + pow(elqc, 20.25) * cos(pi4 * 4.0 * 9.0 * u) + pow(elqc, 30.25) * cos(pi4 * 4.0 * 11.0 * u)
            + pow(elqc, 42.25) * cos(pi4 * 4.0 * 13.0 * u) + pow(elqc, 56.25) * cos(pi4 * 4.0 * 15.0 * u)
            + pow(elqc, 72.25) * cos(pi4 * 4.0 * 17.0 * u) + pow(elqc, 90.25) * cos(pi4 * 4.0 * 19.0 * u));
        tf3[i] = 1.0 + 2.0 * (elqc * cos(pi4 * 8.0 * u) + pow(elqc, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
            + pow(elqc, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elqc, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
            + pow(elqc, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elqc, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
            + pow(elqc, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elqc, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
            + pow(elqc, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elqc, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
        tf0[i] = 1.0 + 2.0 * (-elqc * cos(pi4 * 8.0 * u) + pow(elqc, 4.0) * cos(pi4 * 8.0 * 2.0 * u)
            - pow(elqc, 9.0) * cos(pi4 * 8.0 * 3.0 * u) + pow(elqc, 16.0) * cos(pi4 * 8.0 * 4.0 * u)
            - pow(elqc, 25.0) * cos(pi4 * 8.0 * 5.0 * u) + pow(elqc, 36.0) * cos(pi4 * 8.0 * 6.0 * u)
            - pow(elqc, 49.0) * cos(pi4 * 8.0 * 7.0 * u) + pow(elqc, 64.0) * cos(pi4 * 8.0 * 8.0 * u)
            - pow(elqc, 81.0) * cos(pi4 * 8.0 * 9.0 * u) + pow(elqc, 100.0) * cos(pi4 * 8.0 * 10.0 * u));
    }
    
    // prepare denominator
    for (i = 0; i <= WT; i++) {
        tea[i] = te1[i] * te1[i];
        teb[i] = te0[i] * te0[i];
    }
    for (i = 0; i <= HT; i++) {
        tfa[i] = tf1[i] * tf1[i];
        tfb[i] = tf2[i] * tf2[i];
    }
}

void drawsn1(unsigned char **pp, int i, int j) {
    double sn1, sn2, sn0;
    double z;
    int k;

    sn0 = tea[i] * tfa[j] + teb[i] * tfb[j];
    if (sn0 < anp) { // pole
        *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
        return;
    }
    sn1 = elkc / elk * te1[i] * te0[i] * tf3[j] * tf0[j];
    sn2 = elkc / elk * te2[i] * te3[i] * tf1[j] * tf2[j];
    if ((sn1 * sn1 + sn2 * sn2) / sn0 / sn0 < asp) { // zero point
        *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
        return;
    }
    if (sn1 * (te1[i + 1] * te0[i + 1] * tf3[j] * tf0[j]) <= 0) { // pure imaginary
        *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
        return;
    }
    if (sn1 * (te1[i] * te0[i] * tf3[j + 1] * tf0[j + 1]) <= 0) { // pure imaginary
        *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
        return;
    }
    if (sn2 * (te2[i + 1] * te3[i + 1] * tf1[j] * tf2[j]) <= 0) { // pure real
        *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
        return;
    }
    if (sn2 * (te2[i] * te3[i] * tf1[j + 1] * tf2[j + 1]) <= 0) { // pure real
        *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
        return;
    }
    z = atan(sqrt(sn1 * sn1 + sn2 * sn2) / sn0) - pi4; // -pi/4 <= z < pi/4
    if ((z > 0.0) && (fmod(z, DD) <= DE * DD)) { // north edge
        *(*pp)++ = cl0[5][0]; *(*pp)++ = cl0[5][1]; *(*pp)++ = cl0[5][2];
        return;
    }
    if ((z <= 0.0) && (fmod(-z, DD) <= DE * DD)) { // south edge
        *(*pp)++ = cl0[6][0]; *(*pp)++ = cl0[6][1]; *(*pp)++ = cl0[6][2];
        return;
    }
     k = (int)((atan2(sn2, sn1) / pi4) * 2.999 + 12.0); // hard
    *(*pp)++ = cl1[k][0]; *(*pp)++ = cl1[k][1]; *(*pp)++ = cl1[k][2];

}

void drawcn1(unsigned char** pp, int i, int j) {
    double cn1, cn2, cn0;
    double z;
    int k;

    cn0 = tea[i] * tfa[j] + teb[i] * tfb[j];
    if (cn0 < anp) { // pole
        *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
        return;
    }
    cn1 = elkc / elk * te2[i] * te0[i] * tf2[j] * tf0[j];
    cn2 = -elkc / elk * te1[i] * te3[i] * tf1[j] * tf3[j];
    if ((cn1 * cn1 + cn2 * cn2) /cn0 / cn0 < asp) { // zero point
        *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
        return;
    }
    if (cn1 * (te2[i + 1] * te0[i + 1] * tf2[j] * tf0[j]) <= 0) { // pure imaginary
        *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
        return;
    }
    if (cn1 * (te2[i] * te0[i] * tf2[j + 1] * tf0[j + 1]) <= 0) { // pure imaginary
        *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
        return;
    }
    if (cn2 * (-te1[i + 1] * te3[i + 1] * tf1[j] * tf3[j]) <= 0) { // pure real
        *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
        return;
    }
    if (cn2 * (-te1[i] * te3[i] * tf1[j + 1] * tf3[j + 1]) <= 0) { // pure real
        *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
        return;
    }
    z = atan(sqrt(cn1 * cn1 + cn2 * cn2) / cn0) - pi4; // -pi/4 <= z < pi/4
    if ((z > 0.0) && (fmod(z, DD) <= DE * DD)) { // north edge
        *(*pp)++ = cl0[5][0]; *(*pp)++ = cl0[5][1]; *(*pp)++ = cl0[5][2];
        return;
    }
    if ((z <= 0.0) && (fmod(-z, DD) <= DE * DD)) { // south edge
        *(*pp)++ = cl0[6][0]; *(*pp)++ = cl0[6][1]; *(*pp)++ = cl0[6][2];
        return;
    }
    k = (int)((atan2(cn2, cn1) / pi4) * 2.999 + 12.0); // hard
    *(*pp)++ = cl1[k][0]; *(*pp)++ = cl1[k][1]; *(*pp)++ = cl1[k][2];
}

void drawdn1(unsigned char** pp, int i, int j) {
    double dn1, dn2, dn0;
    double z;
    int k;

    dn0 = tea[i] * tfa[j] + teb[i] * tfb[j];
    if (dn0 < anp) { // pole
        *(*pp)++ = cl0[1][0]; *(*pp)++ = cl0[1][1]; *(*pp)++ = cl0[1][2];
        return;
    }
    dn1 = elkc * te3[i] * te0[i] * tf2[j] * tf3[j];
    dn2 = -elkc * te1[i] * te2[i] * tf1[j] * tf0[j];
    if ((dn1 * dn1 + dn2 * dn2) / dn0 / dn0 < asp) { // zero point
        *(*pp)++ = cl0[2][0]; *(*pp)++ = cl0[2][1]; *(*pp)++ = cl0[2][2];
        return;
    }
    if (dn1 * (te3[i + 1] * te0[i + 1] * tf2[j] * tf3[j]) <= 0) { // pure imaginary
        *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
        return;
    }
    if (dn1 * (te3[i] * te0[i] * tf2[j + 1] * tf3[j + 1]) <= 0) { // pure imaginary
        *(*pp)++ = cl0[3][0]; *(*pp)++ = cl0[3][1]; *(*pp)++ = cl0[3][2];
        return;
    }
    if (dn2 * (-te1[i + 1] * te2[i + 1] * tf1[j] * tf0[j]) <= 0) { // pure real
        *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
        return;
    }
    if (dn2 * (-te1[i] * te2[i] * tf1[j + 1] * tf0[j + 1]) <= 0) { // pure real
        *(*pp)++ = cl0[4][0]; *(*pp)++ = cl0[4][1]; *(*pp)++ = cl0[4][2];
        return;
    }
    z = atan(sqrt(dn1 * dn1 + dn2 * dn2) / dn0) - pi4; // -pi/4 <= z < pi/4
    if ((z > 0.0) && (fmod(z, DD) <= DE * DD)) { // north edge
        *(*pp)++ = cl0[5][0]; *(*pp)++ = cl0[5][1]; *(*pp)++ = cl0[5][2];
        return;
    }
    if ((z <= 0.0) && (fmod(-z, DD) <= DE * DD)) { // south edge
        *(*pp)++ = cl0[6][0]; *(*pp)++ = cl0[6][1]; *(*pp)++ = cl0[6][2];
        return;
    }
    k = (int)((atan2(dn2, dn1) / pi4) * 2.999 + 12.0); // hard
    *(*pp)++ = cl1[k][0]; *(*pp)++ = cl1[k][1]; *(*pp)++ = cl1[k][2];

}

void draw1(void) {
    int i, j;
    unsigned char *p2, *p3, *p4;

    p2 = b2; p3 = b3; p4 = b4;
    for (j = 0; j < HT; j++) {
        for (i = 0; i < WT; i++) {
            drawsn1(&p2, i, j);
            drawcn1(&p3, i, j);
            drawdn1(&p4, i, j);
        }
    }
}

int main()
{
    FILE* ftp;
    errno_t err;
    unsigned int ui;

    // writing hedder
    b1[0x00] = 'B'; b1[0x01] = 'M';
    b1[0x02] = FZ & 0xff; b1[0x03] = FZ >> 8 & 0xff; b1[0x04] = FZ >> 16 & 0xff; b1[0x05] = FZ >> 24;
    b1[0x06] = 0; b1[0x07] = 0;  b1[0x08] = 0; b1[0x09] = 0;
    b1[0x0a] = 0x36; b1[0x0b] = 0; b1[0x0c] = 0; b1[0x0d] = 0;
    b1[0x0e] = 0x28; b1[0x0f] = 0; b1[0x10] = 0; b1[0x11] = 0;
    b1[0x12] = WT & 0xff; b1[0x13] = WT >> 8 & 0xff; b1[0x14] = WT >> 16 & 0xff; b1[0x15] = WT >> 24;
    b1[0x16] = HT & 0xff; b1[0x17] = HT >> 8 & 0xff; b1[0x18] = HT >> 16 & 0xff; b1[0x19] = HT >> 24;
    b1[0x1a] = 0x1; b1[0x1b] = 0;  b1[0x1c] = 0x18; b1[0x1d] = 0;
    b1[0x1e] = 0; b1[0x1f] = 0; b1[0x20] = 0; b1[0x21] = 0;
    b1[0x22] = BZ & 0xff; b1[0x23] = BZ >> 8 & 0xff; b1[0x24] = BZ >> 16 & 0xff; b1[0x25] = BZ >> 24;
    b1[0x26] = 0; b1[0x27] = 0; b1[0x28] = 0; b1[0x29] = 0;
    b1[0x2a] = 0; b1[0x2b] = 0; b1[0x2c] = 0; b1[0x2d] = 0;
    b1[0x2e] = 0; b1[0x2f] = 0; b1[0x30] = 0; b1[0x31] = 0;
    b1[0x32] = 0; b1[0x33] = 0; b1[0x34] = 0; b1[0x35] = 0;

    // initialize
    init1();

    // making bitmap data b0
    draw1();

    // output file
    err = fopen_s(&ftp, "sn.bmp", "wb");
    if (err != NULL) return err;
    for (ui = 0; ui < 0x36; ui++) {
        fputc(b1[ui], ftp);
    }
    for (ui = 0; ui < BZ; ui++) {
        fputc(b2[ui], ftp);
    }
    err = fclose(ftp);
    if (err != NULL) return err;
    err = fopen_s(&ftp, "cn.bmp", "wb");
    if (err != NULL) return err;
    for (ui = 0; ui < 0x36; ui++) {
        fputc(b1[ui], ftp);
    }
    for (ui = 0; ui < BZ; ui++) {
        fputc(b3[ui], ftp);
    }
    err = fclose(ftp);
    if (err != NULL) return err;
    err = fopen_s(&ftp, "dn.bmp", "wb");
    if (err != NULL) return err;
    for (ui = 0; ui < 0x36; ui++) {
        fputc(b1[ui], ftp);
    }
    for (ui = 0; ui < BZ; ui++) {
        fputc(b4[ui], ftp);
    }
    err = fclose(ftp);
    if (err != NULL) return err;
    err = fopen_s(&ftp, "el.txt", "w");
    if (err != NULL) return err;
    fprintf(ftp, "elliptic function constant\n");
    fprintf(ftp, "k = %10.5g, k^2 = %10.5g, k' = %10.5g, k'^2 = %10.5g\n", elk, elk * elk, elkc, elkc * elkc);
    fprintf(ftp, "K = %10.5g, K' = %10.5g\n", ellk, ellkc);
    fprintf(ftp, "q  = %10.5g, te2(0) = %10.5g, te3(0) = %10.5g, te0(0) = %10.5g\n", elq, te2a, te3a, te0a);
    fprintf(ftp, "q' = %10.5g, tf2(0) = %10.5g, tf3(0) = %10.5g, tf0(0) = %10.5g\n", elqc, tf2a, tf3a, tf0a);
    fprintf(ftp, "horizontal: WT = %6i, CT = %10.5g, D = %10.5g\n", WT, UC, UD);
    fprintf(ftp, "vertical:   HT = %6i, CT = %10.5g, D = %10.5g\n", WT, VC, VD);
    fprintf(ftp, "n pole (infinity) area = %10.5g, s pole (zero point) area = %10.5g\n", anp, asp);
    err = fclose(ftp);
    if (err != NULL) return err;
    printf("complete\n");
    return err;
}

　実行させるためにはコンパイル作業が必要です。Microsoft Visual Studio Community 2022 (64 ビット) を使いました。

　なお、アイデア自体には著作権はありませんので、改変後のプログラムなどはご自由に発表などなさって下さい。適当に数値をアレンジなどして出てきた図の著作権者は私では無く、あなたです。いずれも私への報告や著作権表示の必要はありません。ただし、本ブログに掲載されたプログラムや図をそのままの形で自分のものであると主張するのはお控え下さい。

4136. 楕円関数、複素平面編

　ヤコビの楕円関数の複素平面表示、はとても印象的な図が出てくるので、ネットに散見されると思います。これもその一つ。プログラムが一応の完成を見たと判断したので公開します。数学愛好家の趣味の範囲でご活用下さい。

　まず、図から。sn(u, k) (k = √(1/2))です。
(sn)

　これが最も典型的なヤコビの楕円関数の複素平面表示の結果と思います。

　数学的意味を知りたい方は、「森口繁一他。岩波 数学公式集 III。岩波書店(1960/1987)」の40ページの「sn(u, k)」の図をご覧下さい。本図では平行四辺形では無く2:1の長方形になっていて、中央が(0)、少し右の白線の十字路が(K)、その右の青系の灰色の楕円が集中している所が(2K)、画面左側の同様の絵の所がこれも(2K)で、(3K)を経て中央に戻ります。つまり、実数軸方向には4Kの周期関数であることが分かります。
　今度は中央(0)から上に行くと、赤系の楕円が集中している所が(iK') (iは虚数単位)で、その上に再び中央と全く同じに青系の楕円が集中している所が(2iK')で、虚数軸方向には2K'の周期関数であることが分かります。つまり、関数sn(u, k)は定義域が複素平面で4Kと2K'の二重周期を持つ関数です。

　ちなみに、今回上げる3図の母数(modulus) kの値はいずれも√(1/2)です。k値を変えると等高線などが違う様相を見せますが、それは近々ご覧いただけると思います。

　次が、今回最も印象的だったcn(u, k) (k = √(1/2))の図。
(cn)

　cn(u, k)の周期は、実数軸方向(中央から左右へ)は4Kで、虚数方向は右上45°方向に見られ、(2K+2iK')です。中央と同じ図柄が見えるでしょうか。こちらは周期平行四辺形(群論で言う基本領域)が横が4K、縦が2K'の高さの平行四辺形になっています。

　この図cn(u, √(1/2))で、図の読み方を解説します。中央の白線がクロスしている所がガウス平面(複素平面/複素数平面)の中心、(0+0i)です。水平線が実軸(純実数)です。k = √(1/2)の時のK値は1.8541なので、右の青系の楕円が集中してる(K)の座標は(1.8541+0i)です。中心から上下に虚軸が走っていて、この白線の部分が定義域の純虚数です。中央上の赤系の楕円が集中している所が(iK')で、K'値は1.8541なので座標は(0+1.8541i)です。

　K=K'なのはk=√(1/2) (0.70711)の場合だけです。他のk(0≦k≦1。楕円関数の母数。楕円の離心率に相当)ではK≠K'となります。本プログラムでは周期を固定しているので、他のkを与えた場合も図の上ではKとK'が等長に表示(図の例は次回以降に)されます。

　まず、中央付近の図が気になると思いますので先に解説します。白線は関数値が純実数(虚数=0)の場所です。斜め45°方向に黒線がクロスしている場所が見えると思います。黒線は純虚数(実数=0)の場所です。ヤコビの楕円関数ではどちらも水平線と垂直線になってしまうようです。
　さっきから赤系・青系の灰色の楕円と言っているのは関数値の絶対値の等高線です。中央では隣接していて、境界部は絶対値1の場所です。中央のcn((0+0i), k)の関数値は(1+0i)で、そのことを示します。中央右の青系の灰色の楕円が集中している(K)の場所は(K+0i)で、関数値は0 (0+0i)です。少しだけ明るい青系の丸が見られて、ここが零点、つまり関数値(0+0i)付近ということを示します。この青丸は等高線とは別に打っています。
　中央から上、(0+iK')の場所の赤系灰色の丸は極、つまり無限大というか分母が0になる点の付近です。

　ですから、青系の灰色の等高線図は絶対値が0～1の、リーマン球面で言えば南半球の部分です。赤系は1～∞の北半球です。色はプログラムで適当に変更できます。この赤系・青系は学習用棒磁石の北極がなぜか赤色が多いようなので、そうしました。無限大(極)が赤色点で、0 (零点)が青色点です。

　背景の虹色は関数値の複素数の偏角を示します。マゼンタと赤の中間のピンクの部分が0度(正の実数)付近、黄が90度(正の虚数)付近、シアンと緑の中間の草色の部分が180度(負の実数)付近、青が270度(負の虚数)付近です。

　最後がdn(u, √(1/2))のグラフです。
(du)

　何となくdn(u, k)とsn(u, k)のグラフが同様の感じがするのはその通りだと思います。3図とも極の位置が同じなのは注目点です。

　最後に、3図と同時にプログラムが出力する定数のテキストファイルを示します。

elliptic function constant
k =    0.70711, k^2 =        0.5, k' =    0.70711, k'2 =        0.5
K =     1.8541, K' =     1.8541
q  =   0.043214, te2(0) =    0.91358, te3(0) =     1.0864, te0(0) =    0.91358
q' =   0.043214, tf2(0) =    0.91358, tf3(0) =     1.0864, tf0(0) =    0.91358
horizontal: WT =   1024, CT =    0.00125, D =     0.0025
vertical:   HT =   1024, CT =    0.00125, D =     0.0025
n pole (infinity) area =      0.001, s pole (zero point) area =      0.001

　k^2はk値の自乗です。k'は補母数。qは計算途中で出てくる楕円テータ関数の入力変数で、本来はk値から算出されるものです。te…とtf…は計算に必要な楕円テータ関数の定数です。
　WTとHTは図の横と縦のドット数。CTは中央の値で、Dは隣接するドットとの差分値です。多少の理由があって、半ドット分中央をずらせています。0に指定することはできます。最下行は極と零点のドットの大きさです(単位は独特で、グラフを見た感じで調整する値です)。

4135. 複素数編へ、続き^3

　本日は土曜日で私は休日です。ちょっと奮起して、ヤコビの楕円関数の複素平面表示のプログラムを作りました。それなりの絵が出てきていますが、もう少し確認して踏ん切りが付いたら公開します。
　適当に作ってしまったので、整理が付いていないプログラムの感じがするのです。しかし公開しない理由は無いので、ある時点で妥協すると思います。

　級数が必要なのはわずかに8ライン分なのですが、それでも現代の計算機の速度には感銘を受けます。良い時代になりました。