给一个小数值，能否找到最接近于它的$a+b\sqrt{c}$的有理数$a,b,c$

郝酒

RT，有时候用MMA，算最值的时候不一定能够化简，求数值解是会得到一个近似值，现在的问题是：有没有一个函数或算法，能找出固定结构，$a+b\sqrt{c}$的最接近的近似值.

hbghlyj

下面给出一个 “找 \(a+b\sqrt{c}\) 型最优二次近似” 的完整方案，并说明怎样把参数调到能让
\[x=\pi\;\longmapsto\;\frac{198261+\sqrt{105043517473}}{166274}\]

---

思路

• 若 \(x\approx a+b\sqrt{c}\)，则
  \[(b^2c)\;x^{2}-2ab\,x+(a^{2}-b^{2}c)=0 .\]
  这说明 \(x\) 是某个整数二次多项式的近似根。
• 因此可把问题化成 整数关系 搜索：找到不全为零的整数
  \[(A,B,C)\quad\text{s.t.}\quad A\,x^{2}+B\,x+C\approx0 .\]

整数关系最稳健的实用算法是 PSLQ (1992 Ferguson–Bailey)(en.wikipedia.org)；Python 里 mpmath / SymPy 已内置实现(omz-software.com)。

---

算法步骤
用 Python-mpmath/SymPy 描述；其他语言只要有 PSLQ 或 LLL 子程序即可照搬。

• 输入
  近似实数 \(x\)；工作精度 prec(十进制位)；系数上界 M.
• 设置高精
  mp.dps = prec；y1 = mp.mpf(x); y2 = y1**2.
• PSLQ
  A,B,C = pslq([y2, y1, 1], maxcoeff=M) —— 若返回 None 则提高 prec 或 M.
• 约简
  设 g = gcd(A,B,C)，令 A/=g, B/=g, C/=g；若 A<0 同时变号保证首项正。
• 判别式\(\Delta = B^{2}-4AC\)
  若 \(\Delta<0\) 或不是完全平方，说明取的 M 太小或 prec 不足；回到(1).
• 构造根
  \(\displaystyle x\approx\frac{-B+\sqrt{\Delta}}{2A}\) （符号选能逼近原 \(x\) 的那支）。
• 规范成 \(a+b\sqrt c\)
  把分数约分到最简，写成 \((a+b\sqrt c)/d\).
• 输出
  四元组 \((a,b,c,d)\) 或直接返回字符串。

---

关键参数调优

参数	经验值	说明
prec	\(\lceil 3\log_{10}M \rceil +10\)	PSLQ 需要“每位系数≈3 位十进制精度”作为安全冗余。
M (maxcoeff)	从 \(10^{4}\) 起指数级递增	系数上界过小会找不到关系；过大会拖慢时间。先小后大逐次尝试最稳妥。
容差	tol = 10**(-prec*0.8)	判 Ax²+Bx+C<tol；过松会误判，过紧会漏解。
判别式平方性	先用 is_square(D)	若不是完全平方，可尝试把 \(\gcd(A,B,C)\) 再提取平方因子后重测。

Pi 的特例
目标式中的最大系数 ≈ \(2\times10^5\)，因此取

1. prec = 80
2. maxcoeff = 3*10**5

Copy the Code

就足以让 PSLQ 找到 \((-166274)\,\pi^{2}+\,396522\,\pi-105043517473=0\)。 判别式即 \(\Delta = 396522^{2}-4\cdot(-166274)\cdot(-105043517473)=1.74365\ldots\times10^{11}\) 正好是 \(2\cdot105043517473\) 的平方，故推出
\[  \pi = \frac{198261+\sqrt{105043517473}}{166274}.  \]

---

参考实现（Python 3 + mpmath / SymPy）

1. from mpmath import mp, pslq, mpf
2. from math import isqrt, gcd
4. def quad_approx(x, prec=80, maxcoeff=3*10**5):
5.     mp.dps = prec
6.     x = mp.mpf(x)
7.     rel = pslq([x**2, x, 1], maxcoeff=maxcoeff)
8.     if rel is None:
9.         raise ValueError("Relation not found – raise prec or maxcoeff")
10.     A, B, C = map(int, rel)
11.     g = gcd(gcd(abs(A), abs(B)), abs(C))
12.     A//=g; B//=g; C//=g
13.     if A < 0: A, B, C = -A, -B, -C
14.     D = B*B - 4*A*C
15.     s = isqrt(D)
16.     if s*s != D:
17.         raise ValueError("Discriminant not a perfect square")
18.     num = -B + s if (-B + s)/(2*A) > 0 else -B - s
19.     den = 2*A
20.     g = gcd(num, den)
21.     num//=g; den//=g
22.     a = num // den
23.     b_num = num - a*den
24.     return a, b_num//den, s, den   # (a,b,c,den)  with b=1 here

Copy the Code

测试

1. quad_approx(mp.pi)

Copy the Code

(198261, 1, 105043517473, 166274)

---

进一步改进

• 分层搜索 先用连分数找最佳有理数 \(p/q\) （如 \(355/113\)）；若误差仍在目标容差之外，再上 PSLQ 找二次关系，可明显加快收敛。
• 绝对值筛选 对候选 \((A,B,C)\) 先计算粗略残差 \(|Ax^{2}+Bx+C|\)；若 > 1e-(prec/2) 即可丢弃，减少验证次数。
• 多线程/向量化 对不同上界或不同随机扰动并行启动 PSLQ，提高大常数时的成功率。
• 高阶推广 换成向量 $x^3, x^2, x, 1$ 即可搜索立方型近似；再搭配判别式平方-立方检测即可得到 \(a+b\sqrt[3]{c}\) 等形式。