problem

luogu-P4456

solution

预处理阶乘和阶乘的逆元，枚举 $1$ 出现次数 $i$，$\sum \left(\genfrac{}{}{0ex}{}{n-i+1}{i}\right)(n-i{)}^a{i}^b$。

• $\left(\genfrac{}{}{0ex}{}{n-i+1}{i}\right)$ 如何推出来?

  从 $n$ 个中选 $i$ 个 $\left(\genfrac{}{}{0ex}{}{n}{i}\right)$。容斥不太可能。

  隔板法。分成 $i$ 堆需要插 $i-1$ 块板。

  考虑 $1$ 不能连续也就是说两块板之间至少要隔两个盒子。

  隔板法经典强制每个至少有 $1$ 个的方法类比过来。

  发现只需要拿走 $i-1$ 个盒子，剩下的选的位置就可以相邻。

但这样是不行的，快速幂部分耗时大，且 $n,m$ 之间大小关系不定，逆元可能不存在，无法预处理。

考虑解决组合数计算部分。

• 用 $\text{lucas}$ 定理算。

  不能在要求时间内通过。

• 算一个组合数，可以将分子分母质因数分解，然后在指数位置进行加减运算。最后在把所有质因数乘起来，就可以巧妙避免计算逆元的问题。

  但本题组合数要求计算多个，时间开销依然未能缩减。

考虑解决快速幂计算部分。

• $x^a,y^b$ 都是完全积性函数。

  $f(i)=i^a\Rightarrow f(xy)=(xy{)}^a=f(x)f(y)=x^a{y}^a$。

  所以可以 $O(n)$ 线性筛。

  这是解决幂运算指数固定的常见方法。

这种做法是不能通过本题的，洛谷题解有一篇能过完全是数据问题。这里只是想记录一些 $\text{trick}$。

---

考虑计算贡献 $x^a{y}^b=(n-y{)}^a{y}^b={\sum }_{i=0}^a(\genfrac{}{}{0ex}{}{a}{i})n^i(-1{)}^{a-i}{y}^{a-i}{y}^b$。

发现当枚举 $i$ 时 $\left(\genfrac{}{}{0ex}{}{a}{i}\right)n^i(-1{)}^{a-i}$ 均为常数，唯一随序列不同而变化的是 $y^{a+b-i}$，准确来说应该是 $y$。

我们可以计算所有序列中 $1$ 的个数的 $a+b-i$ 次方之和，然后就可以计入答案。

设 $f(i,j,k):$ 考虑前 $i$ 位，第 $i$ 位为 $k\in [0,1]$，所有合法序列的 $1$ 的个数的 $j$ 次方之和，即 $\sum y^j$。

• $i$ 填 $0$，则对前面无限制。没有新增的 $1$ 的个数，贡献不变。

  $f(i,j,0)=f(i-1,j,0)+f(i-1,j,1)$。

• $i$ 填 $1$，则前一位不能为 $1$，只能从 $0$ 转移。

  此时将有 $y^j\to (y+1{)}^j$ 直接二项式展开。

  $(y+1{)}^j={\sum }_{k=0}^j(\genfrac{}{}{0ex}{}{j}{k})y^k{1}^{j-k}$。

  所有序列的 $y^k$ 之和恰恰是 $f(,k,)$ 的定义。

  所以转移为：$f(i,j,1)={\sum }_{k=0}^j\left(\genfrac{}{}{0ex}{}{j}{k}\right)f(i-1,k,0)$。

发现转移压根和 $i$ 这一维没有关系，所以是可以矩阵加速 $n$ 的。

注意我们要计算到 $y^{a+b}$ 次方，且我们将两种转移合并在一起。

构造初始矩阵 $f:[f_{0,0},f_{1,0},...,f_{a+b,0},f_{0,1},f_{1,1},...,f_{a+b,1}]$。形式化为 $[f_{j,0}\mid f_{j,1}],j\in [0,a+b]$。

构造加速矩阵 $g$：分拆为四个部分。

• 左上角为单位矩阵。表示 $f(,0)\to f^{\prime }(,0)$。
• 左下角为单位矩阵。表示 $f(,1)\to f^{\prime }(,0)$。
• 右上角为组合数矩阵。注意是行列交换了的，表示 $f(,0)\to f^{\prime }(,1)$。
• 右下角为全 $0$ 矩阵。表示不合法 $f(,1)\to f^{\prime }(,1)$。

具体可以自己画一下，发现是匹配的。

code

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define maxn 185
int n, a, b, mod, m1, m2;
int C[maxn][maxn];
struct matrix {int c[maxn][maxn];matrix() { memset( c, 0, sizeof( c ) ); }matrix operator * ( matrix &v ) {matrix ans;for( int i = 0;i < m2;i ++ )for( int k = 0;k < m2;k ++ )if( c[i][k] ) //稀疏矩阵经典有效优化for( int j = 0;j < m2;j ++ ) //j,k交换 内存访问连续 优化常数ans.c[i][j] = (ans.c[i][j] + c[i][k] * v.c[k][j]) % mod;return ans;}	
}g, f;signed main() {scanf( "%lld %lld %lld %lld", &n, &a, &b, &mod );m1 = a + b + 1, m2 = m1 << 1;for( int i = 0;i <= m1;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]) % mod;}for( int i = 0;i < m1;i ++ ) {g.c[i][i] = g.c[i + m1][i] = 1;for( int j = i;j < m1;j ++ )g.c[i][j + m1] = C[j][i];}f.c[0][0] = 1; int x = n;while( x ) {if( x & 1 ) f = f * g;g = g * g;x >>= 1;}x = 1; int ans = 0;for( int i = 0;i <= a;i ++, x = x * n % mod )if( (a - i) & 1 )(ans -= (f.c[0][a+b-i] + f.c[0][a+b-i+m1]) % mod * C[a][i] % mod * x) %= mod;else (ans += (f.c[0][a+b-i] + f.c[0][a+b-i+m1]) % mod * C[a][i] % mod * x) %= mod;printf( "%lld\n", (ans + mod) % mod );return 0;
}