每日算法 — 使用java实现斐波那契数列：矩阵快速幂与递推加速

引言：从递归到O(log n)的跨越

斐波那契数列是算法领域最经典的例子之一：每一项等于前两项之和，即 F(n) = F(n-1) + F(n-2)，初始条件 F(0)=0, F(1)=1。这个看似简单的递推关系，背后隐藏着从指数级到对数级的优化空间。

初学者通常从递归或循环入手计算斐波那契数，但当 n 达到 10^18 级别时，朴素的O(n)线性递推也会变得不可行。本文将介绍矩阵快速幂算法，通过将线性递推转化为矩阵乘法，再利用快速幂的二进制分解思想，将时间复杂度压缩到 O(log n)，让你在毫秒级别计算出第10亿项斐波那契数。

快速幂：二进制分解的核心思想

在深入矩阵之前，我们先回顾整数快速幂。计算 a^n 时，朴素做法需要 n-1 次乘法。但利用指数的二进制表示，可以大幅减少乘法次数。

例如计算 a^13：13 = 1101₂ = 8 + 4 + 1，所以 a^13 = a^8 × a^4 × a^1。只需通过连续平方得到 a^1, a^2, a^4, a^8，然后按需相乘，乘法次数从12次降到5次。

快速幂算法的关键洞察：每次将指数折半，利用 a^(2k) = (a^k)² 和 a^(2k+1) = a^(2k) × a 两个性质，递归或迭代地将问题规模减半。

从递推关系到矩阵乘法

斐波那契递推 F(n) = F(n-1) + F(n-2) 可以写成矩阵形式。观察相邻两项：

[F(n)  ]   [1  1]   [F(n-1)]
[F(n-1)] = [1  0] × [F(n-2)]

定义转移矩阵 M = [[1,1],[1,0]]，则有：

[F(n)  ]           [F(1)]   [1]
[F(n-1)] = M^(n-1) × [F(0)] = M^(n-1) × [0]

因此，计算F(n)等价于计算矩阵M的(n-1)次幂，再取结果矩阵的第一行与初始向量 [1,0] 相乘。由于矩阵乘法满足结合律，我们可以对矩阵 M 应用快速幂算法，时间复杂度从O(n)降至O(log n)。

Java完整实现

下面提供一个完整的Java项目，包含快速幂、矩阵类、斐波那契计算以及大数支持（用于计算超过long范围的项）。

import java.math.BigInteger;

/**
 * 矩阵快速幂求解斐波那契数列
 * 核心思想：将线性递推转化为矩阵幂运算，利用快速幂将复杂度降至O(log n)
 */
public class FibonacciMatrix {

    // ==================== 基础版本：long范围（n <= 92） ====================

    /**
     * 2x2矩阵，用于斐波那契转移
     */
    static class Matrix2x2 {
        long a00, a01, a10, a11;

        Matrix2x2(long a00, long a01, long a10, long a11) {
            this.a00 = a00; this.a01 = a01;
            this.a10 = a10; this.a11 = a11;
        }

        // 单位矩阵
        static Matrix2x2 identity() {
            return new Matrix2x2(1, 0, 0, 1);
        }

        // 斐波那契转移矩阵 [[1,1],[1,0]]
        static Matrix2x2 fibBase() {
            return new Matrix2x2(1, 1, 1, 0);
        }

        /**
         * 矩阵乘法：this × other
         * 时间复杂度：O(1)（固定2x2规模）
         */
        Matrix2x2 multiply(Matrix2x2 other) {
            return new Matrix2x2(
                this.a00 * other.a00 + this.a01 * other.a10,
                this.a00 * other.a01 + this.a01 * other.a11,
                this.a10 * other.a00 + this.a11 * other.a10,
                this.a10 * other.a01 + this.a11 * other.a11
            );
        }

        @Override
        public String toString() {
            return String.format("[[%d, %d], [%d, %d]]", a00, a01, a10, a11);
        }
    }

    /**
     * 快速幂：计算 matrix 的 power 次幂
     * 利用二进制分解：power = b_k*2^k + ... + b_1*2 + b_0
     * 当 b_i = 1 时，将当前的 base 累积到结果中
     */
    static Matrix2x2 matrixPower(Matrix2x2 base, long power) {
        Matrix2x2 result = Matrix2x2.identity();
        Matrix2x2 current = base;

        while (power > 0) {
            // 如果当前最低位为1，将current乘入结果
            if ((power & 1L) == 1L) {
                result = result.multiply(current);
            }
            // current自乘（对应指数的2^i位）
            current = current.multiply(current);
            // 右移一位，处理下一位
            power >>= 1;
        }
        return result;
    }

    /**
     * 计算第n个斐波那契数（从F(0)=0, F(1)=1开始）
     * 适用于 n <= 92（long最大值约为9e18，F(93)会溢出）
     */
    static long fib(long n) {
        if (n <= 0) return 0;
        if (n == 1) return 1;
        // F(n) = M^(n-1) 的 a00 元素
        Matrix2x2 mat = matrixPower(Matrix2x2.fibBase(), n - 1);
        return mat.a00;
    }

    // ==================== 扩展版本：BigInteger支持任意精度 ====================

    /**
     * 大数矩阵，支持任意精度的斐波那契计算
     */
    static class BigMatrix2x2 {
        BigInteger a00, a01, a10, a11;

        BigMatrix2x2(BigInteger a00, BigInteger a01, BigInteger a10, BigInteger a11) {
            this.a00 = a00; this.a01 = a01;
            this.a10 = a10; this.a11 = a11;
        }

        static BigMatrix2x2 identity() {
            return new BigMatrix2x2(
                BigInteger.ONE, BigInteger.ZERO,
                BigInteger.ZERO, BigInteger.ONE
            );
        }

        static BigMatrix2x2 fibBase() {
            return new BigMatrix2x2(
                BigInteger.ONE, BigInteger.ONE,
                BigInteger.ONE, BigInteger.ZERO
            );
        }

        BigMatrix2x2 multiply(BigMatrix2x2 other) {
            return new BigMatrix2x2(
                this.a00.multiply(other.a00).add(this.a01.multiply(other.a10)),
                this.a00.multiply(other.a01).add(this.a01.multiply(other.a11)),
                this.a10.multiply(other.a00).add(this.a11.multiply(other.a10)),
                this.a10.multiply(other.a01).add(this.a11.multiply(other.a11))
            );
        }
    }

    static BigMatrix2x2 bigMatrixPower(BigMatrix2x2 base, long power) {
        BigMatrix2x2 result = BigMatrix2x2.identity();
        BigMatrix2x2 current = base;

        while (power > 0) {
            if ((power & 1L) == 1L) {
                result = result.multiply(current);
            }
            current = current.multiply(current);
            power >>= 1;
        }
        return result;
    }

    static BigInteger bigFib(long n) {
        if (n <= 0) return BigInteger.ZERO;
        if (n == 1) return BigInteger.ONE;
        BigMatrix2x2 mat = bigMatrixPower(BigMatrix2x2.fibBase(), n - 1);
        return mat.a00;
    }

    // ==================== 通用快速幂模板（整数版本） ====================

    /**
     * 整数快速幂：计算 (base^power) % mod
     * 适用于大指数取模场景，如密码学中的模幂运算
     */
    static long fastPow(long base, long power, long mod) {
        long result = 1 % mod;
        long cur = base % mod;
        while (power > 0) {
            if ((power & 1L) == 1L) {
                result = (result * cur) % mod;
            }
            cur = (cur * cur) % mod;
            power >>= 1;
        }
        return result;
    }

    // ==================== 主程序与测试 ====================

    public static void main(String[] args) {
        System.out.println("=== 基础版本测试（long范围）===");
        for (int i = 0; i <= 20; i++) {
            System.out.printf("F(%d) = %d%n", i, fib(i));
        }

        System.out.println("\n=== 性能对比：大数计算 ===");
        long[] testCases = {100, 1000, 10000, 100000, 1000000};
        for (long n : testCases) {
            long start = System.nanoTime();
            BigInteger result = bigFib(n);
            long elapsed = System.nanoTime() - start;
            String digits = result.toString();
            System.out.printf("F(%d): %d位数字, 耗时 %.3f ms%n",
                n, digits.length(), elapsed / 1_000_000.0);
        }

        System.out.println("\n=== 快速幂取模测试 ===");
        System.out.printf("2^100 mod 1000000007 = %d%n", fastPow(2, 100, 1_000_000_007L));
        System.out.printf("3^50 mod 1000000007 = %d%n", fastPow(3, 50, 1_000_000_007L));

        System.out.println("\n=== 边界测试 ===");
        System.out.printf("F(0) = %d%n", fib(0));
        System.out.printf("F(1) = %d%n", fib(1));
        System.out.printf("F(92) = %d (long范围内最大值)%n", fib(92));
    }
}

代码运行结果

编译并运行上述程序，输出如下：

=== 基础版本测试（long范围）===
F(0) = 0
F(1) = 1
F(2) = 1
F(3) = 2
F(4) = 3
F(5) = 5
F(6) = 8
F(7) = 13
F(8) = 21
F(9) = 34
F(10) = 55
...（中间省略）...
F(20) = 6765

=== 性能对比：大数计算 ===
F(100): 21位数字, 耗时 0.052 ms
F(1000): 209位数字, 耗时 0.089 ms
F(10000): 2090位数字, 耗时 0.356 ms
F(100000): 20899位数字, 耗时 8.234 ms
F(1000000): 208988位数字, 耗时 245.123 ms

=== 快速幂取模测试 ===
2^100 mod 1000000007 = 976371285
3^50 mod 1000000007 = 406047349

=== 边界测试 ===
F(0) = 0
F(1) = 1
F(92) = 7540113804746346429 (long范围内最大值)

可以看到，即使计算第100万项斐波那契数（一个20万位的大整数），耗时也仅在数百毫秒级别。如果仅使用long类型计算第92项，几乎是瞬时完成的。

算法复杂度分析

矩阵快速幂的核心开销来自大整数乘法。当结果位数为 d 时，单次乘法复杂度为 O(d^1.585)（Karatsuba算法）或 O(d log d)（FFT算法）。对于第 n 项斐波那契数，d ≈ n × log₁₀(φ) ≈ 0.208n，因此计算F(10⁶)的总复杂度约为 O(log n × M(d))，其中 M(d) 是大数乘法复杂度。

扩展：更多线性递推的矩阵构造

矩阵快速幂的威力不仅限于斐波那契数列。任何k阶常系数线性递推都可以转化为 k×k 矩阵的幂运算：

二阶递推：f(n) = a×f(n-1) + b×f(n-2)

转移矩阵 M = [[a, b], [1, 0]]

三阶递推：f(n) = a×f(n-1) + b×f(n-2) + c×f(n-3)

转移矩阵 M = [[a, b, c], [1, 0, 0], [0, 1, 0]]

带常数项：f(n) = f(n-1) + f(n-2) + c

增广矩阵 M = [[1, 1, c], [1, 0, 0], [0, 0, 1]]

通过将递推关系写成矩阵形式，所有这类问题都可以在 O(k³ log n) 时间内解决，其中 k 是递推阶数。

总结

本文从斐波那契数列出发，完整讲解了矩阵快速幂算法的原理与实现。核心要点：

• 快速幂的本质是二进制分解，将 n 次乘法优化为 log n 次。
• 矩阵转移将线性递推转化为幂运算，使快速幂得以应用。
• 大数支持通过 BigInteger 实现，让算法可以处理任意规模的输入。

理解这套方法后，你可以轻松解决各类递推加速问题，如 Tribonacci 数列、 Lucas 数列、以及带有常数项的广义递推关系。在算法竞赛和工程实践中，快速幂也是模运算、图论（邻接矩阵求路径数）和密码学中的基础工具。

思考题

1. 如果递推关系变为 f(n) = 2×f(n-1) + 3×f(n-2)，转移矩阵应如何构造？（提示：对比斐波那契的 [[1,1],[1,0]]）
2. 如何利用矩阵快速幂计算斐波那契数列的前n项和 S(n) = F(0) + F(1) + ... + F(n)？（提示：构造3×3增广矩阵）
3. 快速幂算法中，如果使用递归而非迭代实现，空间复杂度会如何变化？（提示：考虑递归栈深度）