由构造可知$EG+EM=EU+EM=r$，因此G和M就是椭圆的两个焦点，长轴长度为$r=15000$，焦距为$MG=10000$。
椭圆的中心$(3000,1500)$也是整点，因此将坐标原点平移到该点不会影响结果。如此则原问题转换为给定标准方程为
$$\frac{x^2}{7500^2}+\frac{y^2}{7500^2-5000^2}=1,$$
的椭圆，从椭圆外整点向其引两条切线，求满足切线夹角大于45度的整点数量。
由对称性，只需考虑第一象限和两坐标轴正半轴的情况，依次枚举P的横坐标$x_0$，对每个$x_0$，选择一个合适的范围，以下使用的下界是0或椭圆的边界（取决于$(x_0,0)$是否在椭圆内），上界初值是$(r/2)^2$，并随循环的推进而依次减少。
在这一范围内二分查找满足要求的总坐标$y_0$的数量，直到没有满足要求的$y_0$为止。
设$P=(x_0,y_0)$是一个满足要求的点，过该点的直线具有$y-y_0=k(x-x_0)$的形式，代入椭圆方程得
$$(a^2k^2+b^2)x^2+(2a^2ky_0x-2a^2k^2x_0)x+a^2y_0^2-2a^2kx_0y_0+a^2k^2x_0^2-a^2b^2=0,$$
切线对应该方程作为关于x的方程只有一个解的情况，即$\Delta=0$，计算化简得
$$(x_0^2 -a^2)k^2-2x_0y_0k+y_0^2-b^2=0,$$
这一关于k的方程的两个根$k_0,k_1$就是切线的斜率，切线夹角的正切是$(k_1-k_0)/(1+k_0k_1)$，需要该值大于1。
最终结果是$810834388$。
注：以下代码为Python3，因需要使用除法的新特性以及需要math.ceil返回整数（Python2中math.ceil的返回值是浮点数）。使用Python2无法得到正确结果。
import math
a=7500
a2=7500*7500
c2=5000*5000
b2=a2-c2
u2=u4=0
r0=a2
x=0
while True:
y=math.ceil(math.sqrt(b2*(1-x*x/a2))) if x < a else 0
left=y
right=r0
while left < right:
middle=(left+right)//2
if x!=a:
A=x*x-a2
B=-2*x*middle
C=middle*middle-b2
D=B*B-4*A*C
k0,k1=(-B-math.sqrt(D))/(2*A), (-B+math.sqrt(D))/(2*A)
t=(k1-k0)/(1+k0*k1)
else:
t=(2*x*middle)/(middle*middle-b2)
if 0 < t <= 1:
right=middle
else:
left=middle+1
if left==0:
break
if x==0:
u2+=left-y
elif y>0:
u4+=left-y
else:
u4+=left-1
u2+=(x>a)
r0=left
x+=1
print(2*u2+4*u4)
|
