BZOJ5298 CQOI2018 交錯序列 【DP+矩陣快速冪優化】

Description

我們稱一個僅由0、1構成的序列為”交錯序列”，當且僅當序列中沒有相鄰的1（可以有相鄰的0）。例如，000，001，101，都是交錯序列，而110則不是。對於一個長度為n的交錯序列，統計其中0和1出現的次數，分別記為x和y。給定引數a、b，定義一個交錯序列的特徵值為xayb$x^a{y}^b$。注意這裡規定任何整數的0次冪都等於1（包括00=1$0^0=1$）。顯然長度為n的交錯序列可能有多個。我們想要知道，所有長度為n的交錯序列的特徵值的和，除以m的餘數。(m是一個給定的質數)例如，全部長度為3的交錯串為：000、001、010、100、101。當a=1，b=2時，可計算

Input

輸入檔案共一行，包含三個空格分開的整數n，a，b和m。
1≤n≤10000000，0≤a，b≤45，m<100000000

Output

輸出檔案共一行，為計算結果。

Sample Input

3 1 2 1009

Sample Output

10

首先我們看到這題，果斷想到DP，但是xayb$x^a{y}^b$形式的答案十分不好維護，所以我們將它進行轉化
∑xayb$\sum x^a{y}^b$
=∑xa(n−x)b$=\sum x^a(n-x{)}^b$
=∑xa∑bi=0Cibni(−x)b−

#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define N 200
int read(){
    int ans=0,w=1;
    char c=getchar();
    while(c!='-'&&!isdigit(c))c=getchar();
    if(c=='-')w=-1,c=getchar();
    while(isdigit(c))ans=ans*10+c-'0',c=getchar();
    return ans*w;
}
int n,a,b,Mod;
LL J[N],inv[N];
struct Matrix{
    int n;LL a[N][N];
    Matrix(int n):n(n){
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                a[i][j]=0;
    }
    Matrix operator * (const Matrix &c)const{
        Matrix ans=Matrix(n);
        for(int i=0;i<n;i++)
            for(int k=0;k<n;k++)if(a[i][k])
                for(int j=0;j<n;j++)if(c.a[k][j])
                    ans.a[i][j]=(ans.a[i][j]+a[i][k]*c.a[k][j])%Mod;
        return ans; 
    }
    Matrix operator ^ (const int x)const{
        Matrix a=*this,ans=Matrix(n);
        for(int i=0;i<n;i++)ans.a[i][i]=1;
        int b=x;
        while(b){
            if(b&1)ans=ans*a;
            b>>=1;
            a=a*a;
        } 
        return ans;
    }
};
LL C(int n,int m){
    if(n<m)return 0;
    return J[n]*inv[m]%Mod*inv[n-m]%Mod;
}
LL fast_pow(LL a,LL b){
    LL ans=1;
    while(b){
        if(b&1)ans=ans*a%Mod;
        b>>=1;
        a=a*a%Mod;
    }
    return ans;
}
int main(){
    n=read();a=read();b=read();Mod=read();
    int MAX=a+b+1;
    J[0]=1;
    for(int i=1;i<=MAX;i++)J[i]=J[i-1]*i%Mod;
    inv[MAX]=fast_pow(J[MAX],Mod-2);
    for(int i=MAX-1;i>=0;i--)inv[i]=inv[i+1]*(i+1)%Mod;
    Matrix move=Matrix(MAX*2);
    for(int i=0;i<MAX;i++)move.a[i][i+MAX]=1;
    for(int i=0;i<MAX;i++)
        for(int j=0;j<=i;j++)
            move.a[j][i]=move.a[j+MAX][i]=C(i,j);
    move=move^n;
    LL x=1,ans=0;
    for(int i=0;i<=b;i++,x=x*n%Mod)
        ans=(ans+C(b,i)*x%Mod*(move.a[0][a+b-i]+move.a[0][a+b-i+MAX])%Mod*(((b-i)&1)?(-1):(1)))%Mod;
    printf("%lld",(ans+Mod)%Mod);
    return 0;
}

多謝lyw大神指點→$\to$ lyw大神