2020 - 04 - 29 17:39:49

斐波那契之旅

很多算法问题最终都可以归结到斐波那契数列的问题上。

譬如下面这个问题

一个数只用 1,2,3,4这四种数字组成，而且相近的两位数永远相差1，比如1234321, 1212121等。这样的n位数有多少个？

备忘录，尾递归

递归方法的时间复杂度过于惨烈，有两个简化的思路。

其一是备忘录方法，以一个数组避免重复计算，把能剪掉的枝都剪秃了，把递归树剪成了一根竿（？）

res= [1,1] + [0 for _ in range(100)]
def fib(n):
    if not res[n-1]:
        res[n-1] = fib(n-1)+fib(n-2)
    return res[n-1]

另一种方法是尾递归。

没什么意思，就是个循环，但是优雅。

细究的话有滑动窗口的思想。

def fib(n,a=0,b=1):
    if n!=1:
        return fib(n-1,b,a+b)
    return b

注：python这厮本身没有尾递归优化，这种写法虽然避免了直接递归的搜索树，但是仍然在递归过深时会爆栈

数论方法

斐波那契数列是有通项公式的，但是一般并不会有人用这个公式编程求解斐波那契数列。 \[ F_n=\frac{\sqrt{5}}{5}[{(\frac{1+\sqrt{5}}{2})}^n-{(\frac{1-\sqrt{5}}{2})}^n] \]

这是因为$\sqrt5$ 是一个烦人的东西，计算机不能存储参与整形或浮点型运算的无理数

有没有可以绕过$\sqrt5$ 使用通项公式的方法呢？

$\sqrt5$ 的性质在于，$(\sqrt5)^2=5$ ，这好像什么都没有说。因为这是它的定义。

若一个非负数x的平方等于a，即x²=a，则这个数x叫做a的算术平方根

如果在数论中考虑，希望找到整数x，使$x^2=5\ mod\ n$ 。那么在对n取模的时候，也许x可以代替$\sqrt5$

这里为了方便起见，取n为同时满足5^(n-1)/2 ≡ 1 (mod n) 以及n ≡ 3 (mod 4) 的素数，记为p。

复习一下二次剩余的性质

二次剩余
若p是奇质数且p不能整除d, 则方程 x² ≡ d (mod p)有解的充要条件是:
d^(p-1)/2 ≡ 1 (mod p) (欧拉判别准则)
特殊的, 若p ≡ 3 (mod 4)
则可进一步解出来 x ≡ d^(p+1)/4(mod p)
　　事实上这是一个极为精妙的构造解:
　　d^(p-1)/2 ≡ 1 (mod p) ,两边同乘d得:
　　d^(p+1)/2≡ d (mod p)
　　因为 p ≡ 3 (mod 4), 故(p+1)/4是一个整数
　　　　(d^(p+1)/4)² ≡ d (mod p)
　　　　观察得: x ≡ ±d^(p+1)/4 (mod p)

可以快速得出 $x ≡ ± 5^{(p+1)/4} \ mod\ p$

在mod p的情况下，x可以取代 $\sqrt5$

以p = 6655503139 为例，

根据二次剩余的知识求出来 x₁ = 4865890845，x₂=1789612294。

考虑到通项公式里有$\frac{1+\sqrt5}{2}$ 和 $\frac{1-\sqrt5}{2}$ 的形式，我们希望选取x为奇数。取 x = 4865890845。

再将$\sqrt{5}/5$ 写成 $1/\sqrt{5}$的形式，利用广义欧几里得算法求出 x 的逆元即可。 \[ F_n=973178169*(2432945423^n-4222557717^n) \mod 6655503139 \] 以下是随机生成的另外几组斐波那契数列在对质数取模时的通项公式

\[ \begin{align*} F_n&=2274850784*(1949587701^n-5525490819^n) \mod 7475078519\\ F_n&=797768463*(1994421158^n-5676839322^n) \ \ \mod 7671260479\\ F_n&=1062316053*(2655790133^n-4350386407^n) \mod 7006176539\\ F_n&=2101662618*(1392974856^n-6329388524^n) \mod 7722363379\\ F_n&=685265783*(1713164458^n-3687114902^n)\ \ \mod 5400279359\\ F_n&=1930655058*(895691850^n-6966199742^n)\ \ \mod 7861891591\\ F_n&=2950851583*(1758659027^n-3859810905^n) \mod 5618469931\\ F_n&=2324981816*(2844389075^n-3091741857^n) \mod 5936130931\\ F_n&=4172362562*(1853481697^n-3864801443^n) \mod 5718283139\\ F_n&=5096301964*(1800360282^n-5493236138^n) \mod 7293596419\\ \end{align*} \]

放代码，以下是随机选取p并生成满足要求的通项公式的方法

import random

def miller(n, a, m, t):
    b = pow(a, m, n)
    if b == 1 or b == (n - 1):
        return True
    for i in range(t - 1):
        b = pow(b, 2, n)
        if b == (n - 1):
            return True
    return False

def MR(p, k): #Rabin-Miller素性检测
    pAss = False
    block = []
    for i in range(k):
        block.append(random.randint(2, p - 2))
    m = p - 1
    t = 0
    while not m % 2:
        m = m // 2
        t += 1
    for a in block:
        pAss = miller(p, a, m, t)
        if not pAss:
            return False
    return True

def prime(A,B):
    while True:
        p = random.randrange(A,B,4)#A必须为奇数
        if MR(p, 9):
            break
    return p

def egcd(a, b): #扩展欧几里得算法
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)  # g为公因子
        return (g, x - (b // a) * y, y)

def modinv(a, m): #求乘法逆元
    g, x, y = egcd(a, m)
    if (g != 1):
        print('modular inverse does not exist')
    else:
        return x % m
        
def prime2(): #生成符合要求的素数
    while(True):
        p=prime(4*10**9+3,8*10**9+3)
        if (1==pow(5,(p-1)//2,p)): #欧拉判别准则
            return p

p=prime2()
x=pow(5,(p+1)//4,p)    #二次剩余求解x
x=x if x&1 else p-x #x 取奇数

#def fib(n):
#    return (modinv(x,p))*(pow(((1+x)//2),n,p)-pow(((2*p+1-x)//2),n,p))%p

A=modinv(x,p)
B=(1+x)//2
C=(2*p+1-x)//2
#print("Fn={0}*({1}^n-{2}^n) mod {3}".format(A,B,C,p))
print("{0}*(pow({1},n,{3})-pow({2},n,{3})))%{3}".format(A,B,C,p))

Cipolla算法

上面的算法只涵盖了形如4x+3的素数

对于另一部分素数，若存在5的二次剩余，则可用Cipolla算法在$O(log^2n)$的时间复杂度内进行求解

参考链接1

参考链接2

中国剩余定理

一个想法，还没思考过具体应用，当F_n大于模数时，可以考虑用多个模数分别求解，再用中国剩余定理将各自结果进行拼合以求出最终结果

矩阵快速幂

把数列的递推公式写成矩阵形式：

\[ \begin{gathered} \begin{pmatrix} F_{n} \\ F_{n-1} \end{pmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} = \begin{pmatrix} F_{n} + F_{n-1} \\ F_{n} \end{pmatrix}= \begin{pmatrix} F_{n+1} \\ F_{n} \end{pmatrix} \end{gathered} \] 而后可以倒着地把递归写成累乘的形式 \[ \begin{gathered} \begin{pmatrix} F_{n} \\ F_{n-1} \end{pmatrix}= \begin{pmatrix} F_{2} \\ F_{1} \end{pmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2} \end{gathered} \] 注意到F₁和F₂都是1，而且可以定义F₀：F₂ = F₁ + F₀，从而得到F₀ = F₂ - F₁ = 0

把向量拼合成矩阵以进行合并 \[ \begin{align*} \begin{bmatrix} F_{n} & F_{n-1} \\ F_{n-1} & F_{n-2} \end{bmatrix} &= \begin{bmatrix} F_{2} & F_{1} \\ F_{1} & F_{0} \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2} \\&= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2} \\&= \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-1} \end{align*} \]

利用矩阵快速幂，可以做到$O(logn)$级别的求解，很优雅

由对称性，每个矩阵可以少存一个数字

#include<iostream>
#define uint64 unsigned long long
using namespace std;

const uint64 MOD=1000000007;  //MOD不应超出uint32 

struct Matrix{
    uint64 A11,A12,A22;
    Matrix operator * (const struct Matrix &B) const{
        Matrix ans;
        ans.A11=(A11*B.A11%MOD + A12*B.A12%MOD)%MOD;
        ans.A12=(A11*B.A12%MOD + A12*B.A22%MOD)%MOD;
        ans.A22=(A12*B.A12%MOD + A22*B.A22%MOD)%MOD;
        return ans;
    }
};

uint64 Fibonacci(int k){
    k--;
    Matrix ans={1,0,1};
    Matrix p = {1,1,0};
    while(k){      //矩阵快速幂
        if (k & 1) ans = ans * p;
        p = p * p;
        k >>= 1;
    }
    return ans.A11;
}

注：c++的结构体不取别名可直接声明此类型变量

二分递归方法

出处：斐波那契数，计算与分析

我们注意到如下规律： \[ \begin{align*} F_n&=F_{n-1}+F_{n-2}\\ &=1F_{n-1}+1F_{n-2}=F_2F_{n-1}+F_1F_{n-2}\\ &=2F_{n-2}+1F_{n-3}=F_3F_{n-2}+F_2F_{n-3}\\ &=3F_{n-3}+2F_{n-4}=F_4F_{n-3}+F_3F_{n-4}\\ &=5F_{n-4}+3F_{n-5}=F_5F_{n-4}+F_4F_{n-5}\\ &\quad\ ............\\ &=F_kF_{n-k+1}+F_{k-1}F_{n-k} \end{align*} \] 对F_n分奇偶性进行讨论
当n=2m+1：

$F_n=F_{m+1}^2+F_{m}^2$

当n=2m：

$F_n=F_{m}(F_{m+1}+F_{m-1})=F_{m}(F_{m}+2F_{m-1})$

未进行记忆化搜索时，时间复杂度是$O(n)$

记忆化搜索后，时间复杂度可以达到$O(logn)$，

而且比同等级的矩阵快速幂算法具有更小的常数

如下是用Python的lru_cache自动实现记忆化的代码

from functools import lru_cache
@lru_cache(maxsize=128)
def fast_fib_lru_cache(n):
    if n < 2:
        return n
    x = fast_fib_lru_cache((n >> 1) - 1)
    y = fast_fib_lru_cache(n >> 1)
    if isodd(n):
        x += y
        return x * x + y * y
    else:
        return y * (y + 2 * x)

Cassini公式

$F_n^2=F_{n+1}F_{n-1}-(-1)^n$

将以上算法用Cassini公式优化后即是 GMP 大数运算库用来计算Fibonacci数的算法

$$ \[\begin{align*} F_{2k-1}&=F_{k}^2+F_{k-1}^2\\ F_{2k+1}&=4F_{k}^2-F_{k-1}^2+(-1)^k\\ F_{2k}&=F_{2k+1}-F_{2k-1}\\ &=3F_{k}^2-2F_{k-1}^2+(-1)^k\\ \end{align*}\] $$

(懒得搬了)

这里有运用Cassini公式的实现

另一种优化

注意到可以把递推式改写成如下形式

\[ F_{2m+1}=F_{2m}+F_{2m-1}\ \Rightarrow \ F_{2m}=F_{2m+1}-F_{2m-1}\\ \]

且在奇数时只需要算两次乘法

\[ \begin{aligned} F_{2m+1}=F_{m+1}^2+F_{m}^2\\ F_{2m-1}=F_{m}^2+F_{m-1}^2 \end{aligned} \]

利用上述恒等式可把偶数时的情况也简化为两次乘法

\[ \begin{align*} F_{2m}&=F_{2m+1}-F_{2m-1}\\ &=F_{m+1}^2+F_{m}^2-F_{m}^2-F_{m-1}^2\\ &=F_{m+1}^2-F_{m-1}^2 \end{align*} \]

主要是比较优雅，其实没什么用处

dp=[0,1,1,2]+[-1]*100
def fib(n):
	if dp[n]!=-1:
		return dp[n]
	if n&1:
		dp[n]=fib((n>>1)+1)**2 + fib((n-1)>>1)**2
	else:
		dp[n]=fib((n>>1)+1)**2 - fib((n-1)>>1)**2
	return dp[n]

不正经的方法

我们注意到在 不对结果取模 的情况下

64位有符号整型最多求解到 F₉₂

64位无符号整型最多求解到 F₉₃

32位有符号整型最多求解到 F₄₆

32位无符号整型最多求解到 F₄₇

所以一个最自然的思路是事先对这个数列进行计算并存入对应类型的数组中。在程序运行时直接访问数组取对应值，O(1)。

One More Thing

def fib(n):
    if n <= 2:
        return 1
    return  fib(n-1) + fib(n-2)

朴素递归版代码的时间复杂度是多少？

这是一棵未经任何剪枝的递归树。

鉴于没有复杂计算，递归调用次数即是其计算时间复杂度的关键逻辑，可将递归树节点规模视作时间复杂度。

那么递归调用了多少次？

观察代码本身，F_n 是由一个一个 1 加起来的，递归树的每一个叶子节点都会返回一个 1，所以递归树叶子节点的数量即是 F_n 。

观察递归结构，易知 递归树中没有度为1的节点。

经过简单推导，此种二叉树中 非叶子节点数量 = 叶子节点数量 - 1 。

故 递归树节点总数 为: 非叶子节点数量 + 叶子节点数量 = $2F_n-1$ 。

时间复杂度为 $O(F_n)$ 。

仍不够直观

作为数量级而言 F_n 不够直观，我们需要对 $O(F_n)$ 进行一定程度的化简。

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

首先扔掉常数，剩下 ${(\frac{1+\sqrt{5}}{2})}^n-{(\frac{1-\sqrt{5}}{2})}^n$

注意到 $\frac{1-\sqrt{5}}{2}≈-0.618<1$ , 故可以舍弃后一项

$O(F_n)=O((1+\phi)^n)=O((\frac{1+\sqrt{5}}{2})^n)≈O(1.618^n)$