偶斐波那契数列性质与欧拉计划第2题 Properties of Even Fibonacci numbers and Project Euler problems 2

Problem 2 Even Fibonacci numbers

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:
1, 2, 3, 5, 8, 13, 21, 34, 55, 89, …
By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

问题 2 偶数斐波那契数

斐波那契数列中的每个新项都是通过添加前两项来生成的。从 1 和 2 开始，前 10 项将是：
1, 2, 3, 5, 8, 13, 21, 34, 55, 89, …
求其中小于四百万且为偶数的项之和。

目前在网上找到的大部分资料都是通过穷举斐波那契数并求和的, 这个方法针对题中 4M 以内的数据量足够迅速, 但是显然有更高效的方案.

下面将介绍简单介绍斐波那契数列的性质, 较为详细地讲解偶斐波那契数列的性质, 最后给出一个偶斐波那契数列的求和公式; 对于重要的算法, 给出rust代码示例.

如果只需要结论, 请直接看最后一节.

斐波那契数列

斐波那契数列后文记为 \(F(n)\), 定义为

\[F(n)=\begin{cases} 0 & \text{if } n = 0 \\ 1 & \text{if } n = 1 \\ F(n - 1) + F(n - 2) & \text{if } n > 1 \end{cases} \]

欧拉计划第2题中的斐波那契数列由1开始而非从0开始, 但是从数学角度讲从0开始更方便计算.

常用的求 \(F(n)\) 算法包括递归法, 动态规划法, 公式法, 矩阵幂法, 快速倍增法. 相关算法已有很多文章介绍, 以下仅给出示例代码

1. 递归法

fn fib(n: u64) -> u64 {
    match n {
        0 => 0,
        1 => 1,
        _ => fib(n - 1) + fib(n - 2)
    }
}

2. 动态规划法

通过dp数组记录已经计算的结果以减少重复计算, 同时 \(F(n)\) 只需要用 \(F(n-1)\) 和 \(F(n-2)\) 求得, 故可只保留最后两项以减少空间占用. 代码展示用dp数组的方案.

fn fib(n: u64) -> u64 {
    let n = n as usize;
    let mut dp = Vec::with_capacity(n + 1);
    dp.push(0);
    dp.push(1);
    for i in 1..n {
        dp.push(dp[i] + dp[i-1])
    }
    dp[n]
}

3. 公式法

利用通项公式

\[F(n) = \frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2})^n-(\frac{1-\sqrt{5}}{2})^n] \]

其虽然看似复杂度仅为 \(O(1)\), 但求幂过程导致实际复杂度为 \(O(\log n)\), 且存在精度问题, 不适合在n很大时使用. 其中 \((\frac{1-\sqrt{5}}{2})^n\) 恒小于1, n \(\geq\) 2开始小于0.5, 故也可使用

\[F(n) = round(\frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2})^n) \]

注意该公式前几项需要单独处理. 代码只展示原始通项公式.

fn fib(n: u64) -> u64 {
    (1_f32 / 5_f32.sqrt()
        * (((1_f32 + 5_f32.sqrt()) / 2_f32).powi(n as i32)
            - ((1_f32 - 5_f32.sqrt()) / 2_f32).powi(n as i32)))
    .round() as u64
}

4. 矩阵幂法

利用公式

\[\begin{bmatrix} F(n)\\ F(n+1) \end{bmatrix} = \begin{bmatrix} 0 & 1\\ 1 & 1 \end{bmatrix}^{n} \begin{bmatrix} F(0)\\ F(1) \end{bmatrix} \]

或可作

\[\begin{bmatrix} F(n-1) & F(n)\\ F(n) & F(n+1)\\ \end{bmatrix} = \begin{bmatrix} 0 & 1\\ 1 & 1 \end{bmatrix}^{n} \]

可利用矩阵求幂来计算, 同时求幂过程可用快速幂加速至 \(O(\log n)\).

use std::ops::Mul;

#[derive(Copy,Debug, Clone)]
struct Mat2x2(u64, u64, u64, u64);

impl Mul for Mat2x2 {
    type Output = Self;
    fn mul(self, rhs: Self) -> Self::Output {
        Self(
            self.0 * rhs.0 + self.1 * rhs.2,
            self.0 * rhs.1 + self.1 * rhs.3,
            self.2 * rhs.0 + self.3 * rhs.2,
            self.2 * rhs.1 + self.3 * rhs.3,
        )
    }
}

impl Mat2x2 {
    fn pow(mut self, mut exp: u64) -> Self {
        let mut ans = Self(1, 0, 0, 1);
        while exp > 0 {
            if exp & 1 == 1 {
                ans = ans * self;
            }
            self = self * self;
            exp >>= 1;
        }
        ans
    }
}

fn fib(n: u64) -> u64 {
    let ori = Mat2x2(0, 1, 1, 1);
    ori.pow(n).1
}

5. 快速倍增法

通过矩阵幂可推导出以下公式

\[\begin{cases} F(0) &= 0\\ F(1) &= 1\\ F(2n) &= F(n)(F(n) + 2 F(n-1))\text{ if }n>1 \\ F(2n+1) &= F(n)^2 + F(n+1)^2\text{ if } n > 1 \end{cases} \]

fn fib(n: u64) -> u64 {
    match n {
        0 => 0,
        1 => 1,
        _ => match n % 2 {
            0 => fib(n / 2) * (fib(n / 2) + 2 * fib(n / 2 - 1)),
            1 => fib(n / 2 + 1).pow(2) + fib(n / 2).pow(2),
            _ => panic!(),
        },
    }
}

通过优化能消除其中的重复计算和递归:

fn fib_pair(n: u64) -> (u64, u64) {
    if n == 0 {
        (0, 1)
    } else {
        let (mut l, mut h) = (0u64, 1u64);
        let bit_length = n.ilog2();
        let mut mask = 1 << bit_length;
        while mask != 0 {
            (l, h) = (l * (2 * h - l), l * l + h * h);
            if n & mask != 0 {
                (l, h) = (h, l + h);
            }
            mask >>= 1;
        }
        (l, h)
    }
}
 
fn fib(n: u64) -> u64 {
    fib_pair(n).0
}

通过复杂度分析可以发现其与矩阵快速幂一样达到\(O(\log n)\), 但快速倍增法优化后的常数要小于矩阵快速幂. 如果对斐波那契数列算法感兴趣, 可以移步到这篇文章, 这里会进一步介绍优化的方案，并提出了一个新的 \(O(\log n)\) 算法: 扩域法.

偶斐波那契数列

由斐波那契数列中的偶数构成的数列, 下文用\(G(n)\)表示.

\[\begin{align*} &G(0) = 0 = F(0),\\ &G(1) = 2 = F(3),\\ &G(2) = 8 = F(6),\\ &G(3) = 34 = F(9),\\ &G(4) = 144 = F(12)\dots \end{align*} \]

显然有 \(G(n)=F(3n)\), 可以用求和的奇偶性来证明.

同时 \(G(n)\) 满足以下递推公式

\[G(n)=\begin{cases} 0 & \text{if } n = 0 \\ 2 & \text{if } n = 1 \\ 4G(n - 1) + G(n - 2) & \text{if } n > 1 \end{cases} \]

该递推公式可由 \(F(n)=F(n-1)+F(n-2)\) 和 \(F(n)=F(n+1)-F(n-1)\) 推导出来:

\[\begin{align*} F(n)&=F(n-1) + F(n-2)\\ &=F(n-2) + F(n-3) + F(n-2)\\ &=2F(n-2) + F(n-3)\\ &=2(F(n-3) + F(n-4)) + F(n-3)\\ &=3F(n-3) + 2F(n-4)\\ &=3F(n-3) + 2(F(n-3) - F(n-5))\\ &=3F(n-3) + 2(F(n-3) - F(n-3) + F(n-6))\\ &=4F(n-3) + F(n-6) \end{align*} \]

既然\(G(n)\)与\(F(n)\)有相似的递推公式, 两者就有相似的性质, 以上所有算法都适用于\(G(n)\).

1. 递归法

fn even_fib(n: u64) -> u64 {
    match n {
        0 => 0,
        1 => 2,
        _ => even_fib(n - 1) * 4 + even_fib(n - 2),
    }
}

2. 动态规划法

fn even_fib(n: u64) -> u64 {
    let n = n as usize;
    let (mut f1, mut f2) = (0, 2);
    for _ in 0..n {
        let f3 = f1 + f2 * 4;
        f1 = f2;
        f2 = f3;
    }
    f1
}

3. 公式法

所有形如

\[f(n)=af(n - 1) + bf(n - 2) \]

的递推公式都存在通项公式, 通过以下方法求解

1. 建立并求解特征方程
   通过递推式得到特征方程为:$$r^2 - ar - b = 0$$
   假设求得解为 \(r_1\) 与\(r_2\).
2. 通过根确定通项公式
   若存在两实根 \(r_1\neq r_2\) 则通项公式为: $$f(n)=A r_1^n + Br_2^n$$
   若 \(r_1 = r_2\) 则通项公式为: $$f(n)=(A+Bn)r^n$$
   若存在一对共轭复根 \(r_{1,2} = a\pm bi\) 则通项公式仍为: $$f(n)=A r_1^n + Br_2^n$$
   将复数改写为指数形式并通过棣莫弗公式简化, 可得到三角形式: $$f(n) = r^n((A+B)\cos n\theta + i(A-B)\sin n\theta)$$
   其中

\[\begin{align*}r &= \sqrt{a^2 + b^2}\\ \theta &= \arctan\frac{b}{a}\end{align*} \]

3. 确定系数
   注意到2中所有通项公式均存在两个待定系数A和B, 故需要带入两个边界条件求解得到

通过上述三个步骤可得 \(G(n)\) 通项公式为

\[G(n) = \frac{1}{\sqrt{5}}[(2 + \sqrt{5})^n - (2 - \sqrt{5})^n] \]

同样存在精度问题, 其中 \((2 - \sqrt{5})^n\) 恒小于0.5, 也可省略并四舍五入. 代码展示原始通项公式

fn even_fib(n: u64) -> u64 {
    (1_f32 / 5_f32.sqrt()
        * (
            (2_f32 + 5_f32.sqrt()).powi(n as i32)
            - (2_f32 - 5_f32.sqrt()).powi(n as i32))
        ).round() as u64
}

4. 矩阵幂

类似斐波那契数列, 假设有

\[\begin{bmatrix} G(n)\\ G(n+1) \end{bmatrix} = \begin{bmatrix} a & b\\ c & d \end{bmatrix}^{n} \begin{bmatrix} G(0)\\ G(2) \end{bmatrix} \]

可以解得

\[\begin{bmatrix} G(n)\\ G(n+1) \end{bmatrix} = \begin{bmatrix} 0 & 1\\ 1 & 4 \end{bmatrix}^{n} \begin{bmatrix} 0\\ 2 \end{bmatrix} \]

也可作

\[\begin{bmatrix} G(n-1) & G(n)\\ G(n) & G(n+1)\\ \end{bmatrix} =2 \begin{bmatrix} 0 & 1\\ 1 & 4 \end{bmatrix}^{n} \]

故也可用矩阵幂或矩阵快速幂求解

之所以这个矩阵求幂能得到数列, 因为该矩阵同样能够满足递推公式的特征方程:

\[R^2 - 4R - E = O \]

斐波那契数列的矩阵同理:

\[R^2 -R - E = O \]

use std::ops::Mul;

#[derive(Copy, Debug, Clone)]
struct Mat2x2(u64, u64, u64, u64);
impl Mul for Mat2x2 {
    type Output = Self;
    fn mul(self, rhs: Self) -> Self::Output {
        Self(
            self.0 * rhs.0 + self.1 * rhs.2,
            self.0 * rhs.1 + self.1 * rhs.3,
            self.2 * rhs.0 + self.3 * rhs.2,
            self.2 * rhs.1 + self.3 * rhs.3,
        )
    }
}
impl Mul<u64> for Mat2x2 {
    type Output = Self;
    fn mul(self, rhs: u64) -> Self::Output {
        Self(self.0 * rhs, self.1 * rhs, self.2 * rhs, self.3 * rhs)
    }
}
impl Mat2x2 {
    fn pow(mut self, mut exp: u64) -> Self {
        let mut ans = Self(1, 0, 0, 1);
        while exp > 0 {
            if exp & 1 == 1 {
                ans = ans * self;
            }
            self = self * self;
            exp >>= 1;
        }
        ans
    }
}
fn even_fib(n: u64) -> u64 {
    let ori = Mat2x2(0, 1, 1, 4);
    (ori.pow(n) * 2).1
}

5. 快速倍增法

利用矩阵幂公式

\[\begin{bmatrix} G(n-1) & G(n)\\ G(n) & G(n+1)\\ \end{bmatrix} =2 \begin{bmatrix} 0 & 1\\ 1 & 4 \end{bmatrix}^{n} \]

可得

\[\begin{align*} \begin{bmatrix} G(2n-1) & G(2n)\\ G(2n) & G(2n+1)\\ \end{bmatrix} &=2 \begin{bmatrix} 0 & 1\\ 1 & 4 \end{bmatrix}^{2n}\\ &=2 \begin{bmatrix} 0 & 1\\ 1 & 4 \end{bmatrix}^{n}\begin{bmatrix} 0 & 1\\ 1 & 4 \end{bmatrix}^{n}\\ &=2 \begin{bmatrix} G(n-1) & G(n)\\ G(n) & G(n+1) \end{bmatrix}\begin{bmatrix} G(n-1) & G(n)\\ G(n) & G(n+1) \end{bmatrix}\\ \end{align*} \]

按矩阵乘法可得两个不同形式的公式:

\[\left\{\begin{align*} G(0)&=0&\\ G(1)&=2&\\ G(2n)&=2G(n)[G(n)+\frac{G(n-1)}{2}]&n>1\\ G(2n+1)&=\frac{G(n)^2}{2}+\frac{G(n+1)^2}{2} &n>1 \end{align*}\right. \]

\[\left\{\begin{align*} G(0)&=0&\\ G(1)&=2&\\ G(2n)&=2G(n)[\frac{G(n+1)}{2}-G(n)]&n>1\\ G(2n+1)&=\frac{G(n)^2}{2}+\frac{G(n+1)^2}{2} &n>1 \end{align*}\right. \]

利用以上公式得到以下算法:

fn even_fib_pair(n: u64) -> (u64, u64) {
    if n == 0 {
        (0, 2)
    } else {
        let (mut l, mut h) = (0u64, 2u64);
        let bit_length = n.ilog2();
        let mut mask = 1 << bit_length;
        while mask != 0 {
            (l, h) = (2 * l * (h/2 - l), l * l/2 + h * h/2);
            if n & mask != 0 {
                (l, h) = (h, l + 4 * h);
            }
            mask >>= 1;
        }
        (l, h)
    }
}

fn even_fib(n: u64) -> u64 {
    even_fib_pair(n).0
}

求和公式

斐波那契数列和偶斐波那契数列均存在简单的求和公式, 故欧拉计划第2题实际上可以用快速幂或快速倍增法以 \(O(\log n)\) 求特定项, 随后用求和公式以 \(O(1)\) 得到求和结果. 下面给出偶斐波那契数列求和公式证明方法, 斐波那契数列求和公式更简单, 请读者自行推导.

令

\[S(n) = \sum_{i=0}^{n}G(i) \]

有

\[\begin{align*} S(n) &= G(0) + G(1) + G(2) + G(3) + G(4) + \dots + G(n)\\ &= \textcolor{red}{4G(1)} + G(0) + G(1) + G(2) + G(3) + G(4) + \dots + G(n) \textcolor{blue}{- 4G(1)}\\ &= (\textcolor{red}{4G(1)} + G(0)) + G(1) + G(2) + G(3) + G(4) + \dots + G(n) \textcolor{blue}{- 4G(1)}\\ &= G(2) + G(1) + G(2) + G(3) + G(4) + \dots + G(n) \textcolor{blue}{- 4G(1)}\\ &= \textcolor{red}{4G(2)} + G(1) + G(2) + G(3) + G(4) + \dots + G(n) \textcolor{blue}{- 4G(1) - 3G(2)}\\ &= G(3) + G(2) + G(3) + G(4) + \dots + G(n) \textcolor{blue}{- 4G(1) - 3G(2)}\\ &= \textcolor{red}{4G(3)} + G(2) + G(3) + G(4) + \dots + G(n) \textcolor{blue}{- G(1) - 3(G(1) + G(2) + G(3))}\\ &= G(4) + G(3) + G(4) + \dots + G(n) \textcolor{blue}{- G(1) - 3(G(1) + G(2) + G(3))}\\ &\dots\\ &=\textcolor{red}{4G(n)}+G(n-1)+G(n) \textcolor{blue}{- G(1) - 3(G(1) + G(2) + G(3) + \dots + G(n-1) + G(n))}\\ &=G(n+1)+G(n) \textcolor{blue}{- G(1) - 3S(n)} \end{align*} \]

解得

\[S(n) =\frac{G(n+1) + G(n)-G(1)}{4} \]

利用求和公式和前文的求数列方法, 可以最快以 \(O(\log n)\) 求0至n项和, 以下为快速幂和快速倍增法的示例代码, 注意其中 even_fib 函数改为even_fib_pair 函数, 返回值从 \(G(n)\) 变成了 \((G(n+1), G(n))\).

1. 快速幂

use std::ops::Mul;

fn s(n: u64) -> u64 {
    let (gn_1, gn) = even_fib_pair(n);
    (gn_1 + gn - 2) / 4
}

#[derive(Copy, Debug, Clone)]
struct Mat2x2(u64, u64, u64, u64);
impl Mul for Mat2x2 {
    type Output = Self;
    fn mul(self, rhs: Self) -> Self::Output {
        Self(
            self.0 * rhs.0 + self.1 * rhs.2,
            self.0 * rhs.1 + self.1 * rhs.3,
            self.2 * rhs.0 + self.3 * rhs.2,
            self.2 * rhs.1 + self.3 * rhs.3,
        )
    }
}

impl Mul<u64> for Mat2x2 {
    type Output = Self;
    fn mul(self, rhs: u64) -> Self::Output {
        Self(self.0 * rhs, self.1 * rhs, self.2 * rhs, self.3 * rhs)
    }
}

impl Mat2x2 {
    fn pow(mut self, mut exp: u64) -> Self {
        let mut ans = Self(1, 0, 0, 1);
        while exp > 0 {
            if exp & 1 == 1 {
                ans = ans * self;
            }
            self = self * self;
            exp >>= 1;
        }
        ans
    }
}

fn even_fib_pair(n: u64) -> (u64, u64) {
    let ori = Mat2x2(0, 1, 1, 4);
    let mat = ori.pow(n) * 2;
    (mat.3, mat.1)
}

2. 快速倍增法

fn s(n: u64) -> u64 {
    let (gn_1, gn) = even_fib_pair(n);
    (gn_1 + gn - 2) / 4
}

fn even_fib_pair(n: u64) -> (u64, u64) {
    if n == 0 {
        (2, 0)
    } else {
        let (mut h, mut l) = (2u64, 0u64);
        let bit_length = n.ilog2();
        let mut mask = 1 << bit_length;
        while mask != 0 {
            (l, h) = (2 * l * (h / 2 - l), l * l / 2 + h * h / 2);
            if n & mask != 0 {
                (l, h) = (h, l + 4 * h);
            }
            mask >>= 1;
        }
        (h, l)
    }
}

虽然这两种方法都可以由 \(O(\log n)\) 求和, 但是无法提前预知哪一项开始大于 4M, 如果仅针对欧拉计划第2题, 采用偶斐波那契数列的动态规划法, 每项判断其是否小于 4M 也很高效, 从数学角度看最好的方法则是用公式法计算

\[\frac{1}{\sqrt{5}}[(2 + \sqrt{5})^n - (2 - \sqrt{5})^n] \approx \frac{1}{\sqrt{5}}(2 + \sqrt{5})^n \leq 4,000,000 \]

得到

\[n \leq \frac{\ln(4000000\sqrt{5})}{\ln(2 + \sqrt{5})} \approx 11.08764 \]

即可得

\[S(11)=4613732 \]