これは最初に提起された質問への回答ではありませんが、同様の問題に取り組んでいるプログラマーであれば、この回答が役立つかもしれません。
認識された困難がどこにあるのか、私には本当にわかりません。厳密な IEEE-754 binary64 セマンティクスを提供しながら、80387 浮動小数点演算に制限され、80 ビット長の double 計算を保持することは、GCC-4.6.3 と clang-3.0 (LLVM に基づく) の両方で明確に指定された C99 キャスト規則に従っているようです。 3.0)。
編集して追加: それでも、Pascal Cuoq は正しい: gcc-4.6.3 または clang-llvm-3.0 のどちらも、実際には '387 浮動小数点演算に対してこれらの規則を正しく適用しません。適切なコンパイラ オプションを指定すると、ルールはコンパイル時に評価される式に正しく適用されますが、実行時の式には適用されません。以下の休憩の後にリストされている回避策があります。
私は分子動力学シミュレーションのコードを作成しており、再現性/予測可能性の要件と、可能な限り最大の精度を維持したいという欲求に精通しているため、ここで話していることを知っていると主張します. この回答は、ツールが存在し、使いやすいことを示す必要があります。問題は、これらのツールを認識していない、または使用していないことから生じます。
(私が気に入っている好ましい例は、Kahan 加算アルゴリズムです。C99 と適切なキャスト (Wikipedia のサンプル コードなどにキャストを追加する) を使用すると、トリックや追加の一時変数はまったく必要ありません。実装は、コンパイラの最適化レベルに関係なく機能します。-O3と-Ofast。）
C99 では、キャストと代入の両方がすべての余分な範囲と精度を削除することを (たとえば 5.4.2.2 で) 明示的に述べています。これは、long double計算中に使用される一時変数を として定義しlong double、入力変数をその型にキャストすることで算術を使用できることを意味します。IEEE-754 binary64 が必要なときはいつでも、 にキャストするだけdoubleです。
'387 では、キャストは上記の両方のコンパイラで代入とロードを生成します。これにより、80 ビット値が IEEE-754 binary64 に正しく丸められます。このコストは、私の意見では非常に妥当です。かかる正確な時間は、アーキテクチャと周囲のコードによって異なります。通常はそうであり、他のコードとインターリーブして、コストを無視できるレベルに下げることができます。MMX、SSE、または AVX が使用可能な場合、それらのレジスタは 80 ビットの 80387 レジスタから分離され、キャストは通常、値を MMX/SSE/AVX レジスタに移動することによって行われます。
(一時変数に特定の浮動小数点型などを使用する製品コードを好むので、必要なアーキテクチャと速度/精度のトレードオフに応じてtempdouble定義できますdouble。 )long double
手短に：
すべての変数とリテラル定数が正確だからといって(expression)、が正確であると仮定しないでください。正確な結果が必要であるかのようにdouble記述します。(double)(expression)double
これは複合式にも当てはまり、多くのレベルのキャストを持つ扱いにくい式になることがあります。
80 ビットの精度で計算したいが、最初に 64 ビットに丸めたそれぞれの積も必要な場合は、次を使用しますexpr1。expr2
long double expr1;
long double expr2;
double product = (double)(expr1) * (double)(expr2);
product2 つの 64 ビット値の積として計算されることに注意してください。80 ビット精度で計算されず、切り捨てられます。80 ビット精度で積を計算してから切り捨てると、次のようになります。
double other = expr1 * expr2;
または、何が起こっているかを正確に伝える説明的なキャストを追加します。
double other = (double)((long double)(expr1) * (long double)(expr2));
productとotherがしばしば異なることは明らかです。
C99 のキャスト規則は、32 ビット/64 ビット/80 ビット/128 ビットの浮動小数点値が混在している場合に習得しなければならないもう 1 つのツールです。float実際、binary32 と binary64 の float を混在させると (そしてdoubleほとんどのアーキテクチャで)まったく同じ問題が発生します!
Pascal Cuoq の探索コードを書き直して、キャスト規則を正しく適用すると、これがより明確になるのでしょうか?
#include <stdio.h>
#define TEST(eq) printf("%-56s%s\n", "" # eq ":", (eq) ? "true" : "false")
int main(void)
{
double d = 1.0 / 10.0;
long double ld = 1.0L / 10.0L;
printf("sizeof (double) = %d\n", (int)sizeof (double));
printf("sizeof (long double) == %d\n", (int)sizeof (long double));
printf("\nExpect true:\n");
TEST(d == (double)(0.1));
TEST(ld == (long double)(0.1L));
TEST(d == (double)(1.0 / 10.0));
TEST(ld == (long double)(1.0L / 10.0L));
TEST(d == (double)(ld));
TEST((double)(1.0L/10.0L) == (double)(0.1));
TEST((long double)(1.0L/10.0L) == (long double)(0.1L));
printf("\nExpect false:\n");
TEST(d == ld);
TEST((long double)(d) == ld);
TEST(d == 0.1L);
TEST(ld == 0.1);
TEST(d == (long double)(1.0L / 10.0L));
TEST(ld == (double)(1.0L / 10.0));
return 0;
}
GCC と clang の両方を使用した出力は次のとおりです。
sizeof (double) = 8
sizeof (long double) == 12
Expect true:
d == (double)(0.1): true
ld == (long double)(0.1L): true
d == (double)(1.0 / 10.0): true
ld == (long double)(1.0L / 10.0L): true
d == (double)(ld): true
(double)(1.0L/10.0L) == (double)(0.1): true
(long double)(1.0L/10.0L) == (long double)(0.1L): true
Expect false:
d == ld: false
(long double)(d) == ld: false
d == 0.1L: false
ld == 0.1: false
d == (long double)(1.0L / 10.0L): false
ld == (double)(1.0L / 10.0): false
ただし、最近のバージョンの GCC では、 の右側がld == 0.1最初に long double に (つまり にld == 0.1L) 昇格され、 が生成されます。trueまた、SSE/AVX ではlong double128 ビットになります。
純粋な '387 テストでは、使用しました
gcc -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test
clang -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test
、、、、、、...など-fomit-frame-pointerのさまざまな最適-O0化フラグ-O1の組み合わせ。-O2-O3-Os
他のフラグまたは C99 コンパイラを使用しても、long doubleサイズ (およびld == 1.0現在の GCC バージョン) を除いて、同じ結果になるはずです。違いがある場合は、お知らせいただければ幸いです。そのようなコンパイラ (コンパイラのバージョン) についてユーザーに警告する必要があるかもしれません。Microsoft は C99 をサポートしていないので、私にはまったく興味がないことに注意してください。
Pascal Cuoq は、以下の一連のコメントで興味深い問題を提起していますが、私はすぐには気づきませんでした。
式を評価するとき、GCC と clang の両方で-mfpmath=387、すべての式が 80 ビット精度を使用して評価されるように指定します。これは、例えば
7491907632491941888 = 0x1.9fe2693112e14p+62 = 110011111111000100110100100110001000100101110000101000000000000
5698883734965350400 = 0x1.3c5a02407b71cp+62 = 100111100010110100000001001000000011110110111000111000000000000
7491907632491941888 * 5698883734965350400 = 42695510550671093541385598890357555200 = 100000000111101101101100110001101000010100100001011110111111111111110011000111000001011101010101100011000000000000000000000000
バイナリ結果の中央にある 1 の文字列は、53 ビットと 64 ビットの仮数 (それぞれ 64 ビットと 80 ビットの浮動小数点数) の差にあるため、誤った結果が得られます。したがって、期待される結果は
42695510550671088819251326462451515392 = 0x1.00f6d98d0a42fp+125 = 100000000111101101101100110001101000010100100001011110000000000000000000000000000000000000000000000000000000000000000000000000
だけで得られた結果-std=c99 -m32 -mno-sse -mfpmath=387は
42695510550671098263984292201741942784 = 0x1.00f6d98d0a43p+125 = 100000000111101101101100110001101000010100100001100000000000000000000000000000000000000000000000000000000000000000000000000000
理論的には、オプションを使用して正しい C99 丸めルールを強制するように gcc と clang に指示できるはずです。
-std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard
ただし、これはコンパイラが最適化する式にのみ影響し、387 の処理をまったく修正していないようです。たとえばclang -O1 -std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard test.c -o test && ./test、Pascal Cuoq のサンプル プログラムtest.cを使用すると、IEEE-754 ルールに従って正しい結果が得られますが、これはコンパイラが 387 をまったく使用せずに式を最適化するためです。
簡単に言えば、計算の代わりに
(double)d1 * (double)d2
gcc と clang の両方が実際に '387 に計算を指示します
(double)((long double)d1 * (long double)d2)
これはまさしくこれは、gcc-4.6.3 と clang-llvm-3.0 の両方に影響するコンパイラのバグであり、簡単に再現できるバグだと思います。FLT_EVAL_METHOD=2(Pascal Cuoqは、倍精度引数の操作が常に拡張精度で行われることを意味すると指摘していますがlibm、'387 の一部を書き直さなければならないことを除けば、C99 でそれを行い、 IEEE を考慮する正当な理由は見当たりません- 754 の規則はハードウェアで達成可能です! 結局のところ、'387 制御語を式の精度に合わせて変更することにより、正しい操作はコンパイラで簡単に達成できます. そして、この動作を強制するコンパイラ オプションが与えられた場合-- -std=c99 -ffloat-store -fexcess-precision=standard- - 意味がないFLT_EVAL_METHOD=2適切なコンパイラ フラグが指定されている場合、コンパイル時に評価される式は正しく評価され、実行時に評価される式のみが正しくない結果になることに注意することが重要です。
最も簡単で移植可能な回避策は、fesetround(FE_TOWARDZERO)(from fenv.h) を使用してすべての結果をゼロに丸めることです。
場合によっては、ゼロに丸めることで、予測可能性や病理学的なケースに役立つことがあります。特に、 のような間隔の場合x = [0,1)、ゼロに向かって丸めることは、丸めによって上限に達することがないことを意味します。区分スプラインなどを評価する場合に重要です。
他の丸めモードでは、387 ハードウェアを直接制御する必要があります。
__FPU_SETCW()from#include <fpu_control.h>またはオープンコードのいずれかを使用できます。たとえば、precision.c次のとおりです。
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#define FP387_NEAREST 0x0000
#define FP387_ZERO 0x0C00
#define FP387_UP 0x0800
#define FP387_DOWN 0x0400
#define FP387_SINGLE 0x0000
#define FP387_DOUBLE 0x0200
#define FP387_EXTENDED 0x0300
static inline void fp387(const unsigned short control)
{
unsigned short cw = (control & 0x0F00) | 0x007f;
__asm__ volatile ("fldcw %0" : : "m" (*&cw));
}
const char *bits(const double value)
{
const unsigned char *const data = (const unsigned char *)&value;
static char buffer[CHAR_BIT * sizeof value + 1];
char *p = buffer;
size_t i = CHAR_BIT * sizeof value;
while (i-->0)
*(p++) = '0' + !!(data[i / CHAR_BIT] & (1U << (i % CHAR_BIT)));
*p = '\0';
return (const char *)buffer;
}
int main(int argc, char *argv[])
{
double d1, d2;
char dummy;
if (argc != 3) {
fprintf(stderr, "\nUsage: %s 7491907632491941888 5698883734965350400\n\n", argv[0]);
return EXIT_FAILURE;
}
if (sscanf(argv[1], " %lf %c", &d1, &dummy) != 1) {
fprintf(stderr, "%s: Not a number.\n", argv[1]);
return EXIT_FAILURE;
}
if (sscanf(argv[2], " %lf %c", &d2, &dummy) != 1) {
fprintf(stderr, "%s: Not a number.\n", argv[2]);
return EXIT_FAILURE;
}
printf("%s:\td1 = %.0f\n\t %s in binary\n", argv[1], d1, bits(d1));
printf("%s:\td2 = %.0f\n\t %s in binary\n", argv[2], d2, bits(d2));
printf("\nDefaults:\n");
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nExtended precision, rounding to nearest integer:\n");
fp387(FP387_EXTENDED | FP387_NEAREST);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nDouble precision, rounding to nearest integer:\n");
fp387(FP387_DOUBLE | FP387_NEAREST);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nExtended precision, rounding to zero:\n");
fp387(FP387_EXTENDED | FP387_ZERO);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nDouble precision, rounding to zero:\n");
fp387(FP387_DOUBLE | FP387_ZERO);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
return 0;
}
clang-llvm-3.0 を使用してコンパイルおよび実行すると、正しい結果が得られます。
clang -std=c99 -m32 -mno-sse -mfpmath=387 -O3 -W -Wall precision.c -o precision
./precision 7491907632491941888 5698883734965350400
7491907632491941888: d1 = 7491907632491941888
0100001111011001111111100010011010010011000100010010111000010100 in binary
5698883734965350400: d2 = 5698883734965350400
0100001111010011110001011010000000100100000001111011011100011100 in binary
Defaults:
Product = 42695510550671098263984292201741942784
0100011111000000000011110110110110011000110100001010010000110000 in binary
Extended precision, rounding to nearest integer:
Product = 42695510550671098263984292201741942784
0100011111000000000011110110110110011000110100001010010000110000 in binary
Double precision, rounding to nearest integer:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
Extended precision, rounding to zero:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
Double precision, rounding to zero:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
fp387()つまり、 を使用して精度と丸めモードを設定することで、コンパイラの問題を回避できます。
欠点は、一部の数学ライブラリ ( libm.a、libm.so) が、中間結果が常に 80 ビット精度で計算されるという前提で記述されている可能性があることです。少なくともfpu_control.hx86_64 の GNU C ライブラリには、「libm には拡張精度が必要です」というコメントがあります。libm幸いなことに、機能が必要な場合は、GNU C ライブラリなどから '387 実装を取得して、それらをヘッダー ファイルに実装するか、動作することがわかっている を記述することができますmath.h。実際、私はそこを助けることができるかもしれないと思います。