51nod 1341 混合序列 （矩阵快速幂）

https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1341

1341 混合序列

基准时间限制：1 秒 空间限制：131072 KB 分值: 80 难度：5级算法题

当给定p,q,r时，我们可以定义

对于给定的p,q,r,n，请计算

对于p=2 q=2 r=1 n=1这组数据，

所以答案是

Input

共1行，4个整数数p, q, r, n中间用空格分隔(1 <= p, q, r, n<=1000000000)。

Output

对于每一个数据，在一行中输出答案。

Input示例

2 2 1 1

Output示例

3

题意：a[0]=0,a[n]=a[n-1]*p+r;  b[0]=3,b[n]=b[n-1]*q;   求

思路：a[0]*b[n]+a[1]*b[n-1]+......a[n]*b[0] .   因为a[0]和b[n]相乘,可以想到，能不能把b[n]变成b[0],b[0]变成b[n]

这样就能a[0]*b[0]了，b[n]=3*q^n,b[n-1]=b[n]/q; 因为%mod里不能出现出发，现在假设Y是q对mod的逆元，假设

X=3*q^n,那么b数列的表达式就变成了b[0]=XY^0,b[n]=b[n-1]*Y;

现在我假设b[0]=Y^0,b[n]=b[n-1]*Y;

现在求X*a[0]*b[0]+X*a[1]*b[1]+...X*a[n]*b[n];(把X提出来也一样怎么处理都行）

设sum[n]表示前n项和，

那么sum[n]=sum[n-1]+a[n]*b[n];   (a[n]*b[n]这种类型未知）

a[n]*b[n]=X*(a[n-1]*p+r)*(b[n-1]*Y)=a[n-1]*b[n-1]*p*Y*X+b[n-1]*r*Y*X;  （b[n]这种类型未知）

b[n]=b[n-1]*Y;   (现在左右的类型都已知）所以第n-1个矩阵，系数矩阵，和第n个矩阵如下

sum[n-1]   a[n-1]*b[n-1]   b[n-1]   

#########################

1           0            0

p*Y        p*Y         0

X*Y*r     X*Y*r      Y

#########################

sum[n]        a[n]*b[n]         b[n]  

然后矩阵快速幂轻松解决。

#include<stdio.h>#include<algorithm>#include<queue>#include<iostream>#include<fstream>#include<cstring>#include<math.h>#define LL long longconst LL mod=1e9+7LL;using namespace std;struct node{    LL a[3][3];    void mem()    {        memset(a,0,sizeof(a));    }};node operator *(node a,node b){    node c;    c.mem();    for(int i=0;i<3;i++)        for(int j=0;j<3;j++)        {            for(int k=0;k<3;k++)                c.a[i][j]+=a.a[i][k]*b.a[k][j];            c.a[i][j]%=mod;        }        return c;}LL poww(LL a,LL b){    LL ans=1LL;    while(b)    {        if(b&1) ans=ans*a%mod;        a=a*a%mod;        b/=2;    }    return ans;}LL inv(LL a)//求a对mod的逆元{    return a==1?1:inv(mod%a)*(mod-mod/a)%mod;}LL pow(LL p,LL q,LL r,LL n){    LL Y=inv(q),X=3*poww(q,n)%mod;    node ans,a;    ans.mem();    a.mem();    ans.a[0][2]=1;    a.a[0][0]=1;    a.a[1][0]=a.a[1][1]=p*Y%mod;    a.a[2][0]=a.a[2][1]=X*Y%mod*r%mod;a.a[2][2]=Y;    while(n)    {        if(n&1) ans=ans*a;        a=a*a;        n/=2;    }    return ans.a[0][0];}int main(){    LL p,q,r,n;    scanf("%I64d%I64d%I64d%I64d",&p,&q,&r,&n);    printf("%I64d\n",pow(p,q,r,n));}