I387 FSQRT命令で正しい四捨五入を得る方法はありますか？

StackOverflow https://stackoverflow.com//questions/9678224

質問

i387 FSQRT命令で正しい四捨五入を得る方法はありますか？...

... x87コントロールワードの精密モードを変更することをたどり取られています。 SQRT操作が中断された場合は間違っています。

私が扱っている問題は次のとおりです。x87 fsqrt Opcodeは、FPUレジスタの精度で正しく丸みを帯びた（IEEE 754あたり）平方根操作を実行します。 ）精度。ただし、効率的なシングルと倍精度の平方根関数を正しく丸めて（現在の丸めモードごとに）効率的な単一の正確な正方形の根関数を実装したいです。結果が過度の精度を有するので、結果を1つまたは2倍の精度に変換する第2ステップは再び正しく丸みを帯びた結果を残す。

いくつかの操作では、これをバイアスで回避することが可能です。例えば、2つの倍数の形でバイアスを追加することで、63ビット拡張精度の仮定の最後の52ビットに強制的にバイアスを加えることで、追加の精度を避けることができます。 。しかし、私は平方根でそのようなトリックをするための明らかな方法を見ません。

任意の賢い考え？

（意図されたアプリケーションはC sqrtおよびsqrtf関数の実装であるため、cもタグ付けされています。）

役に立ちましたか？

解決

最初に、邪魔にならないように明らかにしましょう：x87の代わりにsseを使用する必要があります。 SSE sqrtssおよびsqrtsdの指示は、まさにあなたが望むものを正確にしていますが、最近のすべてのx86システムでサポートされており、もっと速く速いです。

今、あなたがx87を使うことを主張するならば、私は良いニュースから始めます：あなたはフロートのために何かをする必要はありません。 P-BIT浮動小数点形式で正しく丸みを帯びた正方根を計算するには、2p + 2ビットが必要です。 80 > 2*24 + 2、単精度への追加の四捨五入は常に正しく丸くなり、正しく丸みを帯びた平方根があります。

今悪いニュース：80 < 2*53 + 2、倍精度のような運もできません。私はいくつかの回避策を提案することができます。これが私の頭の上から素敵な簡単なものです。

  1. y = round_to_double(x87_square_root(x));
  2. Dekker（ヘッドテール）製品を使用して、aとbを正確に計算します。
  3. 残差y*y = a + bを計算します。
  4. r = x - a - b
  5. if (r == 0) return y、if (r > 0)を使用し、y1 = y + 1 ulp、a1 S.T。 b1。 y1*y1 = a1 + b1とr1 = x - a1 - b1を比較し、残差が小さい（または残差が大きさが大きい場合は、ゼロのゼロが大きい場合）。
  6. rは、yに同じことをします。

    7. この手順はデフォルトの丸めモードのみを処理します。ただし、指示された丸めモードでは、目的のフォーマットに丸めるだけで正しいものがあります。

他のヒント

OK、私は私がより良い解決策を持っていると思います：

  1. 拡張精度（y=sqrt(x)）でfsqrtを計算します。
  2. 最後の11ビットが0x400ではない場合は、単に倍精度に変換して戻ります。
  3. 拡張精度表現の低ワードに0x100-(fpu_status_word&0x200)を追加します。
  4. 倍精度に変換して戻る

    5. ステップ3は、fsqrtの結果が切り上げられた場合に限り、ステータスワードのC1ビット（0x200）が1であるという事実に基づいています。ステップ2のテストにより、xが完全な正方形ではなかったため有効です。完璧な正方形だった場合、yは倍精度を超えてビットがないでしょう。

    は、ビット表現とリロードに取り組むのではなく、コンディションの浮動小数点でステップ3を実行するのが速い場合があります。

    これはコード（すべての場合でも機能するようです）：

    sqrt:
    fldl 4(%esp)
    fsqrt
    fstsw %ax
    sub $12,%esp
    fld %st(0)
    fstpt (%esp)
    mov (%esp),%ecx
    and $0x7ff,%ecx
    cmp $0x400,%ecx
    jnz 1f
    and $0x200,%eax
    sub $0x100,%eax
    sub %eax,(%esp)
    fstp %st(0)
    fldt (%esp)
1:  add $12,%esp
    fstpl 4(%esp)
    fldl 4(%esp)
    ret
    .

あなたが望むものではないかもしれませんが、387 fsqrt命令を利用していませんが、 glibc 32ビット整数演算で実装されています。それはNAN、INFS、下書きを正しく処理します - 実際のx87命令/ FP制御ワードフラグでこれらのチェックのいくつかを排除することが可能かもしれません。参照：glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

dbl-64/e_sqrt.cコードはそれほどフレンドリーではありません。どの仮定が一目で行われているかを知るのは難しいです。不思議なことに、ライブラリのI386 sqrt[f|l]実装はfsqrtを呼び出しますが、値を異なる方法でロードします。DP用のflds用のfldl

ライセンス： CC-BY-SA と 帰属
所属していません StackOverflow
scroll top