8

私は現在、値の推定のために浮動小数点数値を強化しています。(これp(k,t)は: 興味のある人向けです。) 基本的に、このユーティリティはこの値を過小評価することはありません。可能性のある素数生成のセキュリティは、数値的に堅牢な実装に依存します。出力結果は公開された値と一致しますがDBL_EPSILON、特に除算によって真の値よりも決して小さい結果が得られないようにするために値を使用しました。

検討：double x, y; /* assigned some values... */

評価:r = x / y;頻繁に発生しますが、これらの (有限精度) 結果は真の結果から有効桁数を切り捨てる場合があります - おそらく無限精度の有理展開です。私は現在、分子にバイアスを適用することでこれを軽減しようとしています。つまり、

r = ((1.0 + DBL_EPSILON) * x) / y;

この主題について何か知っている場合、p(k,t)通常、ほとんどの見積もりよりもはるかに小さいですが、この「観察」で問題を却下するのに十分ではありません. もちろん、次のように言えます。

(((1.0 + DBL_EPSILON) * x) / y) >= (x / y)

もちろん、「偏った」結果が「正確な」値以上であることを確認する必要があります。操作またはスケーリングに関係していることは確かDBL_EPSILONですが、「バイアスされた」結果が「正確な」結果を最小値で超えることは明らかです-IEEE-754算術仮定の下で実証可能です。

はい、Goldberg の論文を調べて、堅牢なソリューションを探しました。丸めモードの操作を提案しないでください。理想的には、浮動小数点定理をよく理解しているか、非常によく説明された例を知っている人からの回答を求めています。


編集:明確にすること、(((1.0 + DBL_EPSILON) * x) / y)またはフォーム(((1.0 + c) * x) / y)は前提条件ではありません。これは、確固たる根拠を提供することなく、単に「おそらく十分だ」として使用していたアプローチでした. 分子と分母が特別な値 (NaN、Infs など) になることはなく、分母がゼロになることもありません。

4

1 に答える 1

10

最初：丸めモードを設定したくないことは知っていますが、他の人が指摘したように、精度の観点から、丸めモードを設定すると可能な限り良い答えが得られると本当に言われるべきです。x具体的には、とが両方とも正であると仮定するとy(そのように見えますが、質問では明示的に述べられていません)、以下は望ましい効果を持つ標準的な C スニペットです [1]:

#include <math.h>
#pragma STDC FENV_ACCESS on

int OldRoundingMode = fegetround();
fesetround(FE_UPWARD);
r = x/y;
fesetround(OldRoundingMode);

さて、それはさておき、丸めモードを変更したくない正当な理由があります (一部のプラットフォームでは、丸めモードを変更すると大きなシリアライズ ストールが発生するなど)、およびそうしたくないというあなたの願望は、さりげなく無視されるべきではありません。それで、あなたの質問を尊重して、他に何ができますか?

プラットフォームが融合乗加算をサポートしている場合は、非常に洗練されたソリューションを利用できます。

#include <math.h>
r = x/y;
if (fma(r,y,-x) < 0) r = nextafter(r, INFINITY);

ハードウェア fma をサポートするプラットフォームでは、これは非常に効率的です。fma( ) がソフトウェアで実装されている場合でも、許容される場合があります。このアプローチには、丸めモードを変更した場合と同じ結果が得られるという利点があります。つまり、可能な限り厳しい境界です。

プラットフォームの C ライブラリが旧式で を提供していない場合fmaでも、希望はあります。あなたの主張は正しいです（少なくとも非正規値がないと仮定します-非正規値で何が起こるかについてもっと考える必要があります）。(1.0+DBL_EPSILON)*x/y実際には、常に無限に正確な x/y 以上です。このプロパティの最小値よりも 1 ulp 大きい場合がありますが、これは非常に小さく、おそらく許容できるマージンです。これらの主張の証明はかなり面倒で、おそらく StackOverflow には適していませんが、簡単なスケッチを示します。

  1. デノーマルを無視すると、[1.0, 2.0) の x、y に制限するだけで十分です。

  2. (1.0 + eps)*x >= x + eps > x. これを確認するには、次のことを確認します。

    (1.0 + eps)*x = x + x*eps >= x + eps > x.

  3. P を数学的に正確な x/y とします。我々は持っています：

    (1.0 + eps)*x/y >= (x + eps)/y = x/y + eps/y = P + eps/y

    ここで、y は 2 で制限されているため、次のようになります。

    (1.0 + eps)*x/y > P + eps/2

    これは、結果が値 >= P に丸められることを保証するのに十分です。これは、より厳密な境界への道も示しています。多くの場合、代わりに使用nextafter(x,INFINITY)/yして、より厳密な境界で目的の効果を得ることができます。(nextafter(x,INFINITY)は常に x + ulp ですが、半分の時間は x + 2ulp になります。ライブラリ関数(1.0 + eps)*xの呼び出しを避けたい場合は、代わりに を使用して、正の通常値の作業仮定の下で同じ結果を得ることができます)。nextafter(x + (0.75*DBL_EPSILON)*x)

  1. 厳密に言えば、これはかなり複雑になります。このようなコードを実際に書く人はいませんが、次のようなコードになります。

    #include <math.h>
#pragma STDC FENV_ACCESS on

#if defined FE_UPWARD
    int OldRoundingMode = fegetround();
    if (OldRoundingMode < 0) goto Error;
    if (fesetround(FE_UPWARD)) goto Error;
    r = x/y;
    if (fesetround(OldRoundingMode)) goto TrulyHosed;
    return r;
TrulyHosed:
    // we established the desired rounding mode and did our computation,
    // but now we can't set it back to the original mode.  I have no idea
    // how you handle this gracefully.
Error:
#else
    // we can't establish the desired rounding mode, so fall back on
    // something else.
于 2012-05-10T14:47:59.377 に答える