hdu 4686 Arc of Dream 杭电多校第九场1001 矩阵

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL __int64
const LL mod=1e9+7;
LL n;
struct matrix{
    LL f[5][5];
};
matrix mul(matrix a,matrix b)
{
    matrix c;
    memset(c.f,0,sizeof(c.f));
    LL i,j,k;
    for(k=0;k<5;k++)
    for(i=0;i<5;i++)
    {
        if(!a.f[i][k])continue;
        for(j=0;j<5;j++)
        {
            if(!b.f[k][j])continue;
            c.f[i][j]=(c.f[i][j]+a.f[i][k]*b.f[k][j]%mod)%mod;
        }
    }
    return c;
}
matrix pow_mod(matrix a,LL b)
{
    matrix s;
    LL i,j;
    for(i=0;i<5;i++)
    {
        for(j=0;j<5;j++)
        {
            s.f[i][j]=(i==j?1:0);
            //cout<<s.f[i][j]<<" ";
        }
        //cout<<endl;
    }
    while(b)
    {
        if(b&1)
            s=mul(s,a);
        a=mul(a,a);
        b=b>>1;
    }
   /* for(i=0;i<5;i++)
    {
        for(j=0;j<5;j++)
            cout<<s.f[i][j]<<" ";
        cout<<endl;
    }*/
    return s;
}
int main()
{
    //freopen("D:\\test.txt","r",stdin);
	//freopen("D:\\in.txt","w",stdout);
    LL p,q,x,y,a,b;
    while(cin>>n)
    {
        cin>>a>>x>>y>>b>>p>>q;
        matrix e;
        memset(e.f,0,sizeof(e.f));
        e.f[0][0]=e.f[1][0]=e.f[4][4]=1;
        e.f[1][1]=x*p%mod;e.f[2][1]=x*q%mod;e.f[2][2]=x%mod;e.f[3][1]=p*y%mod;
        e.f[3][3]=p%mod;e.f[4][1]=q*y%mod;e.f[4][2]=y%mod;e.f[4][3]=q%mod;
        /*for(LL i=0;i<5;i++)
        {
            for(LL j=0;j<5;j++)
            cout<<e.f[i][j]<<" ";
            cout<<endl;
        }*/
        e=pow_mod(e,n);
        LL ans;
        ans=(e.f[1][0]*a%mod*b%mod+a*e.f[2][0]%mod+b*e.f[3][0]%mod+e.f[4][0]%mod)%mod;
        cout<<(ans+mod)%mod<<endl;
    }
    return 0;
}
/*
    |1 0   0 0 0 |
    |1 x*p 0 0 0 |
    |0 x*q x 0 0 |*|s[n] f[n] a[n] b[n] 1|=|s[n+1] f[n+1] a[n+1] b[n+1] 1|
    |0 p*y 0 p 0 |
    |0 q*y y q 1 |
    
    a[n+1]=a[n]*x+y;
    b[n+1]=b[n]*p+q;
    f[n]=a[n]*b[n];
    s[n+1]=s[n]+f[n];
*/

我一开始直接把公式推出来了，没想过矩阵，欲哭无泪啊，我还在想大神就是大神，10分钟就推好公式写好代码了。。，下面是公式法：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL __int64
const LL mod=1e9+7;
LL pow_mod(LL a,LL b)
{
    LL s=1;
    a=a%mod;
    while(b)
    {
        if(b&1)
            s=(s*a)%mod;
        a=(a*a)%mod;
        b=b>>1;
       // cout<<a<<endl;
    }
    return s;
}
LL divs(LL a,LL b)
{
    //cout<<a<<" "<<pow_mod(b,mod-2)<<endl;
    LL s;
    a=a%mod;
    s=(a*pow_mod(b,mod-2))%mod;
    return s;
}
void swap(LL &a,LL &b)
{
    LL t=a;
    a=b;
    b=t;
}
int main()
{
	//freopen("D:\\test.txt","r",stdin);
	//freopen("D:\\out.txt","w",stdout);
    LL n,a,x,y,b,p,q,c1,c2,c3,c4,m,z1,z2;
    while(cin>>n>>a>>x>>y>>b>>p>>q)
    {
        if(n==0)
        {
            cout<<0<<endl;
            continue;
        }
        LL i,j,k,ans=0;
        z1=divs((n%mod)*((n-1)%mod)%mod,2);
        z2=divs((((n-1)%mod))*(n%mod)%mod*((n%mod*2-1)%mod)%mod,6);
        if(x==1&&p==1)
        {
            //cout<<divs(n*(n-1),2)<<endl;
            ans=(ans+n%mod*a%mod*b%mod)%mod;
            ans=(ans+(z1*((q*a%mod+b*y%mod)%mod))%mod)%mod;
            ans=(ans+z2*(q*y%mod)%mod)%mod;
        }
        else if(p==1)
        {
            LL t=pow_mod(x,n);
            ans=(ans+divs(n%mod*y%mod*b%mod,1-x))%mod;
            ans=(ans+divs(z1*q%mod*y%mod,1-x))%mod;
            ans=(ans+divs((1-t)*(a*b%mod*(1-x)%mod-y*b%mod)%mod,(1-x)*(1-x)))%mod;
            ans=(ans+divs((divs(x-t,(1-x)*(1-x))-divs((n-1)%mod*t%mod,1-x))*(q*a%mod*(1-x)%mod-y*q%mod),1-x))%mod;
        }
        else if(x==1)
        {
            swap(x,p);swap(y,q);swap(a,b);
            LL t=pow_mod(x,n);
            ans=(ans+divs(n%mod*y%mod*b%mod,1-x))%mod;
            ans=(ans+divs(z1*q%mod*y%mod,1-x))%mod;
            ans=(ans+divs((1-t)*(a*b%mod*(1-x)%mod-y*b%mod)%mod,(1-x)*(1-x)))%mod;
            ans=(ans+divs((divs(x-t,(1-x)*(1-x))-divs((n-1)%mod*t%mod,1-x))*(q*a%mod*(1-x)%mod-y*q%mod),1-x))%mod;
        }
        else
        {
            LL e,f,g,h,d,l1,l2,l3,aa,bb,cc,dd;
            d=x*p;
            e=(1-x)*(1-p);
            f=pow_mod(x*p%mod,n);g=pow_mod(x,n),h=pow_mod(p,n);
            l1=divs(1-f,1-d);l2=divs(1-h,1-p);l3=divs(1-g,1-x);
            aa=y*q%mod;
            bb=y*b%mod;
            cc=q*a%mod;
            ans=(ans+divs(q*y%mod,e)*(n%mod)%mod)%mod;
            ans=(ans+divs(l1*(aa-bb*(1-p)%mod-q*a%mod*(1-x)%mod+a*b%mod*(1-x)%mod*(1-p)%mod)%mod,e)%mod)%mod;
            ans=(ans+divs(l2*(bb*(1-p)%mod-aa)%mod,e)%mod)%mod;
            ans=(ans+divs(l3*(cc*(1-x)%mod-aa)%mod,e)%mod)%mod;
        }
        cout<<(ans+mod)%mod<<endl;
    }
    return 0;
}