固定k，统计在每一个点处能取到45度角的三角形的个数，再减去在该点处能取到90度角的等腰直角三角形的个数，就是固定该k值时满足要求的四元整数组的总数。
用余弦定理可以计算出两种情况的a,b,c（a处作为顶点）都可以表示为$a,x-a,y-a$的形式，不过前者x可以取$k(y-k)/(y+k)$或$k(y+k)/(k-y)$，且a的取值取决于$xy+k^2$的符号或者在区间$[x-a,y-a]$以外或者在其内，此时y是$2k^2$的约数。而后则$x=-k^2/y$，a则可以取$(ky-x^2)/(2(k\pm x))$，且此时x,y都是$k^2$的约数。
因此先用筛法找出$k^2$和$2k^2$的约数列表，再根据上面的结论统计范围内吗满足要求的三元组，最后汇总即得结果$141630459461893728$。
注：为便于整除，以下代码为Pyhthon 3，且代码中打印了进度信息。此外，本题中的约数列表的生成是先维护每个数的最小素因子再递归计算，先可以先维护最大素因子（见第407题）。
target=1000001
bound=1000000000
minP=[0]*target
tick=10000
f=lambda x,y,p,q:0 if q==0 or p%q else (p//q)*2!=x and (p//q)*2!=y and max(abs(p//q),abs(x-p//q),abs(y-p//q))<=bound
g=lambda k,x,y:0 if x>=y else (max(min(min(bound,x+bound),(y-1)//2)-max(max(-bound,y-bound),x//2+1)+1,0) if k*k+x*y<0 else max(min(bound,x+bound)-max(-bound,y-bound)+1,0)-max(min(min(bound,x+bound),y//2)-max(max(-bound,y-bound),(x+1)//2)+1,0))
for i in range(2,target):
if not minP[i]:
for j in range(i,target,i):
if not minP[j]:
minP[j]=i
ans=0
for k in range(1,target):
d=[1]
x=k
while x>1:
j,p=0,minP[x]
while x%p==0:
x//=p
j+=1
t=len(d)
for i in range(t):
q=1
for h in range(2*j):
q*=p
d.append(d[i]*q)
ans-=sum(f(d[i],-k*k//d[i],-(k*k//d[i])*k-d[i]*d[i],2*(k-d[i]))+f(d[i],-k*k//d[i],-(k*k//d[i])*k+d[i]*d[i],2*(k+d[i])) for i in range(len(d)))
t=len(d)
for i in range(t):
if (k*k//d[i])%2:
d.append(2*d[i])
t=len(d)
for i in range(t):
d.append(-d[i])
ans+=sum(g(k,k*(i-2*k)//i,i-k)+g(k,k*(i+2*k)//(-i),i+k) for i in d)
if k%tick==0:
print(k//tick,"percent completed")
print(ans)
|
