大分母连分数收敛项计算方法及ContFrac、SymPy库报错求解

我需要从如下公式中求解Dp和Dq：

我尝试使用ContFrac库完成该任务，对应代码如下：

pip install ContFrac

import contfrac
import math

find_number = 2140e225
value = (1,math.sqrt(find_number))
conv = list(contfrac.convergents(value))
print(conv)

但由于目标数值过大，输出结果异常，内容为：

0 + 1/(my_sqrt)

使用SymPy的相关方法也无法正常运行，返回0，对应代码如下：

print(continued_fraction((1//sqrN)))

更新1 #

我测试SymPy时发现应当使用该库自带的sqrt实现，对应代码如下：

from sympy import sqrt, continued_fraction, flatten, continued_fraction_convergents

lm = continued_fraction((1 / sqrt(number)))
lz = flatten(lm)
print(list(continued_fraction_convergents(lz)))

但运行时抛出内部异常：

PrecisionExhausted: Try simplifying the input, using chop=True, or providing a higher maxn for evalf

经排查我无法为continued_fraction操作调高maxn参数值。另外我尝试传入decimal类型数值也失败，对应代码如下：

import decimal
from sympy import continued_fraction

decimal.getcontext().prec = 125
lm = continued_fraction((1 / sqrt(decimal.Decimal(find_number))))

报错信息：

ValueError: expecting a rational or quadratic irrational, not 2.16152455423408e-125

现寻求可正常计算大分母连分数收敛项的可行方案。

方案：基于整数运算手动实现二次无理数连分数迭代 #

你遇到的报错本质是浮点类型无法承载超大数开根后的有效位数，可直接基于二次无理数的连分数迭代规则，全程用整数运算规避精度损失：

• 先将目标数2140e225做等价转换：2140e225 = 214 * 10^226，其平方根为sqrt(214) * 10^113，对应的倒数就是1/(sqrt(214)*10^113)，只需先求sqrt(214)的连分数收敛项，再做简单数值缩放即可，完全不用处理超大浮点数。
• 实现代码如下：

def quadratic_irrational_cont_frac(n, max_terms=100):
    """计算sqrt(n)的连分数项，n为非完全平方正整数，全程整数运算无精度损失"""
    a0 = int(n ** 0.5)
    if a0 * a0 == n:
        return [a0]
    terms = [a0]
    m, d, a = 0, 1, a0
    for _ in range(max_terms - 1):
        m = d * a - m
        d = (n - m * m) // d
        a = (a0 + m) // d
        terms.append(a)
        if a == 2 * a0:
            break
    return terms

def convergents_from_terms(terms):
    """从连分数项生成对应收敛项"""
    convergents = []
    h0, h1 = terms[0], terms[0] * terms[1] + 1
    k0, k1 = 1, terms[1]
    convergents.append((h0, k0))
    convergents.append((h1, k1))
    for a in terms[2:]:
        h = a * h1 + h0
        k = a * k1 + k0
        convergents.append((h, k))
        h0, h1 = h1, h
        k0, k1 = k1, k
    return convergents

# 代入你的业务场景
n = 214
# 生成sqrt(214)的连分数项，可按需调整max_terms获取更多收敛项
cf_terms = quadratic_irrational_cont_frac(n, max_terms=50)
# 生成收敛项，每个收敛项h/k ≈ sqrt(214)，对应1/(sqrt(214)*10^113) ≈ k/(h*10^113)
convs = convergents_from_terms(cf_terms)
# 输出对应Dp和Dq
for idx, (h, k) in enumerate(convs):
    Dp = k
    Dq = h * (10 ** 113)
    print(f"第{idx+1}组收敛项：Dp={Dp}, Dq={Dq}")

该方案全程使用整数运算，不存在浮点精度溢出问题，支持任意长度的连分数项计算，完全适配超大数场景。

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