计算 (1 - SQRT(0.5)) 的定点表示形式至任意精度级别

Question

我有一个常数

0.29289321881345247559915563789515...

，可以使用方程

(1 - SQRT(0.5))

计算它，然后转换为定点格式，以便与各种大小的无符号整数一起使用。示例请参见下表：

字号	Q格式	十进制值	二进制值
8 位	Q0.3	5	101
16 位	Q0.7	75	1001011
32 位	Q0.15	19195	100101011111011
64 位	Q0.31	1257966796	1001010111110110000110011001100
128 位	Q0.63	5402926248376769403	100101011111011000011001100110000000110001000011001101101111011

问题是计算这些值变得不可行，因为它涉及平方根函数和可能变得任意大的数字。我记得学习过一种“龙头算法”，该算法使用有限的内存量输出 π 的数字，并希望对于如上所述的常量存在类似的东西。 这篇论文在我的搜索过程中出现，但我还没有能够很好地理解它以翻译成代码。

如何使用类似 C 的语言（最好是 C#）来实现这一目标？有没有更好的方法来实现计算 2 的幂字长值的目标？

---

额外背景

我有以下来自另一个问题的 C# 片段：

y -= uint.CreateChecked(value: BinaryIntegerConstants<T>.Size) switch {
    8U => (x * ((y * T.CreateChecked(value: 5UL)) >> 4)),
    16U => (x * ((y * T.CreateChecked(value: 75UL)) >> 8)),
    32U => (x * ((y * T.CreateChecked(value: 19195UL)) >> 16)),
    64U => (x * ((y * T.CreateChecked(value: 1257966796UL)) >> 32)),
    128U => (x * ((y * T.CreateChecked(value: 5402926248376769403UL)) >> 64)),
    _ => throw new NotSupportedException(), // TODO: Research a way to calculate the proper constant at runtime.
}

我已经使用C#代码手动计算了上面的常量，但由于

double

类型的限制，无法进一步计算。

Answer 1

我不知道二进制或十六进制 sqrt(1/2) 的插口算法，但这并不意味着不存在。这不是我的领域，但我对 PLSQ 感兴趣，它可以处理类似的近似问题。

您记得的 pi 的

BPP 算法有一定的特殊性。因为可以使用由复杂的计算机辅助高精度实验数学获得的非常具体的公式来正确计算 pi 的十六进制表示形式的孤立数字。当它在 1995 年被发现时，我们感到非常惊讶。后来有关 PLSQ 和 LLL 的论文仍然使用它的推导作为更新更快的整数关系和格约简算法的威力的例子。

但是，可以使用 Newton-Raphson 选项，每次使用它时，有效位数都会加倍，并且如果您小心的话，可以以缩放整数任意精度形式工作。

想要一个精确的 sqrt(1/2) 在这里特别有帮助。在基数 2 中，乘以 2 和除以 2 是微不足道的。

我认为您可能有一条前进的道路，使用旧的基于牛顿拉夫森的方法的缩放整数版本来计算 1/sqrt(y)，其中 y 在本例中为 2。倒数平方根算法的工作原理如下和日期在此之前，有 FP 硬件支持可以快速完成此操作。一些 cbrt 仍然使用这种方法，巧妙地避免了缓慢的除法。与所有其他二元运算相比，今天的除法仍然相当慢。下面是rsqrt()的简单推导

f(x)  = y - 1/x^2
f'(x) = 2/x^3

令起始猜测为 x0，然后 NR 迭代给出

x1 = x0 - f(x0)/f'(x0)

x1 = x0 - (y - 1/x0^2)*x0^3/2  = x0 - x0*(y*x0^2-1)/2

当 NR 公式进一步简化时，设置 y = 2（您的特殊情况）

x1 = x0 - x0*(x0*x0-1/2)  (#1)

x1 = 1 - 1*(1 - 1/2)    = 1/2
x2 = 1/2 - 1/2*(1/4-1/2) = 5/8
x3 = 3/4 - 5/8*(25/64-1/2) = 355/512

这样您就可以将每次迭代的准确位数加倍。 计算成本相对适中 2 * 和 2 - 以先前结果长度的两倍计算。

但是要小心任何固定深度的最后几位，因为尽管算法是自我校正的，但下一阶段的进位可能会影响当前最低有效部分中的几个位。

我知道（#1）可以重写为，但这种形式不太准确（至少在正常的 FP 工作中）

x1 = x0*(3/2-x0*x0) (#1a) 您也可以使用哈雷方法来完成此操作，但这需要在完全工作精度下进行非平凡的除法，因此这里肯定会更慢。

我倾向于在允许轻松访问快速任意精度浮点（如 Julia 或 Python）的环境中预先计算这些值，并根据需要在 C# 程序中运行时导入它们