大量的numpy精度

如何解决大量的numpy精度

我想使用费马的分解方法分解一个大数。我是这样实现的：

import numpy as np

def fac(n):
    x = np.ceil(np.sqrt(n))
    y = x*x - n
    while not np.sqrt(y).is_integer():
        x += 1
        y = x*x - n
    return(x + np.sqrt(y),x - np.sqrt(y))

使用这种方法，我想将 N 分解到其组件中。请注意 N=p*q，其中 p 和 q 是质数。

我选择了以下值来计算 N：

p = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533846141.0
q = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533845933.0

并定义N

N = p*q

现在我考虑 N：

r = fac(n)

然而，分解似乎不正确：

int(r[0])*int(r[1]) == N

它确实适用于较小的整数：

fac(65537)
Out[1]: (65537.0,1.0)

我很确定原因是某些时候的数值精度。

我尝试使用对象类型在 numpy 中计算 N：

N = np.dot(np.array(p).astype(object),np.array(q).astype(object))

但它没有帮助。尽管如此，numpy 仍需要 sqrt 函数的浮点数。

我还尝试使用 math 库而不是 numpy，该库似乎不需要为其 sqrt 函数提供浮点数，但最终也遇到了精度问题。

解决方法

Python int 是多精度数。但是 numpy 是 C 低级库的包装器，用于加速操作。缺点是它无法处理那些多精度数字。更糟糕的是，如果您尝试对它们使用 np.sqrt，它们将被转换为浮点数（C double 或 numpy float64），其精度约为 15 位十进制数字。

但由于 Python int 类型已经是多精度类型，您可以使用 math.sqrt 获得真实平方根的近似值，然后使用 Newton 找到更接近的值：

def isqrt(n):
    x = int(math.sqrt(n))
    old = None
    while True:
        d = (n - x * x) // (2 * x)
        if d == 0: break
        if d == 1:                 # infinite loop prevention
            if old is None:
                old = 1
            else: break
        x += d
    return x

使用它，您的 fac 函数可以变成：

def fac(n):
    x = isqrt(n)
    if x*x < n: x += 1
    y = x * x - n
    while True:
        z = isqrt(y)
        if z*z == y: break
        x += 1
        y = x*x -n
    return x+z,x-z

演示：

p = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533846141
q = 34058934059834598495823984675767545695711020949846845989934523432842834738974239847294083409583495898523872347284789757987987387543533845933
N = p*q
print(fac(N) == (p,q))

按预期打印True