这段代码中查找 table 的目的是什么？
What purpose is the lookup table serving in this code?
下面我改编了 William Kahan and K.C. Ng (look at the comment block on the bottom) written in 1986 to produce an approximation of 1 / sqrt(x) where x is an IEEE-754 double precision floating point number. It is this code that Cleve Moler and Greg Walsh adapted to become the "fast inverse square root 中因 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 的近似二进制对数 (Mitchell's method)。
我们将这些位右移并从常数中减去它们（对于此比较，低 32 位已被删除并不重要），这是对数域中的近似除法。
下一步我不太清楚，因为查找 table 是由非常少量的原始数字索引的——主要是指数位。所以发生了修正，但我不确定如何解释它。
最后我们将整数位（加上 LS 半部分的 32 个 0）转换为浮点双精度数，得到近似的反对数。
所以我理解（或认为我理解）这里 4 个步骤中的 3 个，但第三步 - 查找 table 正在做什么，如何在双精度范围内建立索引，以及为什么？
xi = *(unsigned long long int*)&x;
这会将表示 double x 的位重新解释为 unsigned long long，这应该是 64 位类型。这种重新解释的方法未由 C 标准定义，但受到许多编译器的支持，可能受命令行开关的影响。
重新解释为 unsigned long long 使这些位易于操作。
/* Extract only the first half of the bits of the integer */
x0 = (xi & 0xffffffff00000000) >> 32;
正如它所说，x0 现在包含 double 的高 32 位。这包括符号位（大概是 0，因为我们应该只取正数的平方根倒数）、11 位编码指数和 20 位主要编码的尾数。
/* Perform bit shift and subtract from an integer constant to get k */
k = 0x5fe80000 - (x0 >> 1);
x0 >> 1 开始计算平方根的近似值。让我们看看为什么：二进制浮点表示是±f•2e，对于一些=77=]f在一定的定点格式和在一定范围内的整数e。对于正数，我们当然有 f•2e 没有符号。令q和r为e除以2的商和余数，所以e = 2q+r。然后 sqrt(f•2e) = sqrt(f•2r)•2q。在 x0 >> 1 中，我们将指数位右移，设置我们在指数字段中具有 q。
但不完全是，因为原始指数字段是有偏差的。它实际上包含e+1023。移位给我们 e/2+511½，该分数移出场。但我们会处理这个。特别是，偏置需要重置为 1023，但它内置于 5FE8000016 常量中。
接下来，- (x0 >> 1) 将其更改为平方根倒数的近似值，仅仅是因为对指数取反取倒数。不过，有效数字仍然是个问题。
那么0x5fe80000 - (x0 >> 1)就是这个非常粗略的平方根倒数近似值加上5FE8000016的double编码。
回想一下，指数编码需要与指数有 1023 的偏差，但我们将它减半并取反，所以 - (x0 >> 1) 有 −511½ 的偏差，需要调整为 +1023 , 所以我们需要加上 1534½。这需要进入指数字段，它位于 64 位 double 编码中的第 62 位到第 52 位，因此在这个 32 位部分中从高半部分开始的第 30 到 20 位。所以我们需要从第 20 位开始添加 1534½，即 1534½•220 = 1,609,039,872 = 5FE8000016.
哒哒！ 0x5fe80000 - (x0 >> 1)的指数位恰好是double2−q的指数位。它的有效位是移出的指数位加上原来的有效位右移一位（移出并丢失一位）。 （请注意，移出的指数位已被反转——它添加了 816 中的那个位。）
现在我们有了正确的指数，但是第 19 位到第 0 位中的有效位有点古怪。它们包含来自指数的一位和来自有效数字编码的一些取反和移位的位。要获得平方根倒数的近似值，我们必须解决这些问题。我们可以采用移位的位并计算出原始有效数是多少（达到我们从中得到的 19 位给出的精度）并乘以 2r，这将为我们提供到目前为止我们准备的指数中未考虑的原始数字部分。然后我们取那部分的平方根倒数，用它代替有效数，就完成了。
这是一个相当大的计算量。相反，我们将使用其他一些软件预先计算结果并将结果记录在 table.
中
63 & (k >> 14) 提取第 19 到 14 位，将它们留在第 5 到 0 位。这些是从指数移出的位和前五个有效位。因此，它们是原始数字中尚未在我们计算的指数字段中考虑的最高有效位。
T2[63 & (k >> 14)] 使用这些位在 table 中查找值。该值是用其他软件预先计算的。它包含从 k 中的位到我们想要的位的调整。然后 k - T2[63 & (k >> 14)] 应用该调整，因此它是平方根倒数的近似值（特别是该近似值的 double 编码的高 32 位）。
您引用的 material 没有说明 table 是如何计算的。由于每个 table 条目仅由 6 位索引，因此它将用于 x 的许多 double 值。该条目可能是通过计算该条目将用于的 x 的最小值的平方根倒数，该条目将用于的最大值 x 的平方根倒数来计算的，并在这两者之间取一些值。然后设置 table 条目，以便在从 k 中减去它时，它会取消用于索引 table 的六位，然后添加对目标值进行编码的位。 （就像指数调整添加 511½ 和 1023 一样，我们希望 table 条目减去不需要的 63 & (k >> 14) 并添加所需的位。）
这涉及到一些技巧，因为 table 条目可以取消用于索引 table 的六位（因为，对于每个 table 条目，我们确切地知道哪六位用于查找它），但它不能取消较低的位（因为它们不是索引的一部分并且可以变化）。因此需要一些额外的工作来 select 使用哪个特定的目标值。此外，使用哪个目标值受我们是否要最小化绝对误差或相对误差或其他因素的影响。你在问题中展示的内容并没有说明这一点，所以我不能准确地说出 table 是如何计算的。
请注意，代码不完整，如果没有进一步的工作或预防措施，不应将其用作平方根倒数。我预计它不能正确处理次正规、无穷大或 NaN。
补充
Norbert Juffa 提供此替代方案 table 以最小化上述代码产生的近似值中的相对误差（这可能只是用于进一步细化步骤的初始近似值）：
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
};
下面我改编了 William Kahan and K.C. Ng (look at the comment block on the bottom) written in 1986 to produce an approximation of 1 / sqrt(x) where x is an IEEE-754 double precision floating point number. It is this code that Cleve Moler and Greg Walsh adapted to become the "fast inverse square root 中因 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 的近似二进制对数 (Mitchell's method)。
我们将这些位右移并从常数中减去它们（对于此比较，低 32 位已被删除并不重要），这是对数域中的近似除法。
下一步我不太清楚，因为查找 table 是由非常少量的原始数字索引的——主要是指数位。所以发生了修正，但我不确定如何解释它。
最后我们将整数位（加上 LS 半部分的 32 个 0）转换为浮点双精度数，得到近似的反对数。
所以我理解（或认为我理解）这里 4 个步骤中的 3 个，但第三步 - 查找 table 正在做什么，如何在双精度范围内建立索引，以及为什么？
xi = *(unsigned long long int*)&x;
这会将表示 double x 的位重新解释为 unsigned long long，这应该是 64 位类型。这种重新解释的方法未由 C 标准定义，但受到许多编译器的支持，可能受命令行开关的影响。
重新解释为 unsigned long long 使这些位易于操作。
/* Extract only the first half of the bits of the integer */
x0 = (xi & 0xffffffff00000000) >> 32;
正如它所说，x0 现在包含 double 的高 32 位。这包括符号位（大概是 0，因为我们应该只取正数的平方根倒数）、11 位编码指数和 20 位主要编码的尾数。
/* Perform bit shift and subtract from an integer constant to get k */
k = 0x5fe80000 - (x0 >> 1);
x0 >> 1 开始计算平方根的近似值。让我们看看为什么：二进制浮点表示是±f•2e，对于一些=77=]f在一定的定点格式和在一定范围内的整数e。对于正数，我们当然有 f•2e 没有符号。令q和r为e除以2的商和余数，所以e = 2q+r。然后 sqrt(f•2e) = sqrt(f•2r)•2q。在 x0 >> 1 中，我们将指数位右移，设置我们在指数字段中具有 q。
但不完全是，因为原始指数字段是有偏差的。它实际上包含e+1023。移位给我们 e/2+511½，该分数移出场。但我们会处理这个。特别是，偏置需要重置为 1023，但它内置于 5FE8000016 常量中。
接下来，- (x0 >> 1) 将其更改为平方根倒数的近似值，仅仅是因为对指数取反取倒数。不过，有效数字仍然是个问题。
那么0x5fe80000 - (x0 >> 1)就是这个非常粗略的平方根倒数近似值加上5FE8000016的double编码。
回想一下，指数编码需要与指数有 1023 的偏差，但我们将它减半并取反，所以 - (x0 >> 1) 有 −511½ 的偏差，需要调整为 +1023 , 所以我们需要加上 1534½。这需要进入指数字段，它位于 64 位 double 编码中的第 62 位到第 52 位，因此在这个 32 位部分中从高半部分开始的第 30 到 20 位。所以我们需要从第 20 位开始添加 1534½，即 1534½•220 = 1,609,039,872 = 5FE8000016.
哒哒！ 0x5fe80000 - (x0 >> 1)的指数位恰好是double2−q的指数位。它的有效位是移出的指数位加上原来的有效位右移一位（移出并丢失一位）。 （请注意，移出的指数位已被反转——它添加了 816 中的那个位。）
现在我们有了正确的指数，但是第 19 位到第 0 位中的有效位有点古怪。它们包含来自指数的一位和来自有效数字编码的一些取反和移位的位。要获得平方根倒数的近似值，我们必须解决这些问题。我们可以采用移位的位并计算出原始有效数是多少（达到我们从中得到的 19 位给出的精度）并乘以 2r，这将为我们提供到目前为止我们准备的指数中未考虑的原始数字部分。然后我们取那部分的平方根倒数，用它代替有效数，就完成了。
这是一个相当大的计算量。相反，我们将使用其他一些软件预先计算结果并将结果记录在 table.
中63 & (k >> 14) 提取第 19 到 14 位，将它们留在第 5 到 0 位。这些是从指数移出的位和前五个有效位。因此，它们是原始数字中尚未在我们计算的指数字段中考虑的最高有效位。
T2[63 & (k >> 14)] 使用这些位在 table 中查找值。该值是用其他软件预先计算的。它包含从 k 中的位到我们想要的位的调整。然后 k - T2[63 & (k >> 14)] 应用该调整，因此它是平方根倒数的近似值（特别是该近似值的 double 编码的高 32 位）。
您引用的 material 没有说明 table 是如何计算的。由于每个 table 条目仅由 6 位索引，因此它将用于 x 的许多 double 值。该条目可能是通过计算该条目将用于的 x 的最小值的平方根倒数，该条目将用于的最大值 x 的平方根倒数来计算的，并在这两者之间取一些值。然后设置 table 条目，以便在从 k 中减去它时，它会取消用于索引 table 的六位，然后添加对目标值进行编码的位。 （就像指数调整添加 511½ 和 1023 一样，我们希望 table 条目减去不需要的 63 & (k >> 14) 并添加所需的位。）
这涉及到一些技巧，因为 table 条目可以取消用于索引 table 的六位（因为，对于每个 table 条目，我们确切地知道哪六位用于查找它），但它不能取消较低的位（因为它们不是索引的一部分并且可以变化）。因此需要一些额外的工作来 select 使用哪个特定的目标值。此外，使用哪个目标值受我们是否要最小化绝对误差或相对误差或其他因素的影响。你在问题中展示的内容并没有说明这一点，所以我不能准确地说出 table 是如何计算的。
请注意，代码不完整，如果没有进一步的工作或预防措施，不应将其用作平方根倒数。我预计它不能正确处理次正规、无穷大或 NaN。
补充
Norbert Juffa 提供此替代方案 table 以最小化上述代码产生的近似值中的相对误差（这可能只是用于进一步细化步骤的初始近似值）：
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
};