ダブルタイプのIEEE754準拠のsqrt（）実装
double __ieee754_sqrt(double x)ハードウェア命令を使用して最初の近似値を取得する関数を実装しようとしています。
double __ieee754_sqrt(double x) {
double z;
/* get reciprocal of the square root (6.75 bits accuracy) */
__asm(" QSEED.DF %0,%1 \n": "=e" (z):"e" (x):);
z = 1 / z;
z = ( z + x / z) / 2; /* 1st Newton-Raphson iteration */
z = ( z + x / z) / 2; /* 2nd Newton-Raphson iteration */
z = ( z + x / z) / 2; /* 3rd Newton-Raphson iteration */
z = ( z + x / z) / 2; /* 4th Newton-Raphson iteration */
return z;
}
ただし、paranoia.c（link、link）テストは次のように文句を言います。
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.
質問：追加のロジックを実装する方法はchopping and correct rounding？
UPD。ハードウェアはネイティブにをサポートしていませんsqrt()。ハードウェアは、平方根の逆数の取得のみをサポートします（6.75ビットの精度）。
UPD2。
- 使用済みnjuffaのソリューション（多くの感謝！）マイナーチェンジ：使用
qseeddf()の代わりにqseedf()=>使用fma()の代わりにfmaf()。どうして？double<=>float変換が省略されるため、より高速になります。 - はい、融合積和命令（FMA）はハードウェアでサポートされています。
- ディスカッションに参加してくれて、そして詳細な答えをくれてありがとう！
- このトピックに関心のあるすべての人に、
sqrt()実装のリストを以下に示します。- Cygwin数学から。ライブラリ（
libm）::cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c著作権で保護されていCopyright (C) 1993 by Sun Microsystemsます。 - GNU Cライブラリから（
glibc）：glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c：タイトル付きIBM Accurate Mathematical Library。glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c：__builtin_fma()関数を使用します。
- Cygwin数学から。ライブラリ（
回答
独自の実装の構築に着手する前に、インターネットを検索して、適切で十分にテストされたオープンソースコードが利用可能かどうかを確認することをお勧めします。
一般的な反復アルゴリズムでは、逆数の平方根に対して除算のない反復を使用して目的の精度にし、引数を逆乗算して平方根を計算し、最後に目的の丸めモードを使用して丸めます。逆平方根の反復では、2次収束を伴うニュートンラプソン反復（正しいビット数を約2倍）または3次収束を伴うハレー反復（正しいビット数を約3倍）のいずれかを使用できます。高次の反復は存在しますが、通常は使用されません。
コードを単純に保つために、バイナリ浮動小数点演算の場合は、引数を2つの連続するビネードで構成される単一の狭い区間に減らすことをお勧めします。指数操作が必要なため、これは通常、最高のパフォーマンスの実装にはならないことに注意してください。パフォーマンス上の理由から、倍精度実装の最初の反復は単精度で実行されることがよくあります。
以下の例示的なISO-C99実装では、正しく丸められた倍精度平方根をこれらの線に沿って実装する方法を示しています。私はと仮定していますdoubleIEEE-754へのマップbinary64とそのfloatIEEE-754にマップをbinary32。sqrtIEEE-754の最も近いまたは偶数の丸めモードで実装されるように制限しています。
非常に重要なことは、プロセスハードウェアが融合された積和命令を提供し、これらが標準の数学ライブラリ関数fmafとを介して正しく公開されていることを前提としていますfma。コメントで、私はFMAの利用可能性についてOPに説明を求めていましたが、フィードバックが利用可能になる前にコードを開始することにしました。 FMAを使用しない実装は可能ですが、はるかに困難であり、十分に完全な処理は、Stackoverflowの回答のスペースを超える可能性があります。
OPはターゲットアーキテクチャを指定せず、開始近似の詳細を提供しなかったため、すべての非例外引数が削減される区間[0.25、1]での多項式ミニマックス近似に基づいて、以下の独自の開始近似を使用しています。qseedf()結果は約7ビットの精度であるため、OPの組み込みハードウェアよりもわずかに優れています。この違いが重要かどうか、私は評価できません。
アルゴリズム、特に丸めロジックは、Peter Marksteinのアイデアに依存しているため、アルゴリズムが構造上正しいと合理的に確信しています。ここでは、非常に基本的なテストのみを実装しました。業界のベストプラクティスは、そのようなアルゴリズムの正しさを数学的に証明することです。たとえば、DavidRussinoffとJohnHarrisonの出版物を参照してください。ピンチでは、すべてのビネードを実行するランダムなパターンベースのテストと組み合わせて、2つの連続するビネードにわたる徹底的なテスト（最近は小さなクラスターを数日間実行することで実現可能）で逃げることができるかもしれません。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
/* Approximate 1/sqrt(a) on [0.25, 1] with an accuracy of about 7 bits */
float qseedf (float a)
{
float r;
r = -2.43845296f;
r = fmaf (r, a, 6.22994471f);
r = fmaf (r, a, -5.91090727f);
r = fmaf (r, a, 3.11237526f);
return r;
}
double my_sqrt (double a)
{
const double QNAN_INDEFINITE = 0.0 / 0.0;
const double half = 0.5;
const double three_eighth = 0.375;
double refined_rsqrt_approx, sqrt_approx, sqrt_residual, result, b;
double rsqrt_approx, rsqrt_approx_err, rsqrt_approx_squared, reduced_arg;
float argf, approxf, approxf_err;
int e, t, f;
/* handle normal cases */
if ((a >= 0) && (a < INFINITY)) {
/* compute exponent adjustments */
b = frexp (a, &e);
t = e - 2*512;
f = t / 2;
t = t - 2 * f;
f = f + 512;
/* map argument into the primary approximation interval [0.25,1) */
reduced_arg = ldexp (b, t);
/* Compute initial low-precision approximation */
argf = (float)reduced_arg;
approxf = qseedf (argf);
/* Apply two Newton-Raphson iterations with quadratic convergence */
approxf_err = fmaf (-argf, approxf * approxf, 1.0f);
approxf = fmaf (0.5f * approxf, approxf_err, approxf);
approxf_err = fmaf (-argf, approxf * approxf, 1.0f);
approxf = fmaf (0.5f * approxf, approxf_err, approxf);
/* rsqrt approximation is now accurate to 1 single-precision ulp */
rsqrt_approx = (double)approxf;
/* Perform a Halley iteration wih cubic convergence. Based on the work
of Peter Markstein. See: Peter Markstein, "IA-64 and Elementary
Functions", Prentice Hall 2000
*/
rsqrt_approx_squared = rsqrt_approx * rsqrt_approx;
rsqrt_approx_err = fma (-reduced_arg, rsqrt_approx_squared, 1.0);
refined_rsqrt_approx = fma (fma (rsqrt_approx_err, three_eighth, half),
rsqrt_approx * rsqrt_approx_err, rsqrt_approx);
sqrt_approx = reduced_arg * refined_rsqrt_approx;
sqrt_residual = fma (-sqrt_approx, sqrt_approx, reduced_arg);
result = fma (sqrt_residual, half * refined_rsqrt_approx, sqrt_approx);
/* map back from primary approximation interval by jamming exponent */
result = ldexp (result, f);
} else {
/* handle special cases */
result = (a < 0) ? QNAN_INDEFINITE : (a + a);
}
return result;
}
/*
https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
From: geo <[email protected]>
Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
Subject: 64-bit KISS RNGs
Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)
This 64-bit KISS RNG has three components, each nearly
good enough to serve alone. The components are:
Multiply-With-Carry (MWC), period (2^121+2^63-1)
Xorshift (XSH), period 2^64-1
Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64 (kiss64_t = (kiss64_x << 58) + kiss64_c, \
kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64 (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
kiss64_y ^= (kiss64_y << 43))
#define CNG64 (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
int main (void)
{
const uint64_t N = 10000000000ULL; /* desired number of test cases */
double arg, ref, res;
uint64_t argi, refi, resi, count = 0;
double spec[] = {0, 1, INFINITY, NAN};
printf ("test a few special cases:\n");
for (int i = 0; i < sizeof (spec)/sizeof(spec[0]); i++) {
printf ("my_sqrt(%22.13a) = %22.13a\n", spec[i], my_sqrt(spec[i]));
printf ("my_sqrt(%22.13a) = %22.13a\n", -spec[i], my_sqrt(-spec[i]));
}
printf ("test %llu random cases:\n", N);
do {
count++;
argi = KISS64;
memcpy (&arg, &argi, sizeof arg);
res = my_sqrt (arg);
ref = sqrt (arg);
memcpy (&resi, &res, sizeof resi);
memcpy (&refi, &ref, sizeof refi);
if (resi != refi) {
printf ("\rerror @ arg=%22.13a res=%22.13a ref=%22.13a\n",
arg, res, ref);
return EXIT_FAILURE;
}
if ((count & 0xfffff) == 0) printf ("\r[%llu]", count);
} while (count < N);
printf ("\r[%llu]", count);
printf ("\ntests PASSED\n");
return EXIT_SUCCESS;
}
上記のプログラムの出力は次のようになります。
test a few special cases:
my_sqrt( 0x0.0000000000000p+0) = 0x0.0000000000000p+0
my_sqrt( -0x0.0000000000000p+0) = -0x0.0000000000000p+0
my_sqrt( 0x1.0000000000000p+0) = 0x1.0000000000000p+0
my_sqrt( -0x1.0000000000000p+0) = -0x1.#IND000000000p+0
my_sqrt( 0x1.#INF000000000p+0) = 0x1.#INF000000000p+0
my_sqrt( -0x1.#INF000000000p+0) = -0x1.#IND000000000p+0
my_sqrt( 0x1.#QNAN00000000p+0) = 0x1.#QNAN00000000p+0
my_sqrt( -0x1.#QNAN00000000p+0) = -0x1.#QNAN00000000p+0
test 10000000000 random cases:
[10000000000]
tests PASSED
z = 1 / z;
z = ( z + x / z) / 2; /* 1st Newton-Raphson iteration */
...
->
z = 1 / z;
z += ( x / z - z) * 0.5; /* 1st Newton-Raphson iteration */
...
これは速いかもしれません。
そして、1回の反復をより早く停止します（私は思います）。
停止したら、比較z*zしてx。z*z（と思う）よりも小さくないだろうx。から1ulpをサブトレースし、vszを確認します。これは「正しい丸め」の完全なチェックではありませんが、との間で決定するのに「十分」である可能性があります。z*zxzz - 1ulp
非常に広範囲のエラーが発生したため、丸め、さらには精度に関しては、浮動小数点の残りの「ハードウェア」がずさんなものになるのではないかと心配しています。
おっと、忘れました。あなたに近似を与える理由がありました1/z-近似1 / zを続けてください。除算の代わりに乗算を使用してそれを行うことができるため、（ほとんどのハードウェアで）大幅に高速になり、場合によっては丸めが少なくなります。
z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;
また、の乗算を行う代わりに、指数をデクリメントする方法があるかどうかを確認してください/ 2。