光学薄膜設計ソフトの設計 その24 - 複素数演算について

具体的な薄膜を表すクラスの実装を始める前に、光学の数値計算では必ず必要になる複素数の扱いについて。

専用の複素数演算ルーチン

Objective-Cでの複素数演算

複素数の演算機能は、Cと同じでObjective-Cには備わっていない。FORTRANや、C++の演算子Overloadを使えば

    complex  z, a, b;
    z = (a - b) / (a + b);

ですむのが、Objective-Cではこれを実数の演算に展開する必要があって、例えば

    complex  z, a, b;
    z = complexDiv(complexSubt(a, b), complexAdd(a, b));

とかのように構造体を返す関数を組み合わせて書かなければいけない。これはデバグを考えるとめちゃ不利である。

しかし、問題はそれぞれ固有の構造を持っていて、プログラミングの過程では計算だけでなく問題の構造を含めて表現できないといけない。今回の場合では、多層膜という材料と厚さが違うだけの繰り返しの表現、最適化のときの変数の管理、マルチコアに対応するための依存関係の分離、などなど、計算そのものから見れば雑用なんだけど、ある程度汎用性があって、デバグしやすいようにするためにはかなりの作業をしなければならない。

そのためにはやはり使い慣れた言語と環境でするのが一番。ということで僕の場合、一番ふさわしいのはMathematicaということになるんだけど、これはコンパイラでハードコードした計算に較べると10進二桁から三桁遅い。計算資源を膨大に持っているお金持ちでない僕は次善の策としてObjective-Cを使う、ということになる。

今回の実装

最終的にはベクタユニットを使いたいので複素数の演算と、2行2列の行列演算は専用のルーチンを書いておく。この演算ルーチンを置き換えることでベクタユニットを使う実装とスカラユニットの実装とを差し替えて効率の評価ができるようにしておく。もちろん、ここ以外でこのレベルの計算は勝手に定義しないようにする必要がある。

めんどくさいけど、複素数演算は普通のCの関数として定義する。

typedef struct {
    double  re;
    double  im;
}   complex;

complex	makeComplex(double r, double i);
complex	complexAdd(complex a, complex b);
complex	complexMult(complex a, complex b);
complex	complexSubt(complex a, complex b);
complex	complexReciprocal(complex a);
complex	complexDiv(complex a, complex b);
complex	complexConj(complex a);
complex	complexMinus(complex a);
complex	complexSin(complex a);
complex	complexCos(complex a);
complex	complexExp(complex a);
complex complexSqrt(complex a);

double  complexSquareAbs(complex a);
double  complexArg(complex a);

BOOL    isPureReal(complex a);
BOOL    isPureImaginary(complex a);
BOOL    isComplexZero(complex a);

というような、どこにでもあるユーティリティ関数群。何をするかは関数の名前から明らかだと思う。これらの関数は戻り値として結果を返すのでスタックを消費することになる。まあ、この程度までなら許せるでしょ。Cの戻り値はレジスタを使うんだっけ？

もちろん、使うやつだけ中身を書くんだけど、ちなみに最初のふたつは

complex	makeComplex(double r, double i)
{
    complex a;
    a.re = r;
    a.im = i;
    return a;
}

complex	complexAdd(complex a, complex b)
{
    complex c;
    c.re = a.re + b.re;
    c.im = a.im + b.im;
    return c;
}

みたいなもの。こないだのことがあるので、このままだときっと簡単な関数はベクタユニットを使うようにコンパイルされてしまうのだろう。あとで試してみないといけないけど。

行列とベクトルは、特性マトリクスとアドミタンスしか使わないので配列としてではなく構造体にしてサイズを固定することにする。

typedef struct {
    complex e00;
    complex e01;
    complex e10;
    complex e11;
}   characteristicMatrix;

typedef struct {
    complex a;
    complex b;
}   admittanceVector;

//  mret = m1 * m2
void    chMatrixMult(characteristicMatrix *mret,
                     characteristicMatrix *m1,
                     characteristicMatrix *m2);
//  mret = m[num - 1] * ... m[1] * m[0]
void    chMatrixSeriesMult(characteristicMatrix *mret,
                           characteristicMatrix **m,
                           unsigned int num);
//  vret = m1 * v
void    chMatrixMultWithVector(admittanceVector *vret,
                               characteristicMatrix *m1,
                               admittanceVector *v);
//  vret = m[num - 1] * ... m[1] * m[0] * v
void    chMatrixSeriesMultWithVector(admittanceVector *vret,
                                     characteristicMatrix **m,
                                     unsigned int num,
                                     admittanceVector *v);

ぐらいを定義しておく。ひとつめはふたつのマトリクスの積、ふたつめはnum個の積、みっつめはベクトルとマトリクスの積、よっつめはベクトルにnum個のマトリクスをかけた結果を返す。これ以外は必要ないだろう。どれも結果は、関数の戻り値としてではなく呼び出し側がメモリを確保するとする。これは今回、計算結果をキャッシュすることにしたから、行列やベクトルを使うどのオブジェクトも自分で持っているので問題ない。

ちなみに最初のふたつの実装は

void    chMatrixMult(characteristicMatrix *mret,
                     characteristicMatrix *m1,
                     characteristicMatrix *m2)
{
    mret->e00 = complexAdd(complexMult(m1->e00, m2->e00),
                           complexMult(m1->e01, m2->e10));
    mret->e01 = complexAdd(complexMult(m1->e00, m2->e10),
                           complexMult(m1->e01, m2->e11));
    mret->e10 = complexAdd(complexMult(m1->e10, m2->e00),
                           complexMult(m1->e11, m2->e10));
    mret->e11 = complexAdd(complexMult(m1->e10, m2->e10),
                           complexMult(m1->e11, m2->e11));
}

void    chMatrixSeriesMult(characteristicMatrix *mret,
                           characteristicMatrix **m,
                           unsigned int num)
{
    characteristicMatrix    tmp;
    int i;
    
    *mret = *(m[0]);
    if (num == 1)
        return;
    for (i = 1 ; i < num ; i ++) {
        chMatrixMult(&tmp, *(m + i), mret);
        *mret = tmp;
    }
}

というようなもの。ああ、うるさい。でもまあ、めんどうだけど簡単。自前のルーチンで呼ばれるところもわかるのでnum引数の範囲チェックははしょる。スタックを食う複素数演算は関数を引数の中に入れ子にできるので便利。

実装はまだ全部を書き下ろさないけどスカラユニットを使うものとベクタユニットを使うものの二通りを用意しておくことにする。complexをはじめからvDoubleで宣言する、とかも試してみたい。それは足し算引き算はいいけどかけ算はperformance hazardになるという話もあるけど。