Project Euler 246

题目

Tangents to an ellipse

A definition for an ellipse is:

Given a circle $c$ with centre $M$ and radius $r$ and a point $G$ such that $d(G,M)<r$, the locus of the points that are equidistant from $c$ and $G$ form an ellipse.

The construction of the points of the ellipse is shown below.

Given are the points $M(-2000,1500)$ and $G(8000,1500)$.

Given is also the circle $c$ with centre $M$ and radius $15000$.

The locus of the points that are equidistant from $G$ and $c$ form an ellipse $e$.

From a point $P$ outside $e$ the two tangents $t_1$ and $t_2$ to the ellipse are drawn.

Let the points where $t_1$ and $t_2$ touch the ellipse be $R$ and $S$.

For how many lattice points $P$ is angle $RPS$ greater than $45$ degrees?

解决方案

$|EG|=|EU|$，$|EU|+|EM|=r$，因此$|EG|+|EM|=r=2a$，其中$2a$是椭圆的长轴长度，焦距$2c=d(G,M)$。由于椭圆的中心也在个点，那么为了方便处理问题，将椭圆的中心挪到原点，那么可以写成椭圆的标准方程：

$$\dfrac{x^2}{a^2}+\dfrac{y^2}{b^2}=1$$

其中$a=7500,c=5000,b^2=a^2-c^2$。

那么，由于椭圆关于两个坐标轴对称，因此只考虑第一象限以及$x,y$轴正半轴上的格点。

对于圆外任意一点$P(x_0,y_0)$，那么该点的点斜式方程为$y-y_0=k(x-x_0)$，那么和椭圆方程联立，得到一个关于$x$的一元二次方程：

$$(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$$

由于所求直线为切线，因此令其$\Delta=0$，得到关于$k$的一元二次方程：

$$(x_0^2 - a^2)k^2-2x_0y_0k+y_0^2-b^2=0$$

也就是说，$k$不同的两个解对应着直线的斜率，不失一般性，假设其分别为$k_0,k_1$。那么这两条直线夹角的正切值为$\left|\dfrac{k_1-k_0}{1+k_0k_1}\right|$。

解决本题时，还采用了一个假设：当$\angle RPS=45°$时，$P$点的轨迹$\Gamma$是一个凸曲线。（如果$\dfrac{a}{b}$值太大，那么$\Gamma$将不会是一个凸曲线，我在这里没有进行太多的研究。）

因此，枚举所有$x_0\ge0$，二分查找$y_0$使其恰好在$\Gamma$外部。那么，$y_0$以下（不包括$y_0$）的所有点都是答案，不过，这需要排除在椭圆内部的点。

代码

from math import ceil, sqrt
from itertools import count

Mx, My = -2000, 1500
Gx, Gy = 8000, 1500
R = 15000
a = R / 2
c = ((Mx - Gx) ** 2 + (My - Gy) ** 2) ** 0.5 / 2
a2 = a ** 2
b2 = a2 - c ** 2
eps = 1e-10
cnt2 = cnt4 = 0
pre_r = int((a + 4) * (a + 4))

for x0 in count(0, 1):
    y0 = ceil(sqrt(b2 * (1 - x0 * x0 / a2)) - eps) if x0 < a else 0
    l = y0
    r = pre_r
    while l < r:
        mid = (l + r) >> 1
        if x0 != a:
            A = x0 * x0 - a2
            B = -2 * x0 * mid
            C = mid * mid - b2
            D = B * B - 4 * A * C
            k = [(-B - sqrt(D)) / (2 * A), (-B + sqrt(D)) / (2 * A)]
            tg = (k[1] - k[0]) / (1 + k[0] * k[1])
        else:
            k = (mid * mid - b2) / (2 * x0 * mid)
            tg = 1 / k
        if 0 < tg <= 1:
            r = mid
        else:
            l = mid + 1
    if l == 0:
        break
    # 统计不同点的个数。
    if x0 == 0:
        cnt2 += l - y0
    elif y0 > 0:
        cnt4 += l - y0
    else:
        cnt4 += l - 1
        if x0 > a:
            cnt2 += 1
    pre_r = l
ans = cnt2 * 2 + cnt4 * 4
print(ans)