题面

求所有长度为\(n\)的、没有相邻的1的01序列中，若0有\(x\)个、1有\(y\)个，\(x^ay^b\)之和（对\(m\)取模）。

\(n \le 10^7, m \le 10^8, 0 \le a, b \le 45\)

题解

本题麻烦的地方在于这个\(x^ay^b\)怎么处理。

所以可以求出对于所有\(i \in [0, a]\)的\(y^i\)之和，然后枚举\(i\)乘上对应的系数，加起来即可。

那么如何求出\(y^i\)之和呢？

设\(f[k][i][0/1]\)表示长度为\(k\)、结尾是0/1的序列中“1的个数”（即\(y\)）的\(i\)次方之和。

0结尾的序列可以从0/1序列转移过来，而1的出现次数不会变。

1结尾的序列只能从0结尾的转移过来，1的出现次数会+1，也就是新的\(y' = (y + 1)^i = \sum_{j = 0}^{i}C_i^j y\)。

然后构建矩阵就可以做了！

#include <cstdio>

#include <cstring>

#include <cmath>

#include <algorithm>

#include <iostream>

#include <set>

#define enter putchar('\n')

#define space putchar(' ')

using namespace std;

typedef long long ll;

template <class T>

void read(T &x){

    char c;

    bool op = 0;

    while(c = getchar(), c < '0' || c > '9')

    if(c == '-') op = 1;

    x = c - '0';

    while(c = getchar(), c >= '0' && c <= '9')

    x = x * 10 + c - '0';

    if(op == 1) x = -x;

}

template <class T>

void write(T x){

    if(x < 0) putchar('-'), x = -x;

    if(x >= 10) write(x / 10);

    putchar('0' + x % 10);

}

const int N = 185;

int n, a, b, P, sze1, sze2;

ll c[N][N], ans;

struct matrix {

    ll g[N][N];

    matrix(){

	memset(g, 0, sizeof(g));

    }

    matrix operator * (const matrix &b) const {

	matrix c;

	for(int i = 0; i < sze2; i++)

	    for(int j = 0; j < sze2; j++){

		for(int k = 0; k < sze2; k++)

		    c.g[i][j] += g[i][k] * b.g[k][j];

		c.g[i][j] %= P;

	    }

	return c;

    }

    friend matrix qpow(matrix a, int x){

	matrix ret;

	for(int i = 0; i < sze2; i++)

	    ret.g[i][i] = 1;

	while(x){

	    if(x & 1) ret = ret * a;

	    a = a * a;

	    x >>= 1;

	}

	return ret;

    }

} op;

int main(){

    read(n), read(a), read(b), read(P);

    sze1 = a + b + 1, sze2 = 2 * sze1;

    c[0][0] = 1;

    for(int i = 1; i <= a + b; i++){

	c[i][0] = 1;

	for(int j = 1; j <= i; j++)

	    c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % P;

    }

    for(int i = 0; i < sze1; i++){

	op.g[i][i] = op.g[i][sze1 + i] = 1;

	for(int j = 0; j <= i; j++)

	    op.g[sze1 + i][j] = c[i][j];

    }

    op = qpow(op, n);

    ll pw = 1;

    for(int i = 0; i <= a; i++){

	ans += (((a - i) & 1) ? -1 : 1) * c[a][i] * pw % P * (op.g[a + b - i][0] + op.g[sze1 + a + b - i][0]) % P;

	pw = pw * n % P;

    }

    write((ans % P + P) % P), enter;

    return 0;

}