将等式两边立方可得
$$2a+3\sqrt[3]{a^2-b^2c}\cdot(\sqrt[3]{a+b\sqrt{c}}+\sqrt[3]{a-b\sqrt{c}})=1,$$
代入原等式，移项后再立方得
$$(1-2a)^3=27(a^2-b^2c),$$
即
$$8a^3+15a^2+6a-1=27b^2c,$$
等式左侧可以进一步分解得
$$(a+1)^2(8a-1)=27b^2c.$$
等式右侧能被3整除，而左侧只有在a模3余2时能被3整除，因此设$a=3k+2$可进一步得
$$(k+1)^2(8k+5)=b^2c$$
令$d=\gcd(k+1,b^2)$，并记$p=(k+1)/d$，$q=b/d$，则有
$$p^2(8k+5)=q^2c,$$
由于p,q互素， 因此存在自然数m满足
$$m=\frac{8k+5}{q^2}=\frac{c}{p^2},$$
即
$$\begin{cases}a=3k+2 \\ b=dq \\ c=mp^2 \\ k+1=dp \\ 8k+5=mq^2 \end{cases}$$
整理后可得这样三个条件：
$$\begin{cases} 8dp-mq^2=3 \\ \gcd(p,q)=1 \\ 3dp+dq+mp^2\le 110000001 \end{cases},$$
枚举p,q，对每组p,q，用辗转相除法计算$8dp-mq^2=3$的一组基本解$(d_0,m_0)$，则其通解的形式为
$$d_t=tq^2+d_0, \quad m_t=8tp+m_0.$$
从而有
$$N+1\ge 3d_tp+d_tq+m_tp^2=3tq^2p+3d_0p+tq^3+d_0q+8tp^3+m_0p^2$$
即
$$t\le\frac{110000001-3d_0p-m_0p^2-d_0q}{3pq^2+q^3+8p^3}.$$
此外，由此丢番图方程的形式还可以看出q需要是奇数，否则$\gcd(8p,q^2)$将是偶数，而等式右侧的3是奇数。
最后还需确定q,p的范围，事实上只需要确定q的范围，因为对每个固定的q，有
$$k=\frac{mq^2-5}{8}$$
以及
$$a+b+c=3k+2+dq+mp^2=\frac{3}{8}(mq^2-5)+2+dq+p^2,$$
取$m=d=1$，并令上式不超过110000000就可得p的上界。
容易验证当$a=3k+2$固定时k也固定，此时在约束条件$b^2c=(k+1)^2(8k+5)$下有$b=2c$时取得最小值（用拉格郎日乘数法，分别求$b^2c$和$b+c$对b,c各自的偏导可得极值点满足$2bc=b^2$）。
在这一条件下令$a+b+c$不超过110000000，就可得k不超过36666665，而$q^2$是$8k+5$的约数，因此q不超过18000（这是个很粗糙的上界，程序运行过程显示q不超过8000）。
最后结果是$18946051$。
注：计算量较大， 代码中打印了进度信息。
def extGCD(a,b):
if b==0:return a,1,0
d,x,y=extGCD(b,a%b)
return d,y,x-(a/b)*y
def bezout(a,b,c):
d,x,y=extGCD(a,b)
if c%d>0:return -1,x,y
y=y*(c/d)%abs(a/d)
x=(c-b*y)/a
return abs(d),x,y
N=110000000
r=0
for q in range(1,18000,2):
p=1
while 3*q*q+8*p*p+8*q+1<=8*N:
x,d,m=bezout(8*p,-q*q,3)
s=N+1-3*d*p-m*p*p-d*q
if x==1 and s>=0:
r+=s/(3*p*q*q+q*q*q+8*p*p*p)+1
p+=1
if (q-1)%200==0: print (q-1)/200,"percent completed."
print r
|
