长度为n的01序列的的一种方案中，保证1不能相邻，记0出现x次，1出现y次，定义其特征值为(x^ay^b),参数a，b已事先给定，求所有方案的特征值(mod) m,m为给定质数，1≤n≤10000000 0≤a,b≤45 m<100000000。
解
对于一种方案分析
[ecause x+y==n herefore x^ay^b=(n-y)^ay^b=
]
[sum_{i=0}^aC_a^in^i(-y)^{a-i}y^b=sum_{i=0}^a(-1)^{a-i}C_a^in^iy^{a+b-i}
]
注意到a，b很小，而前面式子为定值，于是只要求出每种方案的1的个数的次方即可累加即可，所以设(f[i][j][0/1])表示填到第i个数，这个数填0/1，的方案数的所有1的个数的j次方之和，不难有
[f[i][j][0]=f[i-1][j][0]+f[i-1][j][1]
]
[f[i][j][1]=sum_{k=0}^jf[i-1][k][0]C_J^k
]
注意到数据范围很大，必然(O((a+b)n))会超时，考虑到第二维第三维很少，方程较为简单，于是矩阵快速幂优化，压维存下第二维第三维即可，此外考虑到矩阵维数很少，模数较小，于是可以矩阵乘法实现时，加完一个for询问再取模。
参考代码：
#include <iostream>
#include <cstdio>
#include <cstring>
#define il inline
#define ri register
#define ll long long
using namespace std;
ll n,a,b,lsy,mid,li,
c[91][91],pown[91],
qf[91],ans;
struct matrix{
ll num[183][183];
il void clear(){memset(num,0,sizeof(num));}
il void unit(){
clear();for(ri int i(0);i<=li;++i)
num[i][i]|=true;
}
il void print(){
for(ri int i(0),j;i<=li;++i,putchar('
'))
for(j=0;j<=li;++j)printf("%d ",num[i][j]);
putchar('
');
}
il matrix operator*(matrix x){
matrix y;y.clear();
for(ri int i(0),j,k;i<=li;++i)
for(j=0;j<=li;y.num[i][j]%=lsy,++j)
for(k=0;k<=li;++k)
y.num[i][j]+=num[i][k]*x.num[k][j];
return y;
}template<class free>
il matrix operator^(free y){
matrix ans,x(*this);ans.unit();
while(y){
if(y&1)ans=ans*x;
x=x*x,y>>=1;
}return ans;
}
}state,tran;
int main(){
int i,j,k;
scanf("%lld%lld%lld%lld",&n,&a,&b,&lsy);
mid=a+b,li=(mid<<1)+1;
for(i=0;i<=mid;++i){
c[i][0]=1;
for(j=1;j<=i;++j)
c[i][j]=(c[i-1][j]+c[i-1][j-1])%lsy;
}
for(i=0;i<=mid;++i)
tran.num[i][i]=tran.num[i+mid+1][i]=1;
for(j=mid+1;j<=li;++j)
for(k=0;k<=j-mid-1;++k)
tran.num[k][j]=c[j-mid-1][k];
for(i=mid+1;i<=li;++i)state.num[0][i]=1;
state.num[0][0]=1,state=state*(tran^n-1);
for(i=0;i<=mid;++i)state.num[0][i]+=state.num[0][i+mid+1];
for(i=pown[0]=qf[0]=1;i<=mid;++i)pown[i]=pown[i-1]*n%lsy,qf[i]=qf[i-1]*(-1);
for(j=0;j<=a;++j)ans+=(pown[j]*c[a][j]%lsy*qf[a-j]+lsy)%lsy*state.num[0][a+b-j]%lsy,
ans%=lsy;
printf("%lld",ans);
return 0;
}