半径为R的大圆内部有一半径为r的小圆与其相切并沿着大圆圆周不断滚动，小圆上一定点轨迹所形成的曲线称为圆的内摆线。

以原点为大圆圆心，以x轴正半轴与大圆交点作为此内摆线的起点，则其参数方程为： $$x(t)=(R-t)\cos t+r\cos\left(\frac{R-r}{r}\cdot t\right), \quad y(t)=(R-t)\sin t-r\sin\left(\frac{R-r}{r}\cdot t\right),$$ 令$C(R,r)$为此内摆线上所有坐标为整数且所对应的$\sin t$和$\cos t$均为有理数的点所组成的集合，再令 $$S(R,r)=\sum_{(x,y)\in C(R,r)}|x|+|y|,$$ 以及 $$T(N)=\sum_{R=3}^N\sum_{r=1}^{\lfloor\frac{R-1}{2} \rfloor}S(R,r).$$ 可以验证$C(3, 1)=\{(3, 0), (-1, 2), (-1,0), (-1,-2)\}$以及$C(2500, 1000)$包含(2500, 0), (772, 2376), (772, -2376), (516, 1792),(516, -1792), (500, 0), (68, 504), (68, -504),(-1356, 1088), (-1356, -1088), (-1500, 1000), (-1500, -1000)。

注意：$(-625, 0)\notin C(2500, 1000)$，因为该点所对应的$\sin t$和$\cos t$不是有理数。

因此$S(3,1)=(|3|+|0|)+(|-1|+|2|)+(|-1|+|0|)+(|-1|+|-2|)=10$，且可以进一步验证$T(3)=10, T(10)=524, T(100)=580442, T(1000)=583108600$。

求$T(10^6)$。

本题难度:

解答

本题的计算很繁琐，以下仅给出主要步骤，这些步骤系摘自官方论坛。

需要找出能令$\sin(At), \sin(Bt), \cos(At), \cos(Bt)$都为有理数且能令 $$x(t)=Acos(Bt) + Bcos(At), \quad y(t)=Asin(Bt)-Bsin(At), $$ 取到整数并满足$A>B\ge$, $A+B\le N$的参数$A,B,t$。这样的解可以由A,B互素的基础解系生成，并在统计时按t是否是$\pi/2$的整数倍区分计算。

最终结果是$583333163984220940$。

from fractions import gcd

N=10**6
M=10**3
bound=9
    
T=0

for u in range(2,M):
    for v in range(u%2+1,u,2):
        if gcd(u,v)>1:
            continue
        c=u*u+v*v
        if 3*c*c>N:
            break
        a,b=u*u-v*v,2*u*v
        for p in range(2,bound):
            for q in range(1,p):
                if gcd(p,q)>1:
                    continue
                r=(p+q)*c**p
                for z in (complex(a,b),complex(-a,b),complex(b,a),complex(-b,a)):
                    w=p*c**(p-q)*z**q+q*z.conjugate()**p
                    x,y=int(round(w.real)),int(round(w.imag))
                    s=r
                    while x%c==y%c==0:
                        x//=c
                        y//=c
                        s//=c
                    k=N//s
                    T+=k*(k+1)*(abs(x)+abs(y))

f=[1]*(N+1)  
for n in range(2,N+1,2):
    f[n]=2*f[n//2]
    
for R in range(3,N+1):
    rMax=(R-1)//2  
    T+=4*R*rMax    
    k=rMax//f[R]
    s=k*(k+1)//2*f[R]
    T-=2*s  
    if R%2==0:
        x=f[R]//2
        k=rMax//x
        T-=4*(k*(k+1)//2*x-s)  
      
print T