0

下面我改编了William Kahan 和 KC Ng于 1986 年编写的代码（查看底部的注释块），以产生 1 / sqrt(x) 的近似值，其中 x 是 IEEE-754 双精度浮点数。正是这个密码，Cleve Moler 和 Greg Walsh 改编成了以 Quake III 闻名的“快速反平方根”。

#include <stdio.h>
int main() {
    double x = 3.1415926535;
    double y, z;
    unsigned long int x0, y0, k;
    unsigned long long int xi, yi;
    static int T2[64]= {
    0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
    0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
    0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
    0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
    0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
    0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
    0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
    0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd };
    /* Convert double to unsigned long long (64 bits wide)  */
    xi = *(unsigned long long int*)&x;
    /* Extract only the first half of the bits of the integer */
    x0 = (xi & 0xffffffff00000000) >> 32;
    /* Perform bit shift and subtract from an integer constant to get k */
    k = 0x5fe80000 - (x0 >> 1);
    /* T2 is indexed by mostly by exponent bits. 
       Only 5 highest bits from orig 52 in mantissa are used to index T2 */
    y0 = k - T2[63 & (k >> 14)];
    /* Pad with zeros for the LS 32 bits, convert back to long long */
    yi = ((unsigned long long int) y0 << 32);
    /* Get double from bits making up unsigned long long */ 
    y = *(double*)&yi;
    /* 1/sqrt(pi) ~ 0.564 */
    printf("%lf\n", y);
    return 0;
}

看起来它正在做与 Quake 版本相同的事情（减去 newton-raphson 步骤）

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

我们首先将 double 的位转换为整数，这是我们开始使用的 double 的近似二进制对数（米切尔方法）。

我们将这些位右移并从一个常数中减去它们（对于这个比较，丢弃低 32 位并不重要），这是对数域中的近似除法。

下一步我不太清楚，因为查找表由非常少量的原始数字索引 - 主要是指数位。所以发生了更正，但我不确定如何解释它。

最后，我们将整数位（加上 LS 一半的 32 个 0）转换为浮点双精度，给我们一个近似的反对数。

所以我理解（或认为我理解）这里 4 个步骤中的 3 个，但是第三步——查找表在做什么，它是如何在双打范围内索引的，为什么？

4

1 回答 1

3 
xi = *(unsigned long long int*)&x;

这将表示 的位重新解释double x为unsigned long long，它应该是 64 位类型。这种重新解释的方法不是由 C 标准定义的，但被许多编译器支持，可能受命令行开关的影响。

重新解释为 anunsigned long long使这些位易于操作。

/* Extract only the first half of the bits of the integer */
x0 = (xi & 0xffffffff00000000) >> 32;

正如它所说，x0现在包含double. 这包括符号位（大概是 0，因为我们应该只取正数的平方根倒数）、编码指数的 11 位和有效数的主要编码的 20 位。

/* Perform bit shift and subtract from an integer constant to get k */
k = 0x5fe80000 - (x0 >> 1);

x0 >> 1开始一个平方根的近似值。让我们看看为什么：二进制浮点表示是 ± f •2 e，对于某个定点格式的某些f和某个范围内的整数e。对于正数，我们当然有不带符号的f •2 e 。设q和r是e除以 2的商和余数，所以e = 2 q + r。然后 sqrt( f •2 e ) = sqrt( f •2 r )•2 q。在x0 >> 1，我们将指数位右移，设置我们在指数字段中有q 。

不过，不完全是，因为原始指数场是有偏差的。它实际上包含e +1023。移位给我们e /2+511½，该分数移出该字段。但我们会处理这个问题。特别是，需要将偏差重置为 1023，但这是内置于 5FE80000 16常数中的。

接下来，- (x0 >> 1)将其更改为平方根倒数的近似值，仅仅是因为取反指数取倒数。但是，有效数字仍然是一个问题。

然后0x5fe80000 - (x0 >> 1)是添加了 5FE80000 16double的这个非常粗略的反平方根近似值的编码。

回想一下，指数编码需要与指数有 1023 的偏差，但我们将其减半并取反，因此- (x0 >> 1)需要将 -511½ 的偏差调整为 +1023，因此我们需要添加 1534½。这需要进入指数字段，它位于 64 位double编码中的 62 到 52 位，因此在这个 32 位部分的高半部分中的 30 到 20 位。所以我们需要从第 20 位开始添加 1534½，或者 1534½•2 20 = 1,609,039,872 = 5FE80000 16。

达达！的指数位0x5fe80000 - (x0 >> 1)正是double2 −<em>q的指数位。它的有效位是移出的指数位加上原始的有效位右移一位（其中一位移出并丢失）。（除了注意移出的指数位已被反转 - 它添加了 8 16中的那个位。）

现在我们有了正确的指数，但是位 19 到 0 中的有效位是古怪的。它们包含来自指数的一位和来自有效数字编码的一些取反和移位的位。为了得到逆平方根的近似值，我们必须解决这些问题。我们可以取出移位的位并计算出原始有效数是什么（直到我们从中获得的 19 位给出的精度）并乘以 2 r，这将给我们原始数字中未被计入的部分因为在我们迄今准备好的指数中。然后我们取那部分的平方根的倒数并把它放在有效数字的位置，我们就完成了。

这是一个相当大的计算量。相反，我们将使用其他一些软件预先计算结果并将结果记录在表格中。

63 & (k >> 14)提取位 19 到 14，将它们留在位 5 到 0 中。这些是从指数移出的位和前五个有效位。所以它们是我们计算的指数字段中尚未考虑的原始数字的最高有效位。

T2[63 & (k >> 14)]使用这些位在表中查找值。该值是使用其他软件预先计算的。它包含从其中的位k到我们想要的位的调整。然后k - T2[63 & (k >> 14)]应用该调整，因此它是平方根的近似值（特别是double该近似值编码的高 32 位）。

您引用的材料没有说明该表是如何计算的。double由于每个表条目仅由 6 位索引，因此它将用于x. 该条目可能是通过计算将x用于该条目的最小值的平方根的平方根，将用于该条目的最大值的平方根的平方根x，并在这两者之间取一些值来计算的。然后设置表条目，以便当从 中减去它时k，它会取消用于索引表的六位，然后添加对目标值进行编码的位。（就像添加了 511½ 和 1023 的指数调整一样，我们希望表格条目减去不需要63 & (k >> 14)的位并添加所需的位。）

这涉及到一些技巧，因为表条目可以取消用于索引表的六位（因为对于每个表条目，我们确切地知道用于查找它的六位），但它不能取消较低的位（因为它们不是索引的一部分并且可以变化）。因此需要一些额外的工作来选择要使用的特定目标值。此外，使用哪个目标值受我们是否要最小化绝对误差或相对误差或其他因素的影响。您在问题中显示的内容与此无关，因此我无法准确说明该表的计算方式。

请注意，代码是不完整的，如果没有进一步的工作或预防措施，不应将其用作平方根的倒数。我希望它不能正确处理次正规、无穷大或 NaN。

补充

Norbert Juffa提供了这个替代表，以最大限度地减少上述代码产生的近似值中的相对误差（这可能只是用于进一步细化步骤的初始近似值）：

uint32_t T2_NJ[64] =
    {
        0x01289, 0x02efb, 0x04d6a, 0x06b05, 0x087c3, 0x0a39b, 0x0be82, 0x0d86e,
        0x0f153, 0x10927, 0x11fdb, 0x13563, 0x149af, 0x15cb1, 0x16e57, 0x17e91,
        0x18d4a, 0x19a6e, 0x1a5e7, 0x1af9d, 0x1b777, 0x1bd59, 0x1c123, 0x1c2b7,
        0x1c1ef, 0x1bea5, 0x1b8ae, 0x1afdc, 0x1a3fc, 0x194d5, 0x18229, 0x16bb4,
        0x16874, 0x17a1c, 0x18aa4, 0x19a00, 0x1a824, 0x1b501, 0x1c08b, 0x1cab1,
        0x1d365, 0x1da95, 0x1e02e, 0x1e41f, 0x1e652, 0x1e6b1, 0x1e525, 0x1e194,
        0x1dbe4, 0x1d3f8, 0x1c9b0, 0x1bcea, 0x1ad82, 0x19b51, 0x1862c, 0x16de4,
        0x15248, 0x1331f, 0x1102e, 0x0e933, 0x0bde5, 0x08df6, 0x0590c, 0x01ec8
    };
于 2022-01-17T00:11:11.587 回答