algorithm - 是否可以使用 32 位平方根的函数来帮助计算 64 位平方根？

Question

为了扩展这个想法，假设我有 2 个 32 位寄存器，分别代表 64 位 float 的高位和低位。我想计算它们的 64 位平方根。但是，虽然我没有 64 位平方根函数，但我有 32 位平方根函数。

我的问题是:如果我想计算 64 位平方根，我可以使用 32 位平方根对我有帮助吗？类似的东西？

Answer 1

TL;DR 是。

根据您平台的硬件、工具链和数学库的功能和缺陷，这可能不一定是计算 double 平方根的最快或最不痛苦的方法。下面我展示了一种基于 Arnold Schönhage 的平方根和倒数平方根的耦合迭代的直接方法:

从平方根倒数的近似值 r_approx ~= 1/√a 开始，我们计算 s₀ = a * r_approx 和 r₀ = r_approx/2，然后迭代:

s_i+1 = s_i + r_i * (a - s_i * s<子>我 )
r_i+1 = r_i + r_i * (1 - r_i * 2 * s_>i+1)

其中 s_i 是 √a 的近似值，r_i 是 1/(2√a) 的近似值。这个迭代是 Newton-Raphson 迭代巧妙地重新安排，因此具有二次收敛，这意味着每一步将大约加倍正确的位数。从单精度 r_approx 开始，只需两步即可达到 double 精度。

如果我们现在利用由常见的现代处理器支持并且通常可通过函数 fma() 访问的融合乘加运算 (FMA)，则每个半步仅需要两个 FMA。作为一个额外的好处，我们不需要特殊的舍入逻辑，因为 FMA 使用完整的产品 a*b+c 计算 a*b ，而不应用任何截断或舍入。此处以 ISO C99 版本给出的结果代码简短而有趣:

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <fenv.h>
#include <math.h>

double my_sqrt (double a)
{
    double b, r, v, w;
    float bb, rr, ss;
    int e, t, f;

    if ((a <= 0) || isinf (a) || isnan (a)) {
        if (a < 0) {
            r = 0.0 / 0.0;
        } else {
            r = a + a;
        }
    } else {
        /* 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) */
        b = ldexp (b, t);
        bb = (float)b;
        /* compute reciprocal square root */
        ss = 1.0f / bb;
        rr = sqrtf (ss);
        r = (double)rr;
        /* Use A. Schoenhage's coupled iteration for the square root */
        v = 0.5 * r;
        w = b * r;
        w = fma (fma (w, -w, b), v, w);
        v = fma (fma (r, -w, 1), v, v);
        w = fma (fma (w, -w, b), v, w);
        /* map back from primary approximation interval by jamming exponent */
        r = ldexp (w, f);
    }
    return r;
}

/* Professor George Marsaglia's 64-bit KISS PRNG */
static uint64_t xx = 1234567890987654321ULL;
static uint64_t cc = 123456123456123456ULL;
static uint64_t yy = 362436362436362436ULL;
static uint64_t zz = 1066149217761810ULL;
static uint64_t tt;
#define MWC64  (tt = (xx << 58) + cc, cc = (xx >> 6), xx += tt, cc += (xx < tt), xx)
#define XSH64  (yy ^= (yy << 13), yy ^= (yy >> 17), yy ^= (yy << 43))
#define CNG64  (zz = 6906969069ULL * zz + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

int main (void)
{
    volatile union {
        double f;
        unsigned long long int i;
    } arg, res, ref;
    unsigned long long int count = 0ULL;

    do {
        arg.i = KISS64;
        ref.f = sqrt (arg.f);
        res.f = my_sqrt (arg.f);
        if (res.i != ref.i) {
            printf ("\n!!!! arg=% 23.16e %016llx  res=% 23.16e %016llx  ref=% 23.16e %016llx\n",
                    arg.f, arg.i, res.f, res.i, ref.f, ref.i);
        }
        count++;
        if ((count & 0xffffff) == 0) printf ("\rtests = %llu", count);
    } while (1);
    return EXIT_SUCCESS;
}

在两个连续的 binades 上对这段代码进行详尽测试将需要一小群机器大约一周左右的时间，这里我包括了一个使用随机操作数的快速“冒烟”测试。

在不支持 FMA 操作的硬件上，fma() 将基于仿真。这很慢，并且已经证明有几个这样的仿真是错误的。 Schönhage 迭代在没有 FMA 的情况下也能正常工作，但在这种情况下必须添加额外的舍入逻辑。在支持截断(舍入为零)浮点乘法的情况下，最简单的解决方案是使用 Tuckerman rounding 。否则，可能需要将 double 参数和初步结果重新解释为 64 位整数，并借助整数运算执行舍入。