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 long
const 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));
}