15

我想生成两到三百万位数的平方根的位数。

我知道Newton-Raphson，但由于缺乏大整数支持，我不太了解如何在 C 或 C++ 中实现它。有人可以指出我正确的方向吗？

另外，如果有人知道如何在 python 中做到这一点（我是初学者），我也将不胜感激。

4

9 回答 9

7

您可以尝试使用映射：

a/b -> (a+2b)/(a+b)从a= 1, b= 1. 这收敛到 sqrt(2)（实际上给出了它的连分数表示）。

现在关键点：这可以表示为矩阵乘法（类似于斐波那契）

如果 a_n 和 b_n 是步骤中的第 n 个数字，则

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

现在给了我们

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

因此，如果 2x2 矩阵是 A，我们需要计算 A n，这可以通过重复平方来完成，并且只使用整数运算（因此您不必担心精度问题）。

另请注意，您获得的 a/b 将始终采用简化形式（如 gcd(a,b) = gcd(a+2b, a+b)），因此如果您正在考虑使用分数类来表示中间结果，不要！

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

请注意，即使您正在寻找第 360 万项，重复平方也将允许您计算 O(Log n) 乘法和加法中的第 n 项。

此外，与 Newton-Raphson 等迭代式不同，这很容易实现并行化。

于 2011-03-04T16:08:06.807 回答
6

编辑：我比以前更喜欢这个版本。这是一个接受整数和小数的通用解决方案；当 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
于 2011-03-03T23:33:55.470 回答
3

至于任意大数字，您可以查看The GNU Multiple Precision Arithmetic Library (for C/C++)。

于 2011-03-03T22:59:13.893 回答
3

为了工作？使用图书馆！

为了娱乐？对你有益 ：）

编写一个程序来模仿你用铅笔和纸做的事情。从 1 位数字开始，然后是 2 位数字，然后是 3，...，...

不要担心牛顿或其他任何人。就按照你的方式去做吧。

于 2011-03-03T23:04:04.000 回答
3

最好的方法可能是使用[1; 2, 2, ...]二的平方根的连分数展开。

def root_two_cf_expansion():
    yield 1
    while True:
        yield 2

def z(a,b,c,d, contfrac):
    for x in contfrac:
        while a > 0 and b > 0 and c > 0 and d > 0:
            t = a // c
            t2 = b // d
            if not t == t2:
                break
            yield t
            a = (10 * (a - c*t))
            b = (10 * (b - d*t))
            # continue with same fraction, don't pull new x
        a, b = x*a+b, a
        c, d = x*c+d, c
    for digit in rdigits(a, c):
        yield digit

def rdigits(p, q):
    while p > 0:
        if p > q:
           d = p // q
           p = p - q * d
        else:
           d = (10 * p) // q
           p = 10 * p - q * d
        yield d

def decimal(contfrac):
    return z(1,0,0,1,contfrac)

decimal((root_two_cf_expansion())返回所有十进制数字的迭代器。 t1并且t2在算法中是下一位的最小值和最大值。当它们相等时，我们输出那个数字。

请注意，这不处理某些例外情况，例如连分数中的负数。

（此代码改编自 Haskell 代码，用于处理一直浮动的连分数。）

于 2011-03-04T06:05:25.353 回答
2

这是一个用于计算整数a到精度位数的平方根的简短版本。它通过在乘以 10 到 2 x数字后找到a的整数平方根来工作。

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 到 3 百万小数位的平方根可能需要一个小时或更长时间。

我还没有证明循环将永远终止。它可能在最后一位不同的两个值之间振荡。或者它可能不会。

于 2011-03-04T05:00:05.597 回答
2

好了，下面是我写的代码。对我来说，它在大约 60800 秒内为 2 的平方根生成了小数点后一百万位，但是我的笔记本电脑在运行程序时正在休眠，它应该更快。您可以尝试生成 300 万位数字，但可能需要几天时间才能获得。

def sqrt(number,digits_after_decimal=20):
import time
start=time.time()
original_number=number
number=str(number)
list=[]
for a in range(len(number)):
    if number[a]=='.':
        decimal_point_locaiton=a
        break
    if a==len(number)-1:
        number+='.'
        decimal_point_locaiton=a+1
if decimal_point_locaiton/2!=round(decimal_point_locaiton/2):
    number='0'+number
    decimal_point_locaiton+=1
if len(number)/2!=round(len(number)/2):
    number+='0'
number=number[:decimal_point_locaiton]+number[decimal_point_locaiton+1:]
decimal_point_ans=int((decimal_point_locaiton-2)/2)+1
for a in range(0,len(number),2):
    if number[a]!='0':
        list.append(eval(number[a:a+2]))
    else:
        try:
            list.append(eval(number[a+1]))
        except IndexError:
            pass
p=0
c=list[0]
x=0
ans=''
for a in range(len(list)):
    while c>=(20*p+x)*(x):
        x+=1
    y=(20*p+x-1)*(x-1)
    p=p*10+x-1
    ans+=str(x-1)
    c-=y
    try:
        c=c*100+list[a+1]
    except IndexError:
        c=c*100
while c!=0:
    x=0
    while c>=(20*p+x)*(x):
        x+=1
    y=(20*p+x-1)*(x-1)
    p=p*10+x-1
    ans+=str(x-1)
    c-=y
    c=c*100
    if len(ans)-decimal_point_ans>=digits_after_decimal:
            break
ans=ans[:decimal_point_ans]+'.'+ans[decimal_point_ans:]
total=time.time()-start
return ans,total
于 2011-08-09T22:24:20.393 回答
0

Python 已经支持开箱即用的大整数，如果这是在 C/C++ 中阻碍你的唯一因素，你总是可以自己编写一个快速的容器类。

您提到的唯一问题是缺少大整数。如果您不想为此使用库，那么您是否正在寻找编写此类课程的帮助？

于 2011-03-03T23:16:05.180 回答
0

这是一个更有效的整数平方根函数（在 Python 3.x 中），它应该在所有情况下都终止。它从一个更接近平方根的数字开始，因此它需要更少的步骤。请注意，int.bit_length 需要 Python 3.1+。为简洁起见省略了错误检查。

def isqrt(n):
    x = (n >> n.bit_length() // 2) + 1
    result = (x + n // x) // 2
    while abs(result - x) > 1:
        x = result
        result = (x + n // x) // 2
    while result * result > n:
        result -= 1
    return result
于 2011-03-04T11:40:36.327 回答