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

Question

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

我的问题是：如果我想计算 64 位平方根，我可以使用 32 位平方根对我有帮助吗？拉夫森之类的？

Answer 1

可能没有

假设我们有一个 64 位整型变量 i64 的两部分作为 hi 和 lo 那么

sqrt(i64) = sqrt(hi*2³² + lo)

我们没有办法将和的平方根简化为另一个表达式，因此我们无法从 32 位平方根计算 64 位平方根

但是你说你有一个 64 位浮点值。你在没有 FPU 的平台上吗？您的 32 位平方根是浮点函数还是整数函数？无论如何都会出现同样的问题，因为尾数不适合单个寄存器，但如果不需要全精度，您可以获得一些近似值

fast square root optimization?
https://en.wikipedia.org/wiki/Fast_inverse_square_root

Answer 2

您仍然需要对 Newton-Raphson 进行编程，但是您可以通过使用 32 位平方根计算出 32 位近似值并将其用作 Newton-Raphson 的起始值来节省大量迭代，这意味着它将在更少的迭代中收敛到完全正确的答案。这是值得节省的时间 - 硬件平方根有时使用 table 查找来查找 Newton-Raphson 的起点，并且最佳的理论复杂度计算假设您对较早的迭代使用较低的精度以节省时间。

Answer 3

TL;DR 是。

根据您平台的硬件、工具链和数学库的功能和不足，这可能不一定是计算双精度平方根的最快或最不痛苦的方法。下面我展示了一种基于 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_i)
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开始，只需要两步就可以达到双精度精度

如果我们现在利用融合乘加运算 (FMA)，该运算由常见的现代处理器支持并且通常可通过函数访问 fma()，每个半步仅需要两个 FMA。作为额外的好处，我们不需要特殊的舍入逻辑，因为 FMA 使用完整乘积 a*b 计算 a*b+c，而不应用任何截断或舍入。此处以 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;
}

在两个连续的二进制文件中对这段代码进行详尽测试将需要一小群机器大约一周左右的时间，这里我包括了一个使用随机操作数的快速 "smoke" 测试。

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

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

Can a function of a 32-bit square root be used to help calculate a 64-bit square root?

algorithm

math

sqrt