今天听人民日报说算fibonacci数列的是顶尖科学家，我也来算一算。
int fib(ll n){
ll x=0,y=1,t;
for (int d=62-__builtin_clz(++n);d>=0;d--){
t=(x*x+y*y)%P; y=(x+x+y)*y%P; x=t;
if (n>>d&1) x=y,y=(t+y)%P;
}return x;
}
ccsp一道题比赛卡常数算fibonacci，非常难受。
矩阵快速幂谁都会，慢得一匹。当然再暴力三重for循环的时候也记要把阶数写个常数2，不然-O2编译器没法循环展开。一个有效的观察是，发现这个矩阵的幂反对角线上
然后我们知道，计算这种根式$$\frac{1}{\sqrt{5}} ((\frac{\sqrt{5}+1}{2})^n+ (\frac{\sqrt{5}-1}{2})^n )$$的时候，是可以对根式做快速幂的，非常快，比矩阵快5~10倍 (模数固定的话，可以先尝试一下二次剩余，能的话更快) 。美中不足的是，计算这个需要逆元，对模数有要求。
然后就看看怎么得到开头的代码。
几个恒等式：
F[n+m]=F[n-1]F[m]+F[n]F[m+1]
F[2n+1]=[F[n]]^2+[F[n+1]]^2
F[2n]=F[n+1]^2-F[n-1]^2=(2*F[n-1]+F[n])*F[n]
第一条可由数学归纳法得到，然后可推出后面两条。
有了后面两条，可直接递归，set存之前算出来的值记忆化，写起来比矩阵简短一点。但是set比较慢，时间复杂度和实际速度可能不如矩阵。
注意到上述计算有一个特点，就是相邻两个数n,n+1可推出2n,2n+1的值。因此递归的时候每次除以2，就可以只用保留两个数，不需要set。分奇偶讨论。
pair<int,int> f(int n){
if (n==0) {return {0,1};}
if (n==1) {cout<<(n&1); return {1,1};}
auto r=f(n-1>>1);
int a=(1ll*r.first*r.first+1ll*r.second*r.second)%P;
int b=(2ll*r.first+r.second)*r.second%P;
cout<<(n&1); //观察什么时候为奇数
if (n&1) return {a,b};
else return {b,(a+b)%P};
}
但是这个递归写起来比较长，看着不爽。想办法再优化一下，改成直接递推。
注意到递归按n分奇偶讨论。递推如果不能直接知道各个位的奇偶情况，那还是要栈，还不如写成递归。还好这里每次是n-1>>1，比较简单，可以找规律。n-1>>1相当于偶数的时候-1再右移一位，奇数则直接右移。因此连续的0会被取反成1，它之前的一个1会变成0，然后连锁再导致前面的连续1变0，所以效果总结起来就是取反。不放心就再打表测一下~(n+1)与上面是不是相等。
int f1(int n){
int x=0,y=1,d=0,n1=n;
while (n1>>=1) d++;
n1=~(n+1);
if ((n1>>d&1)) x=1;
d--;
for (;d>=0;d--){
cout<<(n1>>d&1);
int tx=(1ll*x*x+1ll*y*y)%P;
int ty=(2ll*x+y)*y%P;
if ((n1>>d)&1) x=tx,y=ty;
else x=ty,y=(tx+ty)%P;
}return x;
}
很好，对上了。最后工作就是对f1化简一下，就是上面的代码了。计算量和根式快速幂差不多，优点是不限制模数。
完整的测试代码：
#include<bits/stdc++.h>using namespace std;
const int P=1000000007;
pair<int,int> f(int n){
if (n==0) {return {0,1};}
if (n==1) {cout<<(n&1); return {1,1};}
auto r=f(n-1>>1);
int a=(1ll*r.first*r.first+1ll*r.second*r.second)%P;
int b=(2ll*r.first+r.second)*r.second%P;
cout<<(n&1);
if (n&1) return {a,b};
else return {b,(a+b)%P};
}int f1(int n){
int x=0,y=1,d=0,n1=n;
while (n1>>=1) d++;
n1=~(n+1);
if ((n1>>d&1)) x=1;
d--;
for (;d>=0;d--){
cout<<(n1>>d&1);
int tx=(1ll*x*x+1ll*y*y)%P;
int ty=(2ll*x+y)*y%P;
if ((n1>>d)&1) x=tx,y=ty;
else x=ty,y=(tx+ty)%P;
}return x;
}typedef long long ll;
int f2(ll n){
ll x=0,y=1,t;
for (int d=62-__builtin_clz(++n);d>=0;d--){
t=(x*x+y*y)%P; y=(x+x+y)*y%P; x=t;
if (n>>d&1) x=y,y=(t+y)%P;
}return x;
}int main(){
int n=100;
for (int i=0;i<=n;i++){
cout<<(bitset<10>(i))<<' '<<(bitset<10>(~(i+1)))<<' ';
int r=f(i).first; cout<<' ';
int r2=f2(i);cout<<'\n';
cout<<r<<' '<<r2<<'\n';
}return 0;
}