顯然答案為Σk^b·(n-k)^a·C(n-k+1,k)。並且可以發現ΣC(n-k,k)=fib_n。但這實際上沒有任何卵用。

　　純組合看起來不太行得通，換個思路，考慮一個顯然的dp，即設f[i][j][0/1]為前i為選了j個1其中第i位是0/1的方案數。這樣當然能求出答案，複雜度O(n²)。

　　注意到ab很小，並且事實上我們並不需要知道所有的方案數，而是隻要求出貢獻就可以了。而又有x^ay^b=x^a(n-x)^b，這個式子顯然只要求出所有Σxⁱ就能求了。再由二項式定理，(k+1)^b=ΣC(b,i)kⁱ。那麼做法就比較顯然了，維護上述矩陣大力矩乘即可。

　　非常卡常。在bzoj排倒數第二2333

#include<iostream> 
#include
 
<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
#define N 190
char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<'0'||c>'9')) c=getchar();return
 
 c;}
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int n,m,A,B,P,ans,C[N][N];
inline 
 
void inc(int &x,int y){x+=y;if (x>=P) x-=P;}
int ksm(int a,int k)
{
    int s=1;
    for (;k;k>>=1,a=1ll*a*a%P) if (k&1) s=1ll*s*a%P;
    return s;
}
struct matrix
{
    int n,a[N][N];
    matrix operator *(const matrix&b) const
    {
        matrix c;c.n=n;memset(c.a,0,sizeof(c.a));
        for (int i=0;i<n;i++)
            for (int j=0;j<b.n;j++)
                for (int k=0;k<b.n;k++)
                inc(c.a[i][j],1ll*a[i][k]*b.a[k][j]%P);
        return c;
    }
}f,a;
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj5298.in","r",stdin);
    freopen("bzoj5298.out","w",stdout);
    const char LL[]="%I64d\n";
#else
    const char LL[]="%lld\n";
#endif
    m=n=read(),A=read(),B=read(),P=read();
    C[0][0]=1;
    for (int i=1;i<N;i++)
    {
        C[i][0]=C[i][i]=1;
        for (int j=1;j<i;j++)
        C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
    }
    f.n=1;f.a[0][0]=1;
    a.n=A+B+1<<1;
    for (int i=0;i<=A+B;i++) a.a[i][i]++,a.a[i+A+B+1][i]++;
    for (int i=0;i<=A+B;i++)
        for (int j=i;j<=A+B;j++)
        inc(a.a[i][j+A+B+1],C[j][i]);
    for (;n;n>>=1,a=a*a) if (n&1) f=f*a;
    for (int i=0;i<=A+B;i++) inc(f.a[0][i],f.a[0][i+A+B+1]);
    for (int i=B;i<=A+B;i++)
    if (i-B&1) inc(ans,P-1ll*f.a[0][i]*ksm(m,A+B-i)%P*C[A][i-B]%P);
    else inc(ans,1ll*f.a[0][i]*ksm(m,A+B-i)%P*C[A][i-B]%P);
    cout<<ans;
    return 0;
}