斐波拉契

Posted 2020-11-08 shine-sky

在OI简单数论中斐波拉契是常常出现的东西

是什么
斐波那契数列，又称黄金分割数列、因数学家列昂纳多·斐波那契以兔子繁殖为例子而引入，故又称为“兔子数列”，指的是这样一个数列：1、1、2、3、5、8、13、21、34、……在数学上，斐波纳契数列以如下被以递归的方法定义：F(0)=1，F(1)=1, F(n)=F(n-1)+F(n-2)（n>=2，n∈N*）在现代物理、准晶体结构、化学等领域，斐波纳契数列都有直接的应用.

 

(百度百科)

模板
1,递归

1 inline int Fib(int n)
2 {
3     if (n == 1 || n == 2)
4     {
5         return 1;
6     }
7     return Fib(n - 1) + Fib(n - 2);
8 }

View Code

 

复杂度分析
由于对于每一个1都是最后一层递归返回上来的故会递归F(n)次 , 由于斐波拉契数列是随着指数上升的 故复杂度约为O(2^n)

2,循环

1 inline int Fib(int n)
2 {
3     F[1] = F[2] = 1;
4     for(int i = 3 ; i <= n ; ++ i)
5         F[i] = F[i - 1] + F[i - 2];
6     return F[n];
7 }

View Code

 

复杂度分析 :O(n)

 

3,矩阵乘法优化
讲数列放在矩阵乘法中会发现

$$\left[\begin{matrix}F[n - 1]\\F[n - 2]\end{matrix}\right] * \left[\begin{matrix}1 & 1\\1 & 0\end{matrix}\right] = \left[\begin{matrix}F[n]\\F[n - 1]\end{matrix}\right]$$

即如果求F[n]
$$\left[\begin{matrix}F[n]\\F[n - 1]\end{matrix}\right] = \left[\begin{matrix}F[n - 1]\\F[n - 2]\end{matrix}\right] * \left[\begin{matrix}1 & 1\\1 & 0\end{matrix}\right] = \left[\begin{matrix}F[n - 2]\\F[n - 3]\end{matrix}\right] * \left[\begin{matrix}1 & 1\\1 & 0\end{matrix}\right] ^ 2 = .... = \left[\begin{matrix}F[1] \\ F[2]\end{matrix}\right] * \left[\begin{matrix}1 & 1\\1 & 0\end{matrix}\right]^{n-1}$$

将$$\left[\begin{matrix}F[1] \\ F[2]\end{matrix}\right]$$
拓展到2 * 2即
$$\left[\begin{matrix}1 & 0 \\ 0 & 1\end{matrix}\right]$$

 1 #include <iostream>
 2 #include <cstdio>
 3 #include <cstring>
 4 #define inc(i) (++ i)
 5 #define int long long
 6 using namespace std;
 7 const int Mod = 1000000000 + 7;
 8 struct JX
 9 {
10     int a[3][3];
11     friend JX operator * (JX a , JX b)
12     {
13         JX c;
14         memset(c.a , 0 , sizeof(c.a));
15         for(int i = 1 ; i <= 2 ; inc(i))
16             for(int j = 1 ; j <= 2 ; inc(j))
17                 for(int k = 1 ; k <= 2 ; inc(k))
18                     c.a[i][j] = (c.a[i][j] + a.a[i][k] * b.a[k][j]) % Mod;
19         return c;
20     }
21 }Ans , A;
22 int n;
23 inline void KSM()
24 {
25     int p = n - 2;
26     while(p)
27     {
28         if(p & 1) Ans = Ans * A;
29         p >>= 1 , A = A * A;
30     }
31 }
32 signed main()
33 {
34     scanf("%ld" , &n);
35     if(n <= 2)  putchar(49);
36     else
37     {
38         Ans.a[1][1] = A.a[1][1] = 1;//直接跳到n = 2时的矩阵
39         Ans.a[1][2] = A.a[1][2] = 1;
40         Ans.a[2][1] = A.a[2][1] = 1;
41         KSM();
42         printf("%lld" , Ans.a[1][1]);
43     }
44     return 0;
45 }

View Code

 

 

4 , 算不上什么优化的优化 2(斐波拉契数列膜p的循环节)
转自[Fib数模n的循环节](https://blog.csdn.net/acdreamers/article/details/10983813)(代码原创)

我们知道Fibonacci数列，现在我们来求一个Fib数模n的循环节的长度.

对于一个正整数n，我们求Fib数模n的循环节的长度的方法如下：

    （1）把n素因子分解，即
    （2）分别计算Fib数模每个的循环节长度，假设长度分别是
    （3）那么Fib模n的循环节长度

从上面三个步骤看来，貌似最困难的是第二步，那么我们如何求Fib模的循环节长度呢？
 
这里有一个优美的定理:
>Fib数模的最小循环节长度等于，其中表示Fib数模素数的最小循环节长度。可以看出我们现在最重要的就是求
 
 
对于求我们利用如下定理：
 
   如果5是模的二次剩余，那么循环节的的长度是的因子，否则，循环节的长度是的因子。

顺便说一句，对于小于等于5的素数，我们直接特殊判断，loop(2)=3,loop(3)=8,loop(5)=20。

那么我们可以先求出所有的因子，然后用矩阵快速幂来一个一个判断，这样时间复杂度不会很大。
 

模板代码：(巨丑 见谅)

1#include<iostream>2#include<cstdio>3#include<cctype>4#include<cstring>5#include<cmath>6#include<algorithm>7#defineintlonglong8#defineinc(i)(++i)9#definedec(i)(--i)10usingnamespacestd;11constintN=100000+7;12intn,Prime[N+7],tot,Pri[N],Cnt[N],cnt,TOT;13intFactor[N],P;14boolNo_Prime[N+7];15structJX16{17inta[4][4];18friendJXoperator*(JXa,JXb)19{20JXc;21memset(c.a,0,sizeof(c.a));22for(inti=1;i<=2;inc(i))23for(intj=1;j<=2;inc(j))24for(intk=1;k<=2;inc(k))25c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%P)%P;26returnc;27}28}A,X;29inlineJXKSM_2(JXa,intq,intP)30{//矩阵优化斐波拉契数列31JXAns=X,x=a;32while(q)33{34if(q&1)Ans=Ans*x;35q>>=1,x=x*x;36}37returnAns;38}39inlinevoidGet_Prime()40{//得出质数便于之后分解质因数41No_Prime[1]=1;42for(inti=2;i<=N;inc(i))43{44if(!No_Prime[i])Prime[inc(tot)]=i;45for(intj=1;j<=tot&&i*Prime[j]<=N;inc(j))46{47No_Prime[i*Prime[j]]=1;48if(i%Prime[j]==0)break;49}50}51}52inlinevoidResolve(intn,intPri[],intCnt[])53{//分解质因数,pri存质因子,cnt存次数54cnt=1;55intt=sqrt(n),TOT;56for(inti=1;Prime[i]<=t;inc(i))57if(n%Prime[i]==0)58{59TOT=0;60Pri[cnt]=Prime[i];61while(n%Prime[i]==0)62{63n/=Prime[i];64inc(TOT);65}66Cnt[cnt]=TOT;67inc(cnt);68}69if(n^1)70{71Pri[cnt]=n,Cnt[cnt]=1;72inc(cnt);73}74dec(cnt);75}7677inlineintKSM_1(intx,intq,intP)78{//普通的快速幂79intAns=1;80x%=P;81while(q)82{83if(q&1)Ans=Ans*x%P;84q>>=1,x=x*x%P;85}86returnAns;87}88inlinevoidWork(intn)89{//求n的因数90TOT=0;91intt=sqrt(n);92for(inti=1;i<=t;inc(i))93{94if(n%i==0)95{96if(i*i==n)Factor[inc(TOT)]=i;97else98{99Factor[inc(TOT)]=i;100Factor[inc(TOT)]=n/i;101}102}103}104}105inlineintGcd(inta,intb)106{107returna%b==0?b:Gcd(b,a%b);108}109inlineintFind_Loop(intn)110{111Resolve(n,Pri,Cnt);112intAns=1,loop;113for(inti=1;i<=cnt;inc(i))114{115loop=1;116P=Pri[i];117switch(P)118{119case2:loop=3;break;120case3:loop=8;break;121case5:loop=20;break;122default:123{124if(KSM_1(5,(P-1)>>1,P)==1)125Work(P-1);126else127Work(2*(P+1));128sort(Factor,Factor+TOT+1);129for(intj=1;j<=TOT;inc(j))130{131JXB=KSM_2(A,Factor[j]-1,P);132intx=(B.a[1][1]+B.a[1][2])%P;133inty=(B.a[2][1]+B.a[2][2])%P;134if(x==1&&y==0)135{136loop=Factor[j];137break;138}139}140}break;141}142loop*=KSM_1(P,Cnt[i]-1,9223372036854775807);143if(loop<Ans)swap(loop,Ans);144Ans=Ans/Gcd(loop,Ans)*loop;145}146returnAns;147}148inlineintMod(char*a,intb)//高精度a除以低精度b149{150intd=0,l=strlen(a);151for(inti=0;i<l;inc(i))d=(d*10+(a[i]-‘0‘))%b;//求出余数152returnd;153}154charpp[30000007];155signedmain()156{157A.a[1][1]=A.a[1][2]=A.a[2][1]=1;158X.a[1][1]=X.a[2][2]=1;159Get_Prime();160scanf("%s%lld",pp,&n);161intloop=Find_Loop(n);162intAns=Mod(pp,loop);163if(!Ans)Ans=loop;164if(Ans<=0)putchar(48);165else166if(Ans<=2)printf("%lld",1ll%n);167else168{169A.a[1][1]=A.a[1][2]=A.a[2][1]=1;170X.a[1][1]=X.a[2][2]=1;171P=n;172X=KSM_2(A,Ans-1,n);173printf("%lld",X.a[1][1]);174}175return0;176}

Ugly Code