光学薄膜設計ソフトの設計 その26 - 複素演算の独自実装のとりやめ

こないだベクタユニットの効率の評価をかねて、低レベルの複素数演算を書いた。けど、それは撤回、もうやめ、という話。

すっかり忘れていたんだけど、最新のCの規格であるC99には複素数型が含まれている。

僕が書き始めた頃は「K&R」が標準だった。かなり使い込んでしまってから「安西さんが標準化するらしい」という噂が聞こえてきた。僕はK&RからANSI Cへの乗り換えはそれほど問題なかったけど、厳しめの型チェック（今から思えばたいしたことない）に悩んだ人もいた。

K&RでもANSI Cでも、複素数型は自分で実装する必要があった。

僕が普段使っている数値計算関連のライブラリ（主にfftwとgsl）はそれぞれ自前で複素数を定義している（fftwはC99の複素数を使うようにも書ける）。fftwとgslの複素数型はどちらも単純な実数部と虚数部の配列なので、これまでキャストするだけでどちらでも使える自前の定義をずっと使っていた。こないだ書いた複素演算も実は結局、その昔書いたソースを丸コピーしてきただけだった。

C99の複素数型は、たとえば虚数単位のIや、実数型にキャストすると実部が返されるとかという、自前の実装では実現できない便利な機能がある。もちろん四則演算の演算子もそのまま使えるし、実数型と混ぜてもちゃんと型変換される。実質的にC++のcomplex標準ライブラリと同じことができるようになっている。

当然効率も十分考えてライブラリとして実装されているはずなので、使った方がいいだろう。ということで、これまでC99は意識して使って来なかったけど、今回から対応することにしよう（とは言ってもこれまでも//で始まるコメント行とか、使う直前での変数宣言とか、結局使っていたけどね）。

まあ、Objective-Cはいい（衝突する機能がないし、それにマイナーだし）として、先に複素数を標準ライブラリとして実装したC++は、C99とは複素数型以外にも衝突している機能があるらしくてすでに互換性はないらしい。メジャーなC++がC99と互換性がないということは規格はどうなるんだろう。C++は一方でどんどん拡張を続けて普段使っていない僕なんかは言語仕様が大きすぎてまったく把握できない。C99は規格としては存在してるけど誰も使わないラテン語みたいなものになるんだろうか。

特性マトリクスのかけ算の例

ということで、こないだ書いた複素数のルーチンは全部C99の複素数型で置き換えよう。今回ベクタユニットを使った複素数計算をスカラユニットを使うようにコンパイルした場合と効率を比較しようと思ってたけど、数値計算の大部分は複素演算が律速するので意味がなくなった。

gccのautovectorizeの件もあるのでその比較はもういいだろう。

試しにC99の複素数を使ってマトリクスのかけ算を書いてみると

#include <complex.h>

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

となって、ずっと読みやすい。あたりまえだけど。