Python中大数计算精度问题的解决方法（欧拉计划100题场景）

问题重现 #

处理超大整数时，用math.sqrt()判断是否为完全平方数会出现精度误判：

示例代码：

import math
i = 292893227695
j = 8*i**2 + 1
k = math.sqrt(j)
print(j)
print(k)
print(k.is_integer())

输出结果：

686291542636760920104201
828427149867.0
True

但实际上k并非整数，浮点数的精度限制导致了误判。

问题背景：欧拉计划第100题 #

题目内容：

一个盒子中有21个彩色圆盘，其中15个蓝色、6个红色。随机取出两个圆盘，取到两个蓝色圆盘的概率为P(BB) = (15/21) × (14/20) = 1/2。
下一个满足随机取两个蓝色圆盘概率恰好为50%的组合是：85个蓝色圆盘和35个红色圆盘。
找到总圆盘数超过10¹²（1000000000000）的第一个组合，求其中蓝色圆盘的数量。

现有代码问题 #

你尝试用decimal模块调整精度，但代码中仍在使用math.sqrt()处理大数，未真正利用decimal的高精度计算能力，同时暴力枚举策略在处理10¹²级别的数时效率极低。

你的完整代码：

import math
from time import perf_counter
from decimal import Decimal, getcontext
getcontext().prec = 50

start = perf_counter()
n = 7 # question specifies 12 (10**12)
limit = 10**n # starts to slow at 7
startno = 1 #292893220000 #2.9*10**11
# squares = set()
# for i in range(1,limit):
#     squares.add(i**2)
candidates = []
for i in range(startno,limit):
    j = (8*i**2 + 1)
    k = math.sqrt(j)
    if k.is_integer() == True:
        b = ((2*i + 1) + k ) / 2
        candidates.append((b,i))
        print(b,i,b/(b+i))
        if b + i > 10**12 and b.is_integer() == True:
            print(b,i)
            break

# print(candidates)
# print(len(squares))
end = perf_counter()
print(end-start, n)

# need to find some way to find the answers. These methods all fail due to rounding in Python

print(math.sqrt(285700000000))
i = 292893227695
j = 8*i**2 + 1
print(j)
k = math.sqrt(j)
print(k)
print(k.is_integer())

解决办法 #

1. 用整数运算验证完全平方数 #

避免依赖浮点数，先计算整数平方根，再验证平方是否等于原数：

import math

i = 292893227695
j = 8 * i**2 + 1
k = math.isqrt(j)  # Python 3.8+，返回不超过sqrt(j)的最大整数
if k * k == j:
    print("是完全平方数")
else:
    print("不是完全平方数")

math.isqrt()专门处理整数，不会有浮点数精度问题。

2. 正确使用decimal模块高精度计算 #

将数值转为Decimal类型后再开平方，判断是否为整数：

from decimal import Decimal, getcontext

getcontext().prec = 50  # 设置足够覆盖大数的精度
i = 292893227695
j = Decimal(8 * i**2 + 1)
k = j.sqrt()
if k == k.to_integral_value():
    print("是完全平方数")
else:
    print("不是完全平方数")

3. 针对欧拉计划第100题的高效解法 #

暴力枚举不可行，题目中的方程8i² + 1 = k²可转化为佩尔方程k² - 8i² = 1，其解可通过递推生成：

• 初始解：k₁=3, i₁=1
• 递推公式：
  kₙ₊₁ = 3kₙ + 8iₙ
  iₙ₊₁ = kₙ + 3*iₙ
• 对应的蓝色圆盘数b=(2i + 1 + k)/2，总圆盘数t=b+i

用递推式可以快速生成所有满足条件的解，直接找到总圆盘数超过10¹²的第一个组合，无需枚举。

示例代码：

# 初始解对应题目中的第一个组合：b=15, i=6（i是红色圆盘数）
k, i = 3, 1
while True:
    b = (2*i + 1 + k) // 2
    total = b + i
    if total > 10**12:
        print("蓝色圆盘数量：", b)
        break
    # 递推生成下一组解
    k, i = 3*k + 8*i, k + 3*i

内容的提问来源于stack exchange，提问作者Dominic