9

i387 fsqrt 命令で正しい丸めを取得する方法はありますか?...

... x87 制御ワードで精度モードを変更することは別として- それが可能であることはわかっていますが、sqrt 操作が中断された場合に精度モードが間違ってしまう厄介な再入可能タイプの問題があるため、合理的な解決策ではありません。

私が扱っている問題は次のとおりです。x87fsqrtオペコードは、拡張 (80 ビット) 精度であると仮定する fpu レジスタの精度で、正しく丸められた (IEEE 754 ごとの) 平方根演算を実行します。ただし、(現在の丸めモードに従って) 結果が正しく丸められた効率的な単精度および倍精度の平方根関数を実装するために使用したいと考えています。結果の精度が過剰になるため、結果を単精度または倍精度に変換する 2 番目のステップで再び丸めが行われ、結果が正しく丸められない可能性があります。

一部の操作では、バイアスを使用してこれを回避できます。たとえば、倍精度値の 52 の有効ビットを 63 ビットの拡張精度仮数の最後の 52 ビットに強制する 2 の累乗の形でバイアスを追加することで、加算結果の過剰な精度を回避できます。 . しかし、平方根を使ってそのようなトリックを行う明確な方法はわかりません。

何か賢いアイデアはありますか？

(意図したアプリケーションが Csqrtとsqrtf関数の実装であるため、C もタグ付けされています。)

4

3 に答える 3

15

まず、明白なことから始めましょう。x87 の代わりに SSE を使用する必要があります。SSEsqrtssとsqrtsd命令はまさにあなたが望むものを実行し、最新のすべての x86 システムでサポートされており、同様に大幅に高速です。

さて、x87 を使用することに固執する場合は、朗報から始めましょう。float については何もする必要はありません。2p + 2p ビット浮動小数点形式で正しく丸められた平方根を計算するには、ビットが必要です。80 > 2*24 + 2、単精度への追加の丸めは常に正しく丸められ、正しく丸められた平方根があるためです。

ここで悪いニュース: 80 < 2*53 + 2、倍精度の場合、そのような運はありません。いくつかの回避策を提案できます。ここに私の頭の上から素敵な簡単なものがあります.

  1. させてy = round_to_double(x87_square_root(x));
  2. Dekker (head-tail) 積を使用して正確に計算aします。by*y = a + b
  3. 残差 を計算しr = x - a - bます。
  4. if (r == 0) return y
  5. if (r > 0)、 let y1 = y + 1 ulp、および compute a1、b1st y1*y1 = a1 + b1。と比較r1 = x - a1 - b1し、残差が小さい方 (残差の大きさが等しい場合は下位ビットがゼロの場合) に応じてまたはrを返します。yy1
  6. if (r < 0)、 についても同じことを行いy1 = y - 1 ulpます。

この手順では、デフォルトの丸めモードのみが処理されます。ただし、有向丸めモードでは、単純に目的の形式に丸めることで正しい処理が行われます。

于 2012-03-13T18:17:35.457 に答える
3

OK、私はより良い解決策があると思います:

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

fsqrtステップ 3 は、の結果が切り上げられた場合にのみ、ステータス ワードの 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
于 2012-03-15T04:00:00.663 に答える
0

fsqrt387命令を利用していないため、希望どおりではない場合がありますが、32ビット整数演算で実装さsqrtf(float)れたglibcには驚くほど効率的です。NaN、Infs、非正規化数も正しく処理します。実際のx87命令/FP制御ワードフラグを使用してこれらのチェックの一部を排除できる可能性があります。見る：glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

dbl-64/e_sqrt.cコードはそれほどフレンドリーではありません。一目でどのような仮定がなされているのか見分けるのは難しい。不思議なことに、ライブラリのi386sqrt[f|l]実装は単にを呼び出すだけfsqrtですが、値のロードは異なります。fldsSPの場合fldl、DPの場合。

于 2012-03-13T17:36:52.230 に答える