[CQOI2018] 交错序列（矩阵加速优化dp）

ikrvxt

于 2022-03-28 18:36:18 发布

阅读量597

点赞数 1

分类专栏： DP # 矩阵加速文章标签： dp 组合数二项式定理矩阵加速卡常

本文链接：https://blog.csdn.net/Emm_Titan/article/details/123802493

版权

DP 同时被 2 个专栏收录

72 篇文章 3 订阅

订阅专栏

矩阵加速

13 篇文章 0 订阅

订阅专栏

problem

luogu-P4456

solution

预处理阶乘和阶乘的逆元，枚举 $1$ 出现次数 $i$ ， $\sum\binom{n-i+1}{i}(n-i)^ai^b$ 。

$\binom{n-i+1}{i}$ 如何推出来?

从 $n$ 个中选 $i$ 个 $\binom ni$ 。容斥不太可能。

隔板法。分成 $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^ay^a$ 。

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

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

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

考虑计算贡献 $x^ay^b=(n-y)^ay^b=\sum_{i=0}^a\binom ain^i(-1)^{a-i}y^{a-i}y^b$ 。

发现当枚举 $i$ 时 $\binom{a}{i}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\rightarrow (y+1)^j$ 直接二项式展开。

$(y+1)^j=\sum_{k=0}^j\binom jky^k1^{j-k}$ 。

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

所以转移为： $f(i,j,1)=\sum_{k=0}^j\binom jkf(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)\rightarrow f'(,0)$ 。
左下角为单位矩阵。表示 $f(,1)\rightarrow f'(,0)$ 。
右上角为组合数矩阵。注意是行列交换了的，表示 $f(,0)\rightarrow f'(,1)$ 。
右下角为全 $0$ 矩阵。表示不合法 $f(,1)\rightarrow f'(,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;
}