- C言語なら俺に聞け 146
851 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 00:06:54.68 ID:N9ICkOCi0 - 10000進数多倍長
超単純なFFT
演算は乗算と加算のみ
誤差の感じから100000進数でも大丈夫そうですね
計算式は基本以下を多倍長にしただけ
多少の無駄は除いてますが
----
uint64_t f(uint64_t n){
n++;
uint64_t a = 1;
uint64_t b = 0;
uint64_t t;
for (int i = 0 ; i < 64 ; i++){
t = b * b;
b = 2 * a * b + t;
a = a * a + t;
if (n & 0x8000000000000000){
t = b;
b = a + b;
a = t;
}
n += n;
}
return a;
}
|
- C言語なら俺に聞け 146
852 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 00:16:27.47 ID:N9ICkOCi0 - コードをアップしようと思ったけど
ideoneだとうまく動かないみたい
|
- C言語なら俺に聞け 146
884 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 21:43:12.30 ID:N9ICkOCi0 - F(1千万)の計算時間
gmp_fib_uiは0.72秒
私のは0.25秒
私が適当に作ったヤツの方が勝っちゃいましたね
コードは以下
https://ideone.com/hklTK1
|
- C言語なら俺に聞け 146
886 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 21:50:41.32 ID:N9ICkOCi0 - ローカル環境では動いてますが、
>>852に書いた通りideoneだと動きません
|
- C言語なら俺に聞け 146
894 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 22:42:41.92 ID:N9ICkOCi0 - >>884のバグを修正しました
https://ideone.com/4a3zE8
F(1千万)のideoneでの計算時間は0.31s
|
- C言語なら俺に聞け 146
895 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 23:13:57.75 ID:N9ICkOCi0 - n log nより遅いように見えるのは
サイズが大きいとキャッシュから外れるため
フーリエ変換はメモリのランダムアクセスが頻発するので
キャッシュに入るかどうかで大きく時間が変わる
適当に並び変えながら変換することで
キャッシュヒット率を上げたり
メインメモリに収まりきらない巨大サイズをHDDを使いながら変換することも出来るのだが、
今回はシンプルさ重視の為見送り
世の中にある高速FFTライブラリに差し替えれば
もっと高速化する
|
- C言語なら俺に聞け 146
896 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 23:25:51.63 ID:N9ICkOCi0 - >>887
環境依存
float, double がそれぞれIEEE754のbinary32とbinary64準拠の場合はおおよそ
-3.4028234663852885981*10^38 ≦ float ≦ 3.4028234663852885981*10^38
-1.7976931348623157081*10^308 ≦ double ≦ 1.7976931348623157081*10^308
正確には以下の範囲
float : ±(2^24-1)*2^104
double : ±(2^53-1)*2^971
これに追加して±∞, NaN などがある
組み込みなどではdoubleも32bitであることも多い
|
- C言語なら俺に聞け 146
897 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 23:33:36.42 ID:N9ICkOCi0 - C言語スレだから ^ は排他的論理和の意味で使うべきでしたかね
-340282346638528859811704183484516925440 ≦ binary32 ≦ 340282346638528859811704183484516925440
-1797693134862315708145274237317043567980705675258449965989174768
0315726078002853876058955863276687817154045895351438246423432132
6889464182768467546703537516986049910576551282076245490090389328
9440758685084551339423045832369032229481658085593321233482747978
26204144723168738177180919299881250404026184124858368
≦ binary64 ≦
1797693134862315708145274237317043567980705675258449965989174768
0315726078002853876058955863276687817154045895351438246423432132
6889464182768467546703537516986049910576551282076245490090389328
9440758685084551339423045832369032229481658085593321233482747978
26204144723168738177180919299881250404026184124858368
|
- C言語なら俺に聞け 146
898 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 23:41:42.36 ID:N9ICkOCi0 - >>888
単純なバグでした
size_t が64bit前提の部分があったのと、
非常に恥ずかしいですが以下のような動作不定なコードがありました
a->data[i++] = b->data[i] + carry;
私のローカル環境だとたまたま動いていたようです
intが16bitだとまだヤバそうですが、
16bit環境で動かすことはないと思うので気にしないことにします
|
- C言語なら俺に聞け 146
900 ：デフォルトの名無しさん (ﾜｯﾁｮｲ 0b50-2km2)[sage]：2018/08/11(土) 23:56:09.61 ID:N9ICkOCi0 - >>899
986.0033 はdoubleのリテラル
986.0033f はfloatのリテラル
|