0

次の関数を検討してください。

#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>

template <typename Type>
inline Type a(const Type dx, const Type a0, const Type z0, const Type b1)
{
    return (std::sqrt(std::abs(2*b1-z0))*dx)+a0;
}

template <typename Type>
inline Type b(const Type dx, const Type a0, const Type z0, const Type a1)
{
    return (std::pow((a1-a0)/dx, 2)+ z0)/2;
}

int main(int argc, char* argv[])
{
    double dx = 1.E-6;
    double a0 = 1;
    double a1 = 2;
    double z0 = -1.E7;
    double b1 = -10;
    std::cout<<std::scientific;
    std::cout<<std::setprecision(std::numeric_limits<double>::digits10);
    std::cout<<a1-a(dx, a0, z0, b(dx, a0, z0, a1))<<std::endl;
    std::cout<<b1-b(dx, a0, z0, a(dx, a0, z0, b1))<<std::endl;
    return 0;
}

私のマシンでは、次のように返されます。

0.000000000000000e+00
-1.806765794754028e-07

(0, 0) の代わりに。2 番目の式には大きな丸め誤差があります。

私の質問は、型を変更せずに各関数の丸め誤差を減らす方法です (これらの 2 つの関数宣言を保持する必要があります (ただし、式は再配置できます): それらはより大きなプログラムからのものです)。

4

2 に答える 2

1

残念なことに、すべての浮動小数点型は丸め誤差で有名です。それなしでは 0.1 を格納することさえできません (これは手で長い除算を使用して証明できます: 2 進数に相当するものは 0b0.0001100110011001100... です)。その pow をハードコードされた乗算に拡張するなどの回避策を試すこともできますが、最終的には、丸め誤差の影響を予測して最小限に抑えるようにプログラムをコーディングする必要があります。ここにいくつかのアイデアがあります：

  • 浮動小数点値が等しいかどうかを比較しないでください。私が見たいくつかの代替比較には次のものがあります。この特定のテストのために。

  • 数値の長い配列をアキュムレータに追加しないでください。アキュムレータが大きくなると、丸め誤差により配列の末尾が完全に失われる可能性があります。Jason Sanders と Edward Kandrot による「Cuda by Example」では、1 要素の配列が得られるまで、各ステップで前のステップの半分のサイズの配列が生成されるように、要素の各ペアを個別に再帰的に追加することを著者は推奨しています。

于 2013-11-14T01:05:32.207 に答える
0

a() では、sqrt()*dx の小さくて不正確な結果に a0 (正確に 1) を追加すると、精度が失われます。

関数 b() は、指定された値を使用して精度を失うことはありません。

2 番目の出力のように、b() の前に a() を呼び出すと、すでに不正確な数値に対して数学演算を実行しているため、エラーが悪化します。

浮動小数点エラーが発生する可能性が低い演算を最初に実行し、浮動小数点エラーが発生する可能性が高い演算を最後に実行するように、数学演算を構造化するようにしてください。

または、関数内で、それらが「long double」値で動作していることを確認してください。たとえば、次の例では、浮動小数点の昇格を使用して、最初の算術演算で double を long double に昇格します (演算子の優先順位に注意してください)。

template <typename Type>
inline Type a(const Type dx, const Type a0, const Type z0, const Type b1)
{
    return (std::sqrt(std::abs(2*static_cast<long double>(b1)-z0))*dx)+a0;
}
于 2013-11-14T00:53:29.147 に答える