移项可得
$$k^2=\sum_{i=1}^m(n+i)^2-(k-i)^2=\sum_{i=1}^m(n+k)(n-k+2i)=m(n^2-k^2)+(n+k)m(m+1),$$
整理后得一关于$k$的二次方程：
$$(m+1)k^2-m(m+1)k-mn^2-nm(m+1)=0$$
由$n\ge k$可得$k\ge 2m(m+1)$，从而得m的上界70710。此方程的判别式是
$$\Delta=m^2(m+1)^2+4(m+1)(mn^2+nm(m+1))=m(m+1)(m^2+4mn+m+4n^2+4n),$$
需要是完全平方数（否则k不能为整数)。此处需要一些处理的技巧，令
$$\begin{pmatrix}k \\n \end{pmatrix}=\begin{pmatrix}m & 1 \\ m+1 & 1 \end{pmatrix}\begin{pmatrix}x \\t \end{pmatrix}+\begin{pmatrix}m \\ -(m+1) \end{pmatrix},$$
则上述方程可以简化为
$$m(m+1)(x^2-1)=t^2$$
令$m=au^2,m+1=bv^2,t=abuvy$（其中a,b都是无平方因子数），就得$x^2-1=aby^2$，即Pell方程$x^2-aby^2=1$。这一方程的解法已在第66题等题中出现过。
枚举m，并分解出相应的a,b后解此Pell方程，检查每一组解所对应的k是否是整数以及是否小于$10^{10}$并去重后即得结果$238890850232021$。
注：为减少码量，以下使用了sympy作质因数分解并求解Pell方程，由于这一依赖关系，以下为Python 3代码。此外，计算量较大，代码中打印了进度信息。
import math,sympy
from sympy.solvers.diophantine.diophantine import diop_DN
target=10**10
bound=70711
f,g=[0],[0]
for n in range(1,bound+1):
x=y=1
for p,e in sympy.factorint(n).items():
if e%2==1:
x*=p
y*=p**(e//2)
f.append(x)
g.append(y)
print("f,g generated")
tick=bound//100
s=set()
for m in range(1,bound):
d=f[m]*f[m+1]
x0,y0=diop_DN(d,1)[0]
x,y=x0,y0
q=f[m]*f[m+1]*g[m]*g[m+1]
k=m*(x+1)+q*y
while k<=target*2:
n=(m+1)*(x-1)+q*y
if n>=k and n%2==k%2==0:
s.add(k//2)
x,y=x0*x+d*y0*y,x0*y+y0*x
k=m*(x+1)+q*y
if m%tick==0:print(f"{m//tick} percent completed")
print(sum(s))
|
