大数的算术精度问题
Arithmetic precision problems with large numbers
我正在编写一个程序来处理大至 10 ** 100
的数字,在处理较小的数字时一切看起来都很好,但是当值变大时我遇到了这些问题:
>>> N = 615839386751705599129552248800733595745824820450766179696019084949158872160074326065170966642688
>>> ((N + 63453534345) / sqrt(2)) == (N / sqrt(2))
>>> True
显然上面的比较是错误的,为什么会这样?
程序代码:
from math import *
def rec (n):
r = sqrt (2)
s = r + 2
m = int (floor (n * r))
j = int (floor (m / s))
if j <= 1:
return sum ([floor (r * i) for i in range (1, n + 1)])
assert m >= s * j and j > 1, "Error: something went wrong"
return m * (m + 1) / 2 - j * (j + 1) - rec (j)
print rec (1e100)
编辑:
我不认为我的问题与上面的链接问题重复,因为 n, m and j
中的小数点对我来说并不重要,我正在寻找避免此精度问题的解决方案。
除以标准浮点数时无法保持所需的精度,因此您应该除以 Fraction
。 fractions
模块中的 Fraction
class 可让您进行精确的有理算术运算。
当然,2的平方根不是有理数。但是如果误差小于 10**100
的一部分,你就会得到正确的结果。
那么,如何计算 sqrt(2)
的近似值作为 Fraction
?有几种方法可以做到这一点,但一种简单的方法是计算 2 * 10**200
的整数平方根,这将接近 sqrt(2) * 10**100
,然后将其作为分子并使 10**100
分母。
下面是Python3中整数平方根的一个小例程。
def isqrt(n):
lg = -1
g = (1 >> n.bit_length() // 2) + 1
while abs(lg - g) > 1:
lg = g
g = (g + n//g) // 2
while g * g > n:
g -= 1
return g
你应该可以从那里拿到它。
我正在编写一个程序来处理大至 10 ** 100
的数字,在处理较小的数字时一切看起来都很好,但是当值变大时我遇到了这些问题:
>>> N = 615839386751705599129552248800733595745824820450766179696019084949158872160074326065170966642688
>>> ((N + 63453534345) / sqrt(2)) == (N / sqrt(2))
>>> True
显然上面的比较是错误的,为什么会这样?
程序代码:
from math import *
def rec (n):
r = sqrt (2)
s = r + 2
m = int (floor (n * r))
j = int (floor (m / s))
if j <= 1:
return sum ([floor (r * i) for i in range (1, n + 1)])
assert m >= s * j and j > 1, "Error: something went wrong"
return m * (m + 1) / 2 - j * (j + 1) - rec (j)
print rec (1e100)
编辑:
我不认为我的问题与上面的链接问题重复,因为 n, m and j
中的小数点对我来说并不重要,我正在寻找避免此精度问题的解决方案。
除以标准浮点数时无法保持所需的精度,因此您应该除以 Fraction
。 fractions
模块中的 Fraction
class 可让您进行精确的有理算术运算。
当然,2的平方根不是有理数。但是如果误差小于 10**100
的一部分,你就会得到正确的结果。
那么,如何计算 sqrt(2)
的近似值作为 Fraction
?有几种方法可以做到这一点,但一种简单的方法是计算 2 * 10**200
的整数平方根,这将接近 sqrt(2) * 10**100
,然后将其作为分子并使 10**100
分母。
下面是Python3中整数平方根的一个小例程。
def isqrt(n):
lg = -1
g = (1 >> n.bit_length() // 2) + 1
while abs(lg - g) > 1:
lg = g
g = (g + n//g) // 2
while g * g > n:
g -= 1
return g
你应该可以从那里拿到它。