生成平方根2的数字

Question

我想生成数字2的平方根，达到300万位数字。

我知道有一种叫做牛顿-拉夫森法的方法，但由于缺少大整数支持，我对如何在C或C++中实现它不是很清楚。有没有人能给我指个方向？

另外，如果有人知道如何用Python实现（我还是个初学者），我也会很感激。

Answer 1

这里有一个简短的方法，可以计算一个整数a的平方根，精确到digits位小数。这个方法的原理是先把a乘以10的2倍digits次方，然后再找出它的整数平方根。

def sqroot(a, digits):
    a = a * (10**(2*digits))
    x_prev = 0
    x_next = 1 * (10**digits)
    while x_prev != x_next:
        x_prev = x_next
        x_next = (x_prev + (a // x_prev)) >> 1
    return x_next

不过，有几点需要注意。

你需要把结果转换成字符串，并在正确的位置加上小数点（如果你想显示小数点的话）。

把一个非常大的整数转换成字符串的速度并不是很快。

在Python中，处理非常大的整数进行除法运算的速度也不快。

根据你电脑的性能，计算2到300万位小数的平方根可能需要一个小时或更长时间。

我还没有证明这个循环一定会结束。它可能会在两个值之间来回摆动，这两个值的最后一位数字不同。也可能不会。

Answer 2

编辑: 我觉得这个版本比之前的更好。它是一个通用的解决方案，可以同时处理整数和小数；当 n = 2 和精度 = 100000 时，大约需要两分钟。感谢 Paul McGuire 的建议，其他建议也欢迎！

def sqrt_list(n, precision):
    ndigits = []        # break n into list of digits
    n_int = int(n)
    n_fraction = n - n_int

    while n_int:                            # generate list of digits of integral part
        ndigits.append(n_int % 10)
        n_int /= 10
    if len(ndigits) % 2: ndigits.append(0)  # ndigits will be processed in groups of 2

    decimal_point_index = len(ndigits) / 2  # remember decimal point position
    while n_fraction:                       # insert digits from fractional part
        n_fraction *= 10
        ndigits.insert(0, int(n_fraction))
        n_fraction -= int(n_fraction)
    if len(ndigits) % 2: ndigits.insert(0, 0)  # ndigits will be processed in groups of 2

    rootlist = []
    root = carry = 0                        # the algorithm
    while root == 0 or (len(rootlist) < precision and (ndigits or carry != 0)):
        carry = carry * 100
        if ndigits: carry += ndigits.pop() * 10 + ndigits.pop()
        x = 9
        while (20 * root + x) * x > carry:
                x -= 1
        carry -= (20 * root + x) * x
        root = root * 10 + x
        rootlist.append(x)
    return rootlist, decimal_point_index

Answer 3

你可以尝试使用这个公式：

a/b -> (a+2b)/(a+b)，从 a= 1, b= 1 开始。这会收敛到平方根2（实际上，它给出了平方根2的连分数表示）。

现在关键点是：这可以用矩阵乘法来表示（类似于斐波那契数列）。

如果 a_n 和 b_n 是第n步的数字，那么

[1 2] [a_n b_n]^T = [a_(n+1) b_(n+1)]^T
[1 1]

这样我们就得到了

[1 2]ⁿ [a_1 b_1]^T = [a_(n+1) b_(n+1)]^T
[1 1]

因此，如果这个2x2的矩阵是A，我们需要计算 An，这可以通过重复平方的方法来完成，而且只用整数运算（所以你不需要担心精度问题）。

还要注意，你得到的 a/b 始终是最简形式（因为 gcd(a,b) = gcd(a+2b, a+b)），所以如果你考虑使用分数类来表示中间结果，那就别这样做！

由于第n个分母像是 (1+sqrt(2))^n，要得到300万位数字，你可能需要计算到第3671656项。

注意，尽管你在寻找大约360万的项，重复平方的方法可以让你在 O(Log n) 的乘法和加法中计算出第n项。

另外，这个方法很容易并行处理，不像牛顿-拉夫森等迭代方法。