Entry

プログラミングメモ - 複素数の計算方法とかちょっと

2008年10月14日

昨日 Amazon から届いた本が,ちとたくさんあるので,順番にやっつけてます。とりあえず,『数値計算の常識』を,プログラミングしながらゴニョゴニョと。1985年初版で,あたしが買ったのが2007年9月15日刷の第26刷ですから,その手の手合いには超おなじみの本なんだと思います。

数値計算の常識
数値計算の常識
posted with amazlet at 08.10.14
伊理 正夫 藤野 和建
共立出版
売り上げランキング: 17864
おすすめ度の平均: 5.0
5 文字どおり「常識」にしたい内容
5 困る前に

数値計算の本は,『C言語による数値計算入門』を使ってたんですけれど,こちらはアルゴリズムや方程式の解き方のように,具体的な数式をコンピュータに計算させる方法が中心になっていたのでした。コンピュータで常微分方程式を解きたいとか,関数近似や補間をしたいとかいった具合に,具体的にやることが決まっている向きには便利なんですけれど,あたしのようにコマコマとつまみ食いしているような似非計算機使いにとっては,使う箇所と使わない箇所の落差が激しかったりします。

C言語による数値計算入門―解法・アルゴリズム・プログラム (UNIX & Information Science)
皆本 晃弥
サイエンス社
売り上げランキング: 28934

一方,『常識』の方は,もっとプリミティブで,桁落ちや速度を考慮した数式の変形方法や,リソースを目いっぱいに使うための基本的なテクニックが紹介されています。いわば,『入門』が各論で,『常識』が総論といったところ。

で,ここでは複素数の計算をしてみました。今,懲りずにフーリエ変換をゴニョゴニョやってるんですけど,フーリエ変換では複素数の演算が欠かせません。この点,実数の計算は,どうしても精度が気になるところ。sqrt(3) や pow(3) なんかをバリバリ挟むと,どうしても「こんなんでいいのか?」と不安になってきます。特に,絶対値の計算なんですけれど,どうも精度を取れていない気がする……。

そんなこんなで,本書を読んだところ,やっぱり公式通り(|z| = √(Re z)2 + (Im z)2)に計算するのはマズかったみたい。

とりあえず比べてみる。処理系は IA-32 Windows xp cygwin gcc 3.4.4 です。

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

double
compabs(double re, double im)
{
    double ret;
    double are = pow(re, 2);
    double aim = pow(im, 2);

    if (are == 0.0) {
        if (aim == 0.0) {
            ret = 0.0;
        } else {
            ret = sqrt(aim);
        }
    } else {
        if (aim == 0.0) {
            ret = sqrt(are);
        } else {
            ret = (are > aim)
                ? (sqrt(are) * sqrt(1 + pow(im / re, 2)))
                : (sqrt(aim) * sqrt(1 + pow(re / im, 2)));
        }
    }

    return ret;
}

double
normalabs(double re, double im)
{
    return sqrt(pow(re, 2) + pow(im, 2));
}

int
main(int argc, char *argv[])
{
    double re = 9999991234.01001;
    double im = 239876543.10005;

    printf("Re : %f\n", re);
    printf("Im : %f\n", im);
    printf("compabs    : %f\n", compabs(re, im));
    printf("normalabs  : %f\n", normalabs(re, im));

    return 0;
}

本書に書かれていた方法を実装したのが compabs(),公式通りの方法を実装したのが normalabs() です。compabs() のポイントは,(Re z)2 や (Im z)2 を直接計算しないで,一度割り算していること。これで指数部あふれを防げます(と書いてある)。もちろん,割り算するので,0割りしないように場合分けしています。

つことで,結果。

aian@IBM-39B80170C23 ~
$ ./a.exe
Re : 9999991234.010010
Im : 239876543.100050
compabs    : 10002867860.579117
normalabs  : 10002867860.579119

少ない桁数だと違いはなかったんですけれど,有効桁数ぎりぎりで使うとそれなりに誤差が出るようです。本書が書かれた時期は1985年なので,実数計算はおそらく倍精度(double 型)ではなく,単精度(float 型)を使うのが一般的だったんだと思います。float 型で計算したら,もっと早くに誤差が出るのかも。

それにしても,桁数ぎりぎりで演算しても,小数点第5位程度の精度は出ているみたい。つことは,あたしんとこのフーリエ変換で出てる不具合は,このせいじゃない感じ(明らかに変な値が出てる)。どこがバグってるんだぁー!

Trackback
Trackback URL:
Ads
About
Search This Site
Ads
Categories
Recent Entries
Log Archive
Syndicate This Site
Info.
クリエイティブ・コモンズ・ライセンス
Movable Type 3.36
Valid XHTML 1.1!
Valid CSS!
ブログタイムズ

© 2003-2012 AIAN