[hoj 2255]Not Fibonacci[矩阵快速幂]

#include <stdio.h>#include <cstring>#include <iostream>#define mod 10000000using namespace std;typedef long long LL;LL a[4][4],tmp[4][4],ans[4][4];//0.00 s        1260 K/**构造矩阵                         列向量                              初始向量a[4][4] = { { 0,0,0,0 },            { 0,p,q,0 },        f(n)                               f(1) = bb            { 0,1,0,0 },        f(n-1)                             f(0) = aa            { 0,1,0,1 } };      S(n-1)(表示从f(0)到f(n-1)的和)       S(0) = f(0) = aa**/void mul(LL a[][4],LL b[][4]){    for(int i=1; i<=3; i++)        for(int j=1; j<=3; j++)        {            tmp[i][j]=0;            for(int k=1; k<=3; k++)            {                tmp[i][j]+=a[i][k]*b[k][j];                tmp[i][j]%=mod;            }        }    memcpy(a,tmp,sizeof(tmp));}void pow(LL ans[][4],LL b[][4],int n)///矩阵快速幂,返回到ans{    memset(ans,0,sizeof(ans));    for(int i=1; i<=3; i++)        ans[i][i]=1;    while(n)    {        if(n&1) mul(ans,b);        mul(b,b);        n>>=1;    }}int main(){    int t,aa,bb,p,q,s,e,a1,a2;    scanf("%d",&t);    while(t--)    {        scanf("%d%d%d%d%d%d",&aa,&bb,&p,&q,&s,&e);        memset(a,0,sizeof(a));        a[1][1]=p;        a[1][2]=q;        a[2][1]=a[3][1]=a[3][3]=1;        s--;        if(s==0) a1=aa;        else if(s<0) a1=0;        else        {            memset(ans,0,sizeof(ans));            pow(ans,a,s);            a1=(bb*ans[3][1]+aa*ans[3][2]+aa*ans[3][3])%mod;        }        if(e==0) a2=aa;        else        {            memset(ans,0,sizeof(ans));            memset(a,0,sizeof(a));            a[1][1]=p;            a[1][2]=q;            a[2][1]=a[3][1]=a[3][3]=1;            pow(ans,a,e);            a2=(bb*ans[3][1]+aa*ans[3][2]+aa*ans[3][3])%mod;        }        printf("%d\n",(a2-a1+mod)%mod);    }    return 0;}

以上代码来源:http://blog.csdn.net/ehi11/article/details/8606751  linyiming学长

题意:

定义头两项为aa,bb的数列,

f[n] = p * f[n-1] + q * f[n-2];

求此数列第s项到第e项间的和.

思路:矩阵快速幂,构造列向量表达递推式,并且可以直接求和!

自己敲一遍:

定义矩阵结构体,重载运算符.

快速幂递归实现.

#include <cstring>#include <cstdio>using namespace std;//0.00 s 728 Kconst int mod = 1e7;typedef long long ll;int aa,bb,p,q;typedef struct Matrix{    ll A[4][4];    Matrix()    {        memset(this,0,sizeof(*this));    }} Matrix;ll mymod(ll x){    return (x%mod + mod) %mod;}Matrix operator*(Matrix m1, Matrix m2){    Matrix ret;    for(int i=1; i<4; i++)        for(int j=1; j<4; j++)        {            ret.A[i][j] = 0;            for(int k=1; k<4; k++)            {                ret.A[i][j] += m1.A[i][k]*m2.A[k][j];                ret.A[i][j] = mymod(ret.A[i][j]);            }        }    return ret;}Matrix mypow(Matrix m, int n)//recursion{    Matrix ret, tmp;    if(!n)    {        for(int i=1; i<4; i++)            ret.A[i][i] = 1%mod;        return ret;    }    tmp = mypow(m, n/2);    if(n & 1)        return tmp * tmp * m;    return tmp * tmp;}ll sum(int n){    ll ans;    if(!n)        ans = aa;    else if(n==-1)        ans = 0;    else    {        Matrix m;        m.A[1][1] = p;        m.A[1][2] = q;        m.A[3][1] = m.A[3][3] = m.A[2][1] = 1;        m = mypow(m,n);        ans = mymod(bb*m.A[3][1] + aa*m.A[3][2] + aa*m.A[3][3]);    }    return ans;}/**v(n+1) = a^n * v(1)    取v(n+1)的第3项即S(n)**/int main(){    int t,s,e;    scanf("%d",&t);    while(t--)    {        scanf("%d %d %d %d %d %d",&aa,&bb,&p,&q,&s,&e);        ll as,ae;        s--;        as = sum(s);        ae = sum(e);        printf("%d\n",(int)mymod(ae-as));    }}

快速幂循环实现.

#include <cstring>#include <cstdio>using namespace std;//0.00 s 728 Kconst int mod = 1e7;typedef long long ll;int aa,bb,p,q;typedef struct Matrix{    ll A[4][4];    Matrix()    {        memset(this,0,sizeof(*this));    }} Matrix;ll mymod(ll x){    return (x%mod + mod) %mod;}Matrix operator*(Matrix m1, Matrix m2){    Matrix ret;    for(int i=1; i<4; i++)        for(int j=1; j<4; j++)        {            ret.A[i][j] = 0;            for(int k=1; k<4; k++)            {                ret.A[i][j] += m1.A[i][k]*m2.A[k][j];                ret.A[i][j] = mymod(ret.A[i][j]);            }        }    return ret;}Matrix mypow(Matrix m, int n)//loop{    Matrix ret;    for(int i=1;i<4;i++)        ret.A[i][i] = 1;    while(n)    {        if(n%2) ret = ret * m;        m = m*m;        n /= 2;    }    return ret;}ll sum(int n){    ll ans;    if(!n)        ans = aa;    else if(n==-1)        ans = 0;    else    {        Matrix m;        m.A[1][1] = p;        m.A[1][2] = q;        m.A[3][1] = m.A[3][3] = m.A[2][1] = 1;        m = mypow(m,n);        ans = mymod(bb*m.A[3][1] + aa*m.A[3][2] + aa*m.A[3][3]);    }    return ans;}/**v(n+1) = a^n * v(1)    取v(n+1)的第3项即S(n)**/int main(){    int t,s,e;    scanf("%d",&t);    while(t--)    {        scanf("%d %d %d %d %d %d",&aa,&bb,&p,&q,&s,&e);        ll as,ae;        s--;        as = sum(s);        ae = sum(e);        printf("%d\n",(int)mymod(ae-as));    }}