题目描述

我们称一个仅由0、1构成的序列为”交错序列“，当且仅当序列中没有相邻的1(可以有相邻的0)。例如，000,001,101,都是交错序列，而110 则不是。

对于一个长度为n 的交错序列，统计其中0 和1出现的次数，分别记为x 和y。给定参数a、b,定义一个交错序列的特征值为(x^ay^b)。注意这里规定任何整数的0 次幂都等于1(包括(0^0=1))。

显然长度为n 的交错序列可能有多个。我们想要知道，所有长度为n 的交错序列的特征值的和，除以m 的余数。(m 是一个给定的质数)

对于30%的数据，1≤n≤15

对于100%的数据，1≤n≤10000000 0≤a,b≤45 m<100000000

1.1 思路

设一个交错序列(lambda)的特征值为(s(lambda))，则我们要求的答案(ans=sum_{|lambda|=n}s(lambda))
(=sum_{|lambda|=n}x^ay^b=sum_{|lambda|=n}(n-y)^ay^b)
(=sum_{|lambda|=n}y^bsum_{k=0}^{a}dbinom{a}{k}n^ky^{a-k}(-1)^{a-k})
(=sum_{|lambda|=n}sum_{k=0}^{a}dbinom{a}{k}n^ky^{a-k+b}(-1)^{a-k})
我们记(f(k)=dbinom{a}{k}n^k(-1)^{a-k}),则原式可以写成(ans=sum_{|lambda|=n}sum_{k=0}^{a}f(k)y^{a-k+b})

注意到，把这个式子展开来看：(ans=sum_{k=0}^{a}f(k)y_{lambda}^{a-k+b}+sum_{k=0}^{a}f(k)y_{lambda'}^{a-k+b}+cdots+sum_{k=0}^{a}f(k)y_{lambda^{(cdots)}}^{a-k+b})

(=sum_{k=0}^{a}(f(k)sum_{|lambda|=n}y^{a+b-k}))

注意到对于每一个序列，(f(k))是一个常数，可以直接把(sum)提到前面去。我们只要维护一下(sum_{|lambda|=n} y_{lambda}^{a+b-k})的值，就可以快速求出答案了。

1.2 贡献的转移

题目中要统计的序列是有形式上的限制的。我们可以考虑使用dp。

我们设(f[l][p][0])表示一个以(0)结尾的，长度为(l)的序列，各序列中含有的“1”的个数的(p)次方和。也就是说，设交错序列(lambda)中含有“1”的个数为(y_{lambda})，则(f[l][p][0]=sum_{lambda_l=0}y_{lambda}^p)。(f[l][p][1])同理。

结尾为0和1的序列都可以直接在后面接0，对答案不产生贡献。由0,1结尾转移到0结尾，得到：(f[l][p][0]=f[l-1][p][0]+f[l-1][p][1])。

只有结尾为0的序列后面可以接1，对答案产生贡献。考虑贡献的转移：
(sum{y^p} ightarrow sum{(y+1)^p})，右边展开得到：(sum{(y+1)^p}=sumsum_{k=0}^{p}dbinom{p}{k}y^k=sum_{k=0}^p{dbinom{p}{k}sum{y^k}}=sum_{k=0}^pdbinom{p}{k}f[l-1][k][0])

这样我们就得到了全部的状态转移方程。

1.3 转移的优化

转移方程基本都是线性递推的，我们可以考虑使用矩阵来加速。

第一步，由于状态(f[l])和状态(f[l-1])是独立的，我们可以考虑删去这一维。记此时的状态为(f[p][0/1])。

第二步，直接将状态(f[p][0])和状态(f[p][1])合并。由于我们最多会用到(sum y^{a+b})，我们令(S=a+b+1)，设(0 leq p < S)，我们可以直接令(sp[p]=f[p][0])，(sp[p+S]=f[p][1])。这样整个状态用(sp[])来表示，被压缩到了一维。

接着考虑构造转移矩阵(T[i][j])。

承袭上面的定义，转移到0结尾的方程：(sp[p]'=sp[p]+sp[p+S]=sum sp[k]T[k][p])。我们只令(T[p][p])和(T[p+S][p])等于1就可以了。

转移到1结尾的方程：(sp[p+S]'=sum_{k=0}^p sp[k]dbinom{p}{k}=sum_{k=0}^p sp[k]T[k][p+S])。令(T[k][p+S]=dbinom{p}{k})就可以了。

初状态：序列长为0。除了(sp[0]=0^0)外，其余都为0。末状态(sp'=sp cdot T^n)。

用线性递推求组合数时初始化为0，则(dbinom{p}{k})在(k > p)时均为0，符合要求。

总代码如下：

#include<cstdio>
#include<cstring>
#include<cctype>
#define RP(i,a,b) for(register int i=a;i<=b;++i)
#define DRP(i,a,b) for(register int i=a;i>=b;--i)
#define fre(z) freopen(z".in","r",stdin),freopen(z".out","w",stdout)
#define mem(a,b) memset(a,b,sizeof(a))
#define prique priority_queue
#define rg register
#include<cstdlib>
#include<ctime>
#ifdef WIN32
#define OT "%I64d"
#else
#define OT "%lld"
#endif
using namespace std;
typedef long long ll;
typedef double db;
template<class T> T qr(T type)
{
    char ch=getchar();T ret=0;T q=1,res=1;
    while(!isdigit(ch))
        q=(ch=='-'?-1:q),ch=getchar();
    while(isdigit(ch))
        ret=ret*10+ch-'0',ch=getchar();
    if(ch=='.')
    {
        ch=getchar();
        while(isdigit(ch))
            ret=ret*10+ch-'0',res=res*10,ch=getchar();
    }
    return q==-1?-ret/res:ret/res;
}

const int maxn=10000002;
const int maxa=50,maxm=100000000;
const int maxS=100;

ll n;
int a,b,S;
ll mod;

struct mat{
    int n;
    ll v[maxS<<1][maxS<<1];
    mat()
    {
        memset(v,0,sizeof(v));
    }
};

inline mat mul(mat a,mat b)
{
    mat ans;
    ans.n=a.n;
    RP(i,0,a.n-1)
    RP(j,0,b.n-1)
    {
        ll sum=0;
        RP(k,0,a.n-1)
        {
            sum+=(a.v[i][k])*b.v[k][j];
            if(sum>=mod)sum%=mod;
        }
        ans.v[i][j]=sum;
    }
    return ans;
}

ll C[maxS<<1][maxS<<1];
ll np[maxa];

inline void init()
{
    np[0]=1;
    RP(i,1,a)
        np[i]=np[i-1]*n%mod;
    
    C[0][0]=1;
    RP(i,1,S-1)
    {
        C[i][0]=1;
        RP(j,1,i)
        C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;
    }
}

inline mat qp(mat s,ll p)
{
    mat ans;
    ans.n=s.n;
    RP(i,0,s.n-1)
        ans.v[i][i]=1;
    while(p)
    {
        if(p&1)
            ans=mul(ans,s);
        s=mul(s,s);
        p>>=1;
    }
    return ans;	
}

inline ll f(ll k)
{
    ll sgn=(a-k)&1?-1:1;
    return sgn*C[a][k]*np[k]%mod;
}

int main()
{
    //fre("B");
    n=qr(1ll),a=qr(1),b=qr(1),mod=qr(1ll);
    S=a+b+1;
    init();
    
    mat T;
    T.n=S<<1;
    
    RP(i,0,S-1)
    {
        T.v[i][i]=T.v[i+S][i]=1;
        RP(j,i,S-1)
        T.v[i][j+S]=C[j][i];
    }
    T=qp(T,n);
    
    ll ans=0;
    RP(k,0,a)
    {
        ans=(ans+f(k)*(T.v[0][a+b-k]+T.v[0][a+b-k+S])%mod)%mod;
    }
    if(ans<0)
        ans=((ans%mod)+mod)%mod;
    printf("%lld",ans);
    return 0;
}