バタフライ演算

質問

現在、c言語でバタフライ演算を行っているのですが、要素が8個のときはうまくいくのですが、16や32個以上になるとうまくいきません。どこが間違っているのかわからないのですが、わかる方よろしくお願いいたします。

プログラムは以下です↓↓↓
http://kissho.xii.jp/1/src/1jyou81037.cpp

回答No.4

Gotthold

> 「log、pow、sin、cosは信じるな」
> こいつら、実は、けっこうな誤差を出してくれる。
> 使うなら「有効数字の最下位桁を四捨五入」とかした方が良い。
四捨五入なんてしたら誤差を大きくするだけじゃないの？
(真の値が十進数表現できりの良い値であると分かっているなら別だけど。)
最終結果を見せるまえに四捨五入するってことならいいんだけど、
演算途中でそういう丸めはしなくて良いと思う。

お礼

解決しました。
やはりlogとって結果を表示しているプログラムですし、
実数をそのまま添字にしているようなところはなかったので
誤差は関係ありませんでした。
バタフライ演算は間違っている点はなく、
ソートのプログラムにバグがありました。
ご回答ありがとうございました。

回答No.3

Tacosan

math.h の何を使えというんだろうかとあえて 1点だけ突っ込む＞#1.

回答No.2

asuncion

＞16や32個以上になるとうまくいきません。

そもそも、「うまくいきません」とは、どういう現象のことを指していますか？
コンパイルができないのですか？
実行時に落ちてしまうのですか？
落ちはしないが、想定と異なる答えが出てしまうのですか？
何かエラーメッセージが出るのですか？

お礼

コンパイルは通ります。要素が８個のとき、つまり
バタフライ演算の３段目までは想定通りの答えがでるのですが、
４段目から異なる答えが出てしまいます。
複素数による符号反転を十分に考慮できていないのか、
原因がわからなくて質問させていただきました。

回答No.1

chie65536（@chie65535）

＞どこが間違っているのかわからないのですが、わかる方よろしくお願いいたします。
間違ってるのは「浮動小数点数の誤差を考慮せず、数学的な数式・理論をそのままプログラムしただけ」ってのが間違い。

間違いを正したいなら
「実数演算は必要最低限に」
「定数を元に決まった演算をして毎回固定の値を求めるなら、最初から定数にしておく」
「実数から整数への変換は行わない。行うなら、桁落ちと丸め誤差を考慮する。場合によっては四捨五入する」
などの改善が必要。

・質問者さんにアドバイス

「log、pow、sin、cosは信じるな」
こいつら、実は、けっこうな誤差を出してくれる。

使うなら「有効数字の最下位桁を四捨五入」とかした方が良い。

「実数から整数への変換は信じるな」
例えば「0.1×1000」のように、論理的に言えば100.0になる数であっても「浮動小数点数の特性」により、この値は99.99999999999999999999999872とかになってしまう可能性がある。

これをそのまま「整数の変数に代入する」とか「配列の添え字に使う」とかで整数に変換したら、結果は100ではなく99になる。

そういう事が起きたら、もう、結果は言わなくても判ると思う。

「どうしても実数じゃなきゃダメな部分以外は、実数を使うな」
理由は言わなくても判ると思う。

・改善すべき点をいくつか

#define PI 3.1415926535
もうね、ハナから「論外」なんだわ。

IEEE 754 のdouble値の有効桁数は何桁か知ってる？約15桁だよ。

貴方が書いたPI、桁が何桁ある？たったの11桁だよ。

安全を考えたら、17桁か18桁は書かないと。

足りない6～7桁が「縞模様として結果画像に浮き出る」って判ってないよね（判ってたら質問なんかしないか）

#define PI 3.1415926535
を
#define PI 3.14159265358979323846
に変えるだけでも、かなり「縞模様が減る」と思うよ（それ以外にもダメダメな点が多過ぎるから、縞模様は完全には消えないけどね）

#define N 256
って定数で書いてるなら
int M=log((double)N)/log(2.0);//Nは2の何乗か変数
みたいな無駄な演算はするなよ。定数で書け定数で！

logと除算に誤差があれば、計算結果の実数を整数に代入した瞬間、結果は想定外な物になる。Mが7になっても責任は取れないぞ。

要素が16や32になって動かないのは、こういう「バカな実数演算」をしているから。

こういうのは、ハナから
#define N 256
#define M 8
と、両方とも定数にしちゃうか
#define M 8
int N = 1 << M;
と書いて「実数演算を排除し、かつ、絶対に誤差の出ない式」を使うこと。

mul=(int)pow(2.0,(double)n);
だから、どうしてこういうのを実数で計算しちゃうかな。

誤差が出て、powが返す「2の9乗」が511.9999999999999999999999999984とかになったら、それをintに変換した瞬間、結果は511だぞ。

2のn乗が欲しいなら
mul = 1 << n;
って書きなさいよ。

cos(2*PI*k/mul)
sin(2*PI*k/mul)
結局、PIは２倍した値しか使ってないよね。

だったら最初から
#define M_2_PI 0.636619772367581343076
とか、PI×2を定義しちゃえば？

つ～か、math.hに定義されてるの使えよ。

FFT[i][j]=pow(pow(Re_countx[i][j],2)+pow(Im_countx[i][j],2),0.5);
どうして、こうも「誤差が大きい関数」ばっか好き好んで使うかな。

0.5乗するならsqrt関数にしようよ。こっちの方がまだマシ。

それと「２乗するのにpowを使う」ってのは「論外中の論外、もう、大論外」だわ。呆れて開いた口が塞がらない。

普通に
Re_countx[i][j]*Re_countx[i][j]
って書いて２乗した方が、精度が１００倍くらい違うぞ。

最後の最後に一言。

「質問者さんに対し、pow関数の使用を禁止する」

この禁止事項を守るだけで、かなり精度が高くなると思うぞ(笑)

お礼

長文でのご回答、本当にありがとうございます。
誤差や精度のお話、とても参考になりました。
今後もっと注意してプログラムを組みたいと思います。
しかし、今回の縞模様はそのような誤差によるものでは無く、
バタフライ演算の過程でのミスによるものである
と思われるのですがどうでしょうか。
周期的に模様が現れていることからそう解釈しています。