algorithm - 如何计算和存储 sqrt(n) 最多 10^6 位小数的数字？

Question

我正在做研究工作。为此，我需要计算并存储 2 到 10^6 位的平方根。我已经用谷歌搜索了这个，但我只有一个 NASA 页面，但我不知道他们是如何计算的。我使用了 c++ 的 set_precision。但这只会使结果最多只有大约 50 个位置。我该怎么办？

NASA 页面链接:https://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil
我也尝试过二进制搜索，但没有成功。

long double ans = sqrt(n);
cout<<fixed<<setprecision(50)<<ans<<endl;

Answer 1

您在这里有多种选择。您可以使用任意精度的浮点库(例如 MPFR 使用 C 或 C++，或 mpmath 或 Python 中的内置 decimal 库)。只要您知道库提供的错误保证，您就可以确保获得正确的十进制数字。例如，MPFR 和 Python 的 decimal保证正确的舍入，但 MPFR 的缺点(对于获取十进制数字的特定用例)它以二进制工作，因此您还需要分析由二进制到十进制转换引起的错误。

您还可以使用纯整数方法，使用任意精度整数库(如 GMP)，或开箱即用支持任意精度整数的语言(例如，Java 及其 BigInteger 类:Java 的最新版本提供 BigInteger.sqrt 方法):规模 2通过 10**2n , 其中 n是您需要的小数点后的位数，取整数平方根(即精确数学平方根的整数部分)，然后按 10**n 缩小.请参阅下面的计算整数平方根的相对简单但有效的算法。

如果您愿意使用另一种语言，这里最简单的开箱即用选项是使用 Python 的 decimal。图书馆。这是您需要的所有代码，假设是 Python 3(不是 Python 2，这会非常慢)。

>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 10**6 + 1  # number of significant digits needed
>>> sqrt2_digits = str(Decimal(2).sqrt())

str(Decimal(2).sqrt())在我的机器上操作不到 10 秒。让我们检查一下长度，以及前百位和后百位(我们显然无法在这里重现整个输出):

>>> len(sqrt2_digits)
1000002
>>> sqrt2_digits[:100]
'1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157'
>>> sqrt2_digits[-100:]
'2637136344700072631923515210207475200984587509349804012374947972946621229489938420441930169048412044'

这样做有一个小问题:保证结果是正确的四舍五入，但这是四舍五入的，而不是截断的。所以这意味着最后的“4”位可能是最后一次向上取整的结果——也就是说，该位置的实际数字可能是“3”，后面是“8”或“9”(例如)它。

我们可以通过计算几个额外的数字，然后截断它们来解决这个问题(在仔细检查这些额外数字的舍入不会影响截断之后)。

>>> getcontext().prec = 10**6 + 3
>>> sqrt2_digits = str(Decimal(2).sqrt())
>>> sqrt2_digits[-102:]
'263713634470007263192351521020747520098458750934980401237494797294662122948993842044193016904841204391'

所以确实小数点后的第 100 位是 3，而不是 4。请注意，如果上面计算的最后 3 位是“400”，我们仍然不知道第 100 位是“3”还是“4”，因为“400”可能再次是向上取整的结果。在这种情况下，您可以计算另外两位数并重试，依此类推，当您有明确的输出时停止。 (如需进一步阅读，请搜索“制表者的困境”。)

(请注意，将 decimal 模块的舍入模式设置为 ROUND_DOWN 在这里不起作用，因为 Decimal.sqrt 方法忽略了舍入模式。)

如果你想使用纯整数算术来做到这一点，Python 3.8 提供了 math.isqrt计算精确整数平方根的函数。在这种情况下，我们将按如下方式使用它:

>>> from math import isqrt
>>> sqrt2_digits = str(isqrt(2*10**(2*10**6)))

这需要更长的时间:在我的笔记本电脑上大约需要 20 秒。其中一半时间用于 str 中隐含的二进制到十进制转换。称呼。但是这一次，我们直接得到了截断的结果，不必担心四舍五入给我们最终数字错误的可能性。

再次检查结果:

>>> len(sqrt2_digits)
1000001
>>> sqrt2_digits[:100]
'1414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572'
>>> sqrt2_digits[-100:]
'2637136344700072631923515210207475200984587509349804012374947972946621229489938420441930169048412043'

这有点作弊，因为(在撰写本文时)Python 3.8 尚未发布，尽管 beta 版本 are available .但是有一个纯 Python 版本的 isqrt CPython source 中的算法, 可以直接复制粘贴使用。这是完整的:

import operator

def isqrt(n):
    """
    Return the integer part of the square root of the input.
    """
    n = operator.index(n)

    if n < 0:
        raise ValueError("isqrt() argument must be nonnegative")
    if n == 0:
        return 0

    c = (n.bit_length() - 1) // 2
    a = 1
    d = 0
    for s in reversed(range(c.bit_length())):
        # Loop invariant: (a-1)**2 < (n >> 2*(c - d)) < (a+1)**2
        e = d
        d = c >> s
        a = (a << d - e - 1) + (n >> 2*c - e - d + 1) // a

    return a - (a*a > n)

该来源还包含对上述算法的解释及其正确性的非正式证明。

您可以检查上述两种方法的结果是否一致(以第一个结果中的额外小数点为模)。它们是通过完全不同的方法计算的，因此可以对这两种方法进行完整性检查。