Google经典编程竞赛题：计算 $(3 + \sqrt{5})^n$ 的小数点前三位数

By Long Luo

安全界大牛 袁哥 在微博上发布了一道 数学挑战题 ：

我的解答

刷微博时看到这道题目时，我觉得很简单啊，于是马上给出了下面的解答：

令 \(y = (3 + \sqrt{5})^n\) ，两边同取对数， \(\log_{10}{y} = n \log_{10}{(3 + \sqrt{5})}\) ， \(y = 10^{n \log_{10}{5.23607}}\) ，\(\log_{10}{5} \approx 0.7\) ，所以 \(y \approx 10^{0.7n}\) 。

但问题没有这么简单，因为上述解答只在 \(n = 1\) 是正确的，\(n = 2\) ， \(y = 10^{1.4} \approx 25\) 就不对了，因为精度不够！

之后根据微博评论中其他人给的构造共轭数思路，分析出 \(3\) 位数是周期性的，于是又提交了下面的答案：

但问题仍然没有这么简单，因为即使循环周期 \(p = 100\) ，而 双精度浮点数 的有效位数也只有 \(15\) 位，而 \(\sqrt{5}\) 是无理数，同时由于舍入误差， \(\log_{10}{(3 + \sqrt{5})}\) 很快就出现精度不够的问题，得到错误的结果。

之后袁哥发布了 解答 ，图片太大，大家可以点开 图片链接 查看详细解答。

袁哥的题解省略了很多东西，对数学不熟悉的人可能看不太明白，我当时也没有完全看明白。根据袁哥解答我重新写了份题解，整理了思路及缺失的步骤，外加证明，有中学数学水平即可看懂，题解第一部分如下：

Google Code Jam 编程竞赛题

这道题其实是 Google Code Jam 2008 Round 1A Problem C: Numbers ¹ 编程竞赛题，是一道极好的编程练习题。原题如下：

如何求解这道题？

这个问题看似简单，但实际上却是一道很难的编程题。对于题目中的小测试集时，比如 \(2 \le n \le 30\) ，虽然 Java 或 Python 等支持任意精度计算，但直接去计算仍然会失败，大概在 \(n > 17\) 之后就无法得到正确的结果了。问题的关键在于需要避免去计算 \(\sqrt{5}\) ，因为使用浮点数去计算 \(\sqrt{5}\) 的会遇到不仅是精度，还有误差累加导致得到错误的结果。

而使用大数据测试集时，则问题更复杂了，当输入值 \(n\) 接近 \(2\,000\,000\,000\) 时，如果使用之前的方法，你注定会超时而无法 AC 。那么这个问题应该怎么解决呢？

利用共轭构造新数列

这个问题的关键在于想到数学上的共轭( \(\textit{Conjugate}\) )² ，注意到对于复数 \(a+bi\) ，其共轭复数为 \(a-bi\) ，那么 \(3 - \sqrt{5}\) 就是 \(3 + \sqrt{5}\) 的共轭。不妨设 \(\alpha = 3 + \sqrt{5}\) , \(\beta = 3 - \sqrt{5}\) ，构造一个新数列 \(X_n\) ：

\[ X_n = \alpha^n + \beta^n , \quad n \in \mathbb{N} \tag{1} \label{1} \]

注意到 \(X_n\) 是一个整数³ ，也很容易证明。我们使用二项式定理( \(\textit{Binomial theorem}\) ) ⁴ 展开 \(X_n\) ，可得：

\[ X_n = 2 \cdot \sum_{k=0}^{n / 2} \binom{n}{2k} \cdot 3^{n-2 k} \cdot 5^k \tag{2} \label{2} \]

可以看出 \(\sqrt5\) 的奇次项相消了，故 \(X_n\) 为整数得证！

显然有 \(0 < \beta < 1\) , 故有 \(0 < \beta^n < 1\) ，那么 \(X_n - 1 < \alpha^n < X_n\) ，所以问题转化为求 \(X_n\) 的小数点前 \(3\) 位即可。

下面我们将给出这个问题的 \(4\) 种解法，由浅入深，最后的解法肯定会让你 A-Ha 一声的！

解法 \(1\): 系数递推法

容易看出 \(\alpha^n = a_n + b_n \sqrt{5}\) ，这里 \(a_n, b_n\) 均为整数。当然有同学会问，哪里容易看出来啊？不急，我们来简单证明下。

\[ \left ( m + n \sqrt{5} \right ) \cdot \left ( p + q \sqrt{5} \right ) = (mp + 5nq) + (mq + pn) \sqrt{5} \]

同理根据二项式定理则有 \(\beta^n = a_n - b_n \sqrt{5}\) ，那么 \(X_n = 2a_n\) 。

不难得到：

\[ \alpha^{n+1} = (3 + \sqrt{5})(a_n + b_n\sqrt{5}) = (3a_n + 5b_n) + (3b_n + a_n)\sqrt{5} \]

于是有：

\[ a_{n+1} = 3a_n + 5b_n, \ b_{n+1} = 3b_n + a_n. \]

这个递推式有没有想到 斐波那契数(Fibonacci Numbers) ，我们将其写成矩阵形式：

\[ \begin{bmatrix} a_{n+1} \\ b_{n+1} \end{bmatrix} = \begin{bmatrix} 3 & 5 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} a_n \\ b_n \end{bmatrix} \]

令矩阵 \(A = \begin{bmatrix} 3 & 5 \\ 1 & 3 \end{bmatrix}\) ，则有：

\[ \begin{bmatrix} a_n \\ b_n \end{bmatrix} = A \begin{bmatrix} a_{n-1} \\ b_{n-1} \end{bmatrix} =A^n \begin{bmatrix} a_0 \\ b_0 \end{bmatrix} \]

而 \(\alpha^0 = 1\) ，有 \(a_0 = 1, b_0 = 0\) ，我们可以使用 快速幂 算法快速求解，算法实现的 Java 参考代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

private static int[][] matrixMultiplication(int[][] A, int[][] B) {
    if (A == null || A[0] == null || B == null || B[0] == null || A[0].length != B.length) {
        return null;
    }

    int rowA = A.length;
    int colA = A[0].length;

    int rowB = B.length;
    int colB = B[0].length;

    int[][] C = new int[rowA][colB];

    for (int i = 0; i < rowA; i++) {
        for (int j = 0; j < colB; j++) {
            int sum = 0;
            for (int k = 0; k < colA; k++) {
                sum += (A[i][k] * B[k][j]) % 1000;
                sum = sum % 1000;
            }

            C[i][j] = sum % 1000;
        }
    }

    return C;
}

private static int[][] fastExponentiation(int[][] A, int n) {
    if (n == 1) {
        return A;
    }

    if (n % 2 == 0) {
        int[][] A_m = fastExponentiation(A, n / 2);
        return matrixMultiplication(A_m, A_m);
    } else {
        return matrixMultiplication(A, fastExponentiation(A, n - 1));
    }
}

private static String findLast3Digits(int n) {
    int[][] A = {{3, 5}, {1, 3}};

    int[][] A_n = fastExponentiation(A, n);

    int result = (2 * A_n[0][0] + 999) % 1000;

    return String.format("%03d", result);
}

解法 \(2\): 数列递推公式

由于 \(\alpha + \beta = 6, \, \alpha \beta = 4\) ，由韦达定理( \(\textit{Vieta's formulas}\) ) ⁵ ，那么 \(\alpha, \beta\) 是方程 \(x^2 - 6x + 4 = 0\) 的 \(2\) 个根。实际上如果数列满足上述形式，很容易证明其递推公式为： \(a_{n+2} = 6a_{n+1} - 4a_n\) ，这里也简单证明下：

\[ \alpha^{n+1} + \beta^{n+1} = (\alpha + \beta ) \left( \alpha^n + \beta^n \right) - \alpha \beta \left( \alpha^{n-1} + \beta^{n-1} \right ) \]

而 \(\alpha \beta = 4\) ，则有：

\[ X_{n+2} = 6X_{n+1} - 4X_n \]

和解法 1 类似，我们令矩阵 \(B = \begin{bmatrix} 6 & -4 \\ 1 & 0 \end{bmatrix}\) ，则有：

\[ \begin{bmatrix} X_{n+1} \\ X_n \end{bmatrix} = B \begin{bmatrix} X_n \\ X_{n-1} \end{bmatrix} =B^n \begin{bmatrix} X_1 \\ X_0 \end{bmatrix} \]

很容易计算 \(X_0 = 2\) ， \(X_1 = 6\) ，那么同样使用 快速幂 算法求解，算法实现的 Java 参考代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

private static int quickMultiply(int n) {
    if (n == 0) {
        return 2;
    } else if (n == 1) {
        return 6;
    }

    n -= 2;

    int[][] mat = {{6, -4}, {1, 0}};
    int[][] smat = mat.clone();

    while (n > 0) {
        if ((n & 1) == 1) {
            mat = matrixMultiplication(mat, smat);
        }
        smat = matrixMultiplication(smat, smat);
        n >>= 1;
    }

    return (6 * mat[0][0] + 2 * mat[0][1]) % 1000;
}

private static String findLast3Digits(int n) {
    int result = (quickMultiply(n) + 999) % 1000;
    return String.format("%03d", result);
}

解法 \(3\): 利用周期性

因为只需要 \(X_n\) 小数点前 \(3\) 位，对结果取模 \(1000\) ，那么最多只有 \(1000^2 = 10^6\) 种可能，所以末尾 \(3\) 位数一定是周期性的，严谨证明可以参考之前我的手写题解。

对于这道题来说，利用前面的解法，我们可以写代码验证输出，周期是 \(100\) ，从第 \(n = 3\) 开始，所以实际上我们只需要计算前面 \(103\) 个数即可。如果 \(n > 103\) 的话，直接计算 \(X_n = X_{(n - 3) \bmod 100 + 3}\) 。

利用周期性，我们可以解决输入值 \(n\) 很大的情况。

解法 \(4\): 剩余定理

下面我们来揭晓这道题最绝妙的解法，利用剩余定理( \(\textit{Chinese remainder theorem}\) ) ⁶ ，详细过程如下图 3 所示：

得到的最终表达式：

\[ X_n = 3^n \cdot \left[ 752 + n(n+1)(200n^2 + 520) \right ] \bmod 1000 \]

对于这个表达式，我们手算都可以。

课后挑战

通过前面的讲解，你学会了这类题的解法了吗？下面是 \(2\) 个挑战：

1. 计算 \((2 + \sqrt3)^{3000}\) 的小数点前 \(3\) 位数；
2. 计算 \((1 + \sqrt2)^{3000}\) 的小数点后第 \(500 - 502\) 位数。

你学会了吗？

参考资料