斐波纳切数列相关算法竞赛知识(一):记忆化搜索、快速幂、矩阵快速幂

众所周知,斐波纳切数列是这样的一个数列,F(0)=0,F(1)=1,F(n)=F(n-1)+F(n-2)。
围绕如何算出F(n),有很多算法竞赛相关知识,总结如下。

一、递归

入门级思维,写一个函数,不断调用自己即可。优点是书写简单,缺点是效率极低,因为有很多重复计算,比如要计算f(6),需计算f(5)+f(4),先算f(5),中途计算一次f(4),再算f(4),又算了一次。

int f(int n){
    return n < 2 ? n : f(n-1) + f(n-2);
}

二、记忆化搜索

记忆化也是算法的一个基本技能,上面的做法既然会重复计算,那么第一次算出来f(i)的时候,把它保存(记忆)起来,下次需要的时候直接取出,加快速度。

int a[40];
int f(int n){
    if(n < 2) return n;
    if(a[n] > 0) return a[n];//已经计算过了,直接返回
    else return a[n] = f(n-1)+f(n-2);//没有算过,先算,然后保存起来,最后返回。
}

三、循环,动态规划启蒙

利用循环,直接根据前两项算出f(n),效率和记忆话搜索一样,去掉了递归。此处蕴含了初级的动态规划思想,用前面已作出的决策,决定当前要采取的决策,当前采取的决策,也不会影响后面决策的正确性。

int a[40] = {0, 1};
int f(int n){
    if (n < 2 ) return n;
    for(int i = 2; i <= n; i++){
        a[i] = a[i-1] + a[i-2];
    }
    return a[n];
}

四、整数范围

上面的程序只对n比较小的时候有效,如果计算f(47),就发生溢出,因为结果比int类型能表示的最大数还大。int是32位的,最大能表示2^31-1,大概是20多亿。可用unsigned long long容纳更大的数,unsigned long long 是64位,并且最高位也用来计数,最大能表示2^64-1,十进制是18446744073709551615,可以用程序输出这个数,体会一下:

#include <cstdio>
#include <climits>//常用数据类型的最大、最小值都在这个头文件里
int main(){
    printf("%llu", ULLONG_MAX);
    return 0;
}

2^64-1对于斐波纳切数列这样的指数级别函数的增长速度来说,也不算很大,大概f(200)就快要超出去了。

五、模

n再大的话怎么办,一般题目要求计算结果对一个很大的数M(比如1e8+7,100000007)求模,求模以后结果总会小于M。关于求模有两个重要公式需要记住:
(a+b)%c = (a %c + b%c)%c
(ab)%c = ((a%c) * (b%c))%c
加入模运算以后,n很大的时候也可以求了,把循环代码改一下即可,不需要用很大的数组记录所有0-n的结果,只需要2个元素就够了:

#include <cstdio>
#define MOD 1000000007
typedef long long ll;
ll a[2] = {0, 1};
ll f(int n){
    if (n < 2 ) return n;
    for(int i = 2; i <= n; i++){
        ll tmp = (a[0] + a[1]) % MOD;
        a[0] = a[1];
        a[1] = tmp;
    }
    return a[1];
}
int main(){
    printf("%lld", f(200000000));
    return 0;
}

六、矩阵乘法辅助转移

上面的代码当n在10亿范围内勉强能算出来,再大的话就很无力了。因为复杂度是O(n),只用简单的手段已经没有优化的空间了,要想加速必须用其它办法。办法就是利用矩阵乘法。矩阵乘法运算规则此处就不展开了,学习资源随处可见。
矩阵象鲁迅先生一样装作很无辜:“你算斐波纳切数列,和我矩阵有什么关系???”
我们先不管矩阵先生的抗议,试定义数列初始值是一个1*2(1行,2列)的矩阵,表示[f(0),f(1)],如果这个矩阵乘以另外一个矩阵,结果能变成[f(1),f(2)],再乘一次,变成[f(2),f(3)],显然乘上n-1次,结果就变成了[f(n-1),f(n)]。把这个式子一般化,即[f(i-2),f(i-1)] xA = [f(i-1),f(i)],根据矩阵乘法定义,A显然是一个2x2的矩阵,第一列是0,1;第二列是1,1,乘完以后,第一列是f(i-2)x0+f(i-1)x1,即f(i-1);第二列是f(i-2)x1+f(i-1)x1,正好是f(i)。
结论是:[f(n-1),f(n)] = [f(0),f(1)] x A^(n-1),[f(n-1),f(n)] 等于[f(0),f(1)] 乘以A的n-1次幂。
矩阵再次装无辜:“这样也没用,算这么多次乘法,还不如原来算加法快呢,你们还是放了我找别人吧......”

七、快速幂

矩阵n次幂,如果能快速算出来,我们的目的就达到了。幸运的是办法很简单,也许对于矩阵来说这是不幸的^ _ ^。
先看看整数的幂怎么算。普通方法用循环乘,当然很慢数据上了1亿就不能秒出了。用初中数学知识把做法稍微改一下即可:
假设n是偶数,A^n = (AA)^(n/2)
若n是奇数,A^n = AxA(n-1),n-1肯定是偶数,所以A^n = Ax(AA)^((n-1)/2)
根据上述原理,一次乘法运算将底数变成原来的平方,就能把幂减少1半,效率非常高,再大的数,也经不起几次减半啊,时间复杂度从O(n)降到O(logn)。
代码也很简单,请运行并体验一下它的威力:

#include <cstdio>
typedef long long ll;
ll fastpow(ll base, ll p, int k){//快速幂,结果对k取模
    ll ret = 1;
    while(p > 0){
        if (p & 1){//奇数,相当于p % 2 == 1
            ret = ret * base % k;
        }
        base = base * base % k;//底数变成原来的平方
        p >>= 1;//相当于p /= 2,指数减半
    }
    return ret % k;
}
int main(){
    printf("%lld", fastpow(3, 8888888880000088, 100000007));
    return 0;
}

八、矩阵快速幂

这下矩阵先生再也躲不掉了,乖乖地干活吧。快速幂同样适用于矩阵,因为行和列大小相等的矩阵(方阵)乘法满足结合率。代码如下:

#include <cstdio>
#include <cstring>
#define MOD 1000000007
using namespace std;
typedef long long ll;
/*
从[f(i-2), f(i-1)]转移到[f(i-1), f(i)]
这是一个1x2的矩阵,乘一个2x2的加速矩阵,就变成一个新的1x2矩阵
加速矩阵是
0 1
1 1
*/

ll n;
struct Mx{
    ll v[2][2];
    int n, m;//行,列
    Mx(int rows, int cols) : n(rows), m(cols){//构造函数
        memset(v, 0, sizeof(v));
    }
    Mx operator*(const Mx &r){//重载乘法运算符,为了代码书写方便
        Mx ret(n, r.m);
        for(int i = 0; i < n; i++){
            for(int j = 0; j < r.m; j++){
                for(int k = 0; k < m; k++){//此处注意k的循环放在最内层,有助于cpu缓存命中,加快运算速度,放在最外层结果正确,但速度慢。
                    ret.v[i][j] = (ret.v[i][j] + (v[i][k] * r.v[k][j] % MOD)) % MOD;
                }
            }
        }
        return ret;
    }
} ans(1, 2), rush(2, 2);
ll calc(ll t){
    if (t == 1 || t == 2){
        return 1;
    }
    t -= 3;//矩阵乘t次,base的t次方
    ans.v[0][0] = ans.v[0][1] = 1;//将结果矩阵初始化为[1, 1],n从1开始算
    rush.v[0][1] = rush.v[1][0] = rush.v[1][1] = 1;//初始化加速矩阵,第一列0,1;第二列1,1
    Mx trans = rush;
    while(t > 0){
        if(t & 1){
            trans = trans * rush;
        }
        rush = rush * rush;
        t >>= 1;
    }
    ans = ans * trans;
    return ans.v[0][1];
}
int main(){
    scanf("%lld", &n);
    printf("%lld\n", calc(n));
    return 0;
}

九、斐波纳切公约数

gcd(f(a),f(b)) = f(gcd(a,b))
gcd表示最大公约数,f(n)表示斐波纳切数列第n位。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,457评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,837评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,696评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,183评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,057评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,105评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,520评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,211评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,482评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,574评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,353评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,213评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,576评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,897评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,174评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,489评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,683评论 2 335

推荐阅读更多精彩内容