IBM Ponder this 每个月会出一道谜题。这个月的题目是求所有的整数，它既是一个平方数，又是一个 triangular number （即可以表示为$1+2+\cdots+m$）。
这个题目挺有意思的，难度适中。可以自己先做一做，锻炼锻炼脑子。
下面给一个解答，以及思考的步骤。
首先问题可以转为求解 $ x^2 = \frac{m(m+1)}2$，两边乘以 8 ，即$8x^2+1=(2m+1)^2$。所以问题等价于求解$8x^2+1=y^2$的整数解。
大多数人到这里就卡住了。显然$(0,1)$，$(1, 3)$ 是两组平凡解。后面就没有思路了。我用 Python 写了个枚举程序：
import math
ans = []
for i in range(100000):
a = 8 * i * i + 1
b = math.isqrt(a)
if b * b == a:
ans.append((i, b))
print(ans)
一下子就找到很多组解：(0, 1), (1, 3), (6, 17), (35, 99), (204, 577), (1189, 3363), (6930, 19601), (40391, 114243)。
然后发挥找规律大法，很容易找到$x$的规律：
然后我们找到通解并验证：
满足：
现在我们需要证明，这些是$8x^2+1=y^2$的所有解了！一个自然的想法是，如果我们有一个解$(n,m)$，可以从这里构造一个更小的解。
继续在(0, 1), (1, 3), (6, 17), (35, 99), (204, 577), (1189, 3363), (6930, 19601)里找规律，我们会发现，如果 $(n, m)$ 是一组解，那么 $(3n-m, 3m-8n)$ 就是它前面的解。简单验证：
同理发现 $(3n+m, 3m+8n)$ 就是它后面那一组解。这样所有的解串成一条线，刚好是前面的$(x_n, y_n)$。这样我们就把方程所有的解都找出来了。
回到 IBM 的原题，求$x^2$小于$10^n$的 x 的个数，显然个数即下面的$f(n)$：
def f(n):
return (math.log(4 * math.sqrt(2)) + n / 2 * math.log(10)) / math.log(3 + 2 * math.sqrt(2))
不过题目还要求$f(10^{100})$的精确解。这个上面的表达式可不够，这时候可以用 Python 的高精度计算库mpmath（对， Python 啥都有），即求解下面的g(100)：
import mpmath as mp
def g(n):
# 设置精确度（位数）
mp.mp.dps = n + 2
# 定义方程中用到的常量和变量
a = mp.mpf(4 * mp.sqrt(2))
b = mp.mpf(10)
c = mp.mpf(3 + 2 * mp.sqrt(2))
# 计算方程的解
n = mp.mpf(n)
n = mp.power(10, n)
return (mp.log(a) + n / 2 * mp.log(b)) / mp.log(c)
Q. E. D.
