Pell方程 | 求特殊的二元二次方程的整数解 | 欧拉项目 | 100题

原题链接
如果一个盒子里有21个有色的碟子，15个蓝色的和6个红色的。从中随机取两个，可知取到两个蓝碟子的几率是 P(BB) = (15/21)×(14/20) = 1/2。
下一个满足此条件（即随机取两个碟子的情况下取到两个蓝色碟子的几率是50%）的情况是85个蓝碟子和35个红碟子。
对于包含超过10 ^ 12 = 1,000,000,000,000个碟子的情况中，满足上述条件的并包含最少碟子的情况，该情况下共需要多少个蓝碟子？

假设碟子总数是T，蓝色盘子总数是B，那么它们需要满足(B/T)((B-1)/(T-1)) = 1/2
可以转化为求方程
(T - 1/2) ^ 2 + 2 * (B - 1/2) ^ 2 = -1/4
的整数解。
两边同时乘以4
(T/2 - 1/4) ^ 2 + 2 * (B/2 - 1/4) ^ 2 = -1
令
T = (x+1)/2
B = (y+1)/2
其中，x和y要是奇数
方程转化为
x ^2 - 2 * y ^ 2 = -1
这是一个Pell方程，最小的正整数解是
x = 1
y = 1
根据Pell方程的通解
$$\begin{cases}
x_n=\frac 1 2\left[ {(x_1+\sqrt d y_1)}^{2n-1} + {(x_1-\sqrt d y_1)}^{2n-1}\right] \
y_n=\frac 1 {2\sqrt d }\left[ {(x_1+\sqrt d y_1)}^{2n-1} - {(x_1-\sqrt d y_1)}^{2n-1}\right]
\end{cases}$$
可以算出整数解，不难看出，x(n)和y(n)都是奇数。

用程序计算得到第一个超过2 * 10 ^ 12的x，这样对应满足题意T，通过对应的y计算得到我们要求的答案B。
公式里的d是2。

double sqrt2 = Math.Sqrt(2);

for (int i = 1; ; i++)
{
    long x = (long)((Math.Pow(1 + sqrt2, 2 * i - 1) + Math.Pow(1 - sqrt2, 2 * i - 1)) / 2) + 1;
    long y = (long)((Math.Pow(1 + sqrt2, 2 * i - 1) - Math.Pow(1 - sqrt2, 2 * i - 1)) / 2 / sqrt2) + 1;

    if (x > 2000000000000)
    {
        return (y + 1) / 2;
    }
}

程序中，在计算x和y的时候都+1了，因为前面的double类型运算可能会损失精度，强转成long的时候少了1，所以最后又加了个1。