MININEC v2 (QB45) — クリーンな版を公開する
ここ数十年で私が手を加えてきたアンテナ計算コードは——MiniNEC x64となったQB64版の書き直しも、今回jMiniNECとして書き上げたJava移植版も——すべて、2003年から2007年頃に私が書いた1本のBASICコードに遡ります。これは1990年代初頭のさらに古いQB45プログラムを、手続き型に書き直したものです。その古い方の非手続き型バージョンは、jMiniNECのエンジンでは今「v1」と呼ばれており、すでに引退しています。こちら——MININEC v2と呼んでいるもの——が、私が実際にその後もずっと開発を続けてきたバージョンであり、私以外の誰かがゼロからMININECを移植したいと思ったときに重要になるのはこちらです。
以下のファイルで変更した点はただ1つ、GOTO文と行番号によるジャンプをなくし、IF/END IF、DO/LOOP、EXIT FORに置き換えたことだけです。これは書き直しではなく、機械的な整理作業です。すべての分岐とループは、元のコードとまったく同じ条件で、まったく同じ回数だけ実行されます。計算内容は一切変わっていません。これが重要なのは、ここにあるのが今でも忠実に、NOSC TD 516のアルゴリズムそのものであって、私が「こう書くべきだった」と考える改変版ではないということだからです。
なぜMiniNEC x64ではなくこちらなのか
MiniNEC x64は私が実際に使っているバージョンで、QB64への移植の中でかなり肥大化しています。CSV出力、追加の診断機能、使い勝手を良くするための入力処理など——これらはどれも、MININEC 2や3の本来の仕様には含まれていないものです。だからこそ、MININECを土台にした私個人のツールではなくMININECそのものを理解したい、あるいは再実装したいという人に渡すには不向きなのです。余計な機能が1つ増えるたびに、その下にある本来のアルゴリズムを見るために、将来の移植者が頭の中で差し引かなければならないものが1つ増えることになります。
このファイルには差し引くべきものが何もありません。8個のSUBと、カーネル計算用の2個のFUNCTIONで構成され、NOSC TD 516に書かれていることをそのまま行い、それ以上のことはしません。MININECをゼロからPython、Rust、Cへ移植しようとしている方、あるいは特定の誰かのバージョンではなくMININECそのものが何であるかを理解したいだけの方にとって、これは私の他の2つのバージョンよりも近道になります。
DECLARE SUB Intro ()
DECLARE SUB EnterData ()
DECLARE SUB WireConnect ()
DECLARE SUB MainEngine ()
DECLARE SUB PrintCurrentData ()
DECLARE SUB Redo ()
DECLARE SUB FUNCTION1 ()
DECLARE SUB FUNCTION2 ()これがプログラム全体の姿を、冒頭でそのまま示したものです。EnterDataとWireConnectが形状を構築し、MainEngineがインピーダンス行列を埋めて解き、PrintCurrentDataが電流・インピーダンス・遠方界パターンを出力し、FUNCTION1とFUNCTION2がMainEngineの行列構築の奥深くから呼び出されるカーネル積分を担います。GUIにも設定ファイルにも出力形式にも何も隠されておらず、8個の手続きと数千行の算術演算だけで、実際の電磁気計算コードに必要な分量ちょうど、それ以上でもそれ以下でもありません。
Wilton–Pocklington分割と、1980年にそれが必要だった理由
ポックリントンは1897年に、薄い導体ワイヤ上の電流に関する厳密な積分方程式を導きました。この方程式自体を書き下すのは難しくありません。難しいのは、それを安価に評価することです。というのも、カーネル——ワイヤの円周方向に平均化した自由空間グリーン関数——は、モーメント法のコードが最も必要とする場所、すなわち観測点がソースセグメントそのものの上にある自己項において、まさに対数特異点を持つからです。これを、特異点近傍の全域にわたる完全な数値処理で正直に評価しようとすれば、注意深い適応求積、あるいは本格的な数値グリーン関数処理が必要になります——まさに、1970年代のローレンス・リバモア研究所で使えたミニコンピュータやメインフレーム上で動くNECがやっていることです。
MININECの著者たちには、そのようなハードウェアを対象にする余裕はありませんでした。1980年のNOSCの取り組み全体の目的は——原論文の言葉を借りれば「アンテナモデリングコードの縮小版につながりうる技術」を調べることであり、実用的なモーメント法ソルバーを、BASICで、現場で、数十キロバイトのRAMしかなく浮動小数点演算コプロセッサも持たないマイコン上で動かすことでした。N×N行列の要素1つひとつについて毎回計算し直す、超特異なカーネルの本格的な適応求積処理など、そもそも選択肢に入りようがありませんでした。
それを可能にした技法——私が「Wilton–Pocklington分割」と呼ぶようになったもので、元のMININEC論文が、当時ミシシッピ大学にいたDonald R. Wilton博士を、Julian、Logan、Rockwayと並ぶ共著者として直接クレジットしているためこう呼んでいます——が、以下のFUNCTION1の裏にあるトリックです。特異なカーネルを直接積分する代わりに、特異点を解析的に(閉形式で積分可能な素の対数として)差し引き、近傍界の残り部分をよく知られた楕円積分の有理式+対数近似(ゼロから計算するのではなく、一度だけ参照する10個の固定定数)で置き換え、なめらかな残りの部分だけを数値求積に渡します。特異で閉形式のもの1つと、なめらかで安価なもの2つ、合わせて3つのパーツに分けることで、特異点付近でサンプリングするのが全く安全でない厄介な1つのパーツを避けています。
これこそが、このプログラムが1980年代のハードウェア上で存在しえた理由のすべてです。NECの完全/厳密な細線カーネルはより精度が高く、特に太いワイヤに対して優れていますが、その分マシン側にも相応の性能が必要でした。MININECの簡略化されたカーネルは、その精度をわずかに犠牲にする(楕円積分近似は良好ですが厳密ではありません)代わりに、各行列要素の評価を、適応数値積分ではなく数十回の浮動小数点演算で済ませます。FPUのないマシンでは、これは些細な最適化ではなく、一晩で終わるプログラムと永遠に終わらないプログラムの違いそのものです。
この分割が、このファイルの中に実際どう残っているかがこちらです——FUNCTION2のカーネル積分から、行列要素ごとに1回呼び出されるFUNCTION1です:
SUB FUNCTION1
Rem Elliptical integral?
X3 = T * D(P4)
Y3 = T * F(P4)
Z3 = F(P4) * T * G(P4)
Den = D0 - 2 * F3 * (X3 * X2 + Y3 * Y2 + Z3 * Z2) + X3 * X3 + Y3 * Y3 + Z3 * Z3
Den = SQR(Den)
IF I6 <> H0 THEN
D3 = Den * Den - A2
B = D3 / (D3 + 4 * A2)
V1 = W(1): Rem A0
V2 = W(6): Rem B0
B1 = 1
B2 = -36
FOR L2 = 2 TO 5
B2 = B2 / (L2 - 1)
IF B < 10 ^ (B2) THEN EXIT FOR
B1 = B1 * B
V1 = V1 + W(L2) * B1
V2 = V2 + W(L2 + 5) * B1
NEXT L2
V = (V1 - V2 * LOG(B)) * SQR(1 - B)
D3 = SQR(D3) / 8 / A(P4)
T3 = T3 + (LOG(D3) + V) / P / A(P4) - 1 / Den
END IF
B1 = Den * Wavelength
T3 = T3 + COS(B1) / Den
T4 = T4 - SIN(B1) / Den
END SUBVは楕円積分近似で、ファイル冒頭近くのDATAとして読み込まれる10個のW()定数から組み立てられる、「素性の良い近傍界の残り」項です。その隣にあるLOG(D3)が閉形式の特異な部分、そして最後の2行(DenのCOS/SIN)が3つ目のなめらかな遠方界の部分です。この分割の3つのパーツすべてが、わずか11行の中に収まっており、厳密なカーネル評価であれば本格的な数値積分を必要とするところを、ほんの数回の浮動小数点演算でこなしています。なお、これは転記上の癖も含め、元のファイルにあったとおりそのまま公開しています——最後のbx行のDATA定数は参照テーブルからごくわずかにずれています(これは私の誤植ではなく桁の入れ替わりです)。気になる方のために、その修正内容と(無視できるほど小さな)影響についてはjMiniNECのページに詳しく書いてあります。
次にこれを移植する人へ、プログラマからのメモ
このファイルのGOTO整理作業を行っていたとき、制御フローが妙に見える箇所が2つ見つかりました。リファクタリングの途中でまだ完全に理解していないものを黙って「修正」してしまうのではなく、そこで手を止め、自分自身——あるいは後から来る誰か——のためにメモを残すことにしました。両方ともそのままの状態で残し、ファイルの末尾に次のように書いておきました:
1) H0 is used both as a plain "0" constant early on and is later set to 1
inside MainEngine, which changes the lower bound of several "FOR I = H0
TO N" loops in the print/pattern routines. That's the original code's
doing, not something introduced by this GOTO cleanup.
2) In MainEngine's pivot-search loop, the pivot-scan block is still guarded
by "IF P = 0 THEN ..." (originally "IF P THEN GOTO ..."). Because P
holds Pi (set near the top of the program) it is never zero, so that
block never runs - exactly as in the original GOTO version. Renaming P
there was presumably meant to be a separate flag; left exactly as-is.そのメモは、実際にどちらの謎も突き止めるまで、このファイルの中で何年もそのままになっていました。どちらも単なる書き方の癖ではなく、本物のバグだったことが判明し、両方ともjMiniNECのページで詳しく解説しています。奇妙な点その1は、H0(Java移植版ではloopBaseに改名)をMainEngineを呼び出すたびに明示的にゼロへリセットしなければならない理由です。放置すると、前回の呼び出しの古い値が行列のゼロクリアに漏れ込み、原点に完全に中心を置いていない形状では結果が壊れてしまいます。奇妙な点その2はさらに厄介で、このコードが明らかに備えているはずだった「完全ピボット選択」——各消去ステップの前に、現在の列だけでなく残りの部分行列全体から最大の要素を探す処理——が、このプログラムのどのBASICバージョンにおいても、一度も実際には実行されていなかったことを意味します。jMiniNECでは、たまたま探索を正しく実装していた、まったく別系統の私のJavaエンジンから移植する形でこれを復元しています。
これは指摘しておくに値することだと思います。上記のコメントは、後知恵を先見の明であるかのように装ったものではありません。「何かがおかしいが、機械的な整理作業の途中でそれを止めて直しはしない」という、まさにその場で書いた本物のメモであり、どちらの問題も実際に理解されるより何年も前に書かれたものです。次にこのファイルを移植する人は、両方の奇妙な点を、記録されたとおりそのまま引き継ぐことになります——これこそ、リファレンスソースが引き継ぎを行うべき正しいやり方だと思います。名前を付け、日付を記し、誰かが実際に直す準備ができるまで手を付けずに残しておくのです。
あえて載せていないこと
このページにはスクリーンショットのセクションもなければ、それに合わせて新たにDOSBoxで実行し直すこともしていません。エミュレータを引っ張り出してこのファイルをそのまま再実行することもできますが、MiniNEC x64のスクリーンショットや、公開されている計算例および独立したPython再実装に対するjMiniNECの検証済みの突き合わせ結果によって、このアルゴリズムそのものの正しさ(上記2つの記録済みの奇妙な点を除く)は、このサイトの他の場所ですでに確立されています。このページはMININEC v2が動くことを証明するためのものではなく、その上に構築された下流のツールとは切り離された形で、ソースコードそのものが恒久的で引用可能な置き場所を持つためのものです。
まとめ
というわけで、このサイトには今、3つのMININECコードベースが存在し、それぞれ別々の理由でその居場所を得ています。MiniNEC x64は実用ツール、jMiniNECは現代的でドキュメント化された、バグ修正済みのエンジン。そしてこれ——MININEC v2、地味なQB45で、GOTOを取り除いた以外は何も手を加えていないもの——はリファレンスです。NOSC TD 516のアルゴリズムが実際どのようなコードなのか、機能を何も追加せず、バグも黙って修正されていない状態で知りたいなら、まずはこのファイルから始めるのが良いでしょう。