【ybtoj高效进阶 21279】排列计数（矩阵乘法）（光速幂）（DP）

排列计数

题目链接：ybtoj高效进阶 21279

题目大意

多次询问，每次问你有多少个长为 n 的排列满足相邻两个的差是 2 一下。

代码

考虑能否 DP，那你想他相差是 \(2\)，你考虑从 \(1\sim n\) 的数组一次取走数，每次去的位置间隔不超过 \(2\)。

考虑可以怎么走。
设 \(f_i\) 为取走前 \(i\) 个的方案数。
\(f_{i}=f_{i-1}+f_{i-3}\)

第一个是直接走，第二个是那个位置走两个过去，走一个回来，在走两个过去。
然后答案是你可能按上面走了那么多，你就两个两个飞过去，再两个两个飞回来。
所以每个 \(f\) 都要加上。

接着考虑优化 DP。
发现可以矩阵乘法。

|\(0\)|\(1\)|\(0\)|\(0\)|
|--|--|--|--|--|

\(ans\)	\(f_{i-1}\)	\(f_{i-2}\)	\(f_{i-3}\)
1	0	0	0
1	1	1	0
0	0	0	1
0	1	0	0

但是你会发现，答案会多了一点，为什么呢？
因为在 \(f_{n-3}\) 的位置直接跳过来是被重复计算的，不过还好，减去就可以了。

然后矩阵乘法不能直接用快速幂的，要用光速幂。（设 \(x\) 为询问的上界）
光速幂就是把它分成 \(\sqrt{x}\) 以内和 \(\sqrt{x}\) 的倍数这些部分预处理出来。
然后到时查询的时候就是分成这两个部分，就只需要乘两次了。

代码

#include<cstdio>
#define ll long long
#define mo 1000000007

using namespace std;

struct matrix {
	int n, m;
	ll a[5][5];
}a, b[100001], one, c[100001];
int T, n, lst;
ll ans;

matrix operator *(matrix x, matrix y) {
	matrix re;
	re.n = x.n; re.m = y.m;
	for (int i = 1; i <= re.n; i++)
		for (int j = 1; j <= re.m; j++)
			re.a[i][j] = 0;
	for (int k = 1; k <= x.m; k++)
		for (int i = 1; i <= re.n; i++)
			for (int j = 1; j <= re.m; j++)
				re.a[i][j] = (re.a[i][j] + x.a[i][k] * y.a[k][j] % mo) % mo;
	return re;
}

int main() {
//	freopen("per.in", "r", stdin);
//	freopen("per.out", "w", stdout);
	
	one.n = one.m = 4;
	one.a[1][1] = 1; one.a[2][2] = 1; one.a[3][3] = 1; one.a[4][4] = 1;
	
	b[0] = one;
	b[1].n = b[1].m = 4;
	b[1].a[1][1] = 1; b[1].a[1][2] = 0; b[1].a[1][3] = 0; b[1].a[1][4] = 0;
	b[1].a[2][1] = 1; b[1].a[2][2] = 1; b[1].a[2][3] = 1; b[1].a[2][4] = 0;
	b[1].a[3][1] = 0; b[1].a[3][2] = 0; b[1].a[3][3] = 0; b[1].a[3][4] = 1;
	b[1].a[4][1] = 0; b[1].a[4][2] = 1; b[1].a[4][3] = 0; b[1].a[4][4] = 0;
	for (int i = 2; i * i <= 1e9; i++)
		b[i] = b[i - 1] * b[1], lst = i;
	c[0] = one; c[1] = b[lst];
	for (int i = 2; i * lst <= 1e9; i++)
		c[i] = c[i - 1] * c[1]; 
	
	
	scanf("%d", &T);
	while (T--) {
		scanf("%d", &n);
		
		if (n <= 1) {
			ans ^= 1;
			continue;
		}
		
		a.n = 1; a.m = 4;
		a.a[1][1] = 0; a.a[1][2] = 1; a.a[1][3] = 0; a.a[1][4] = 0;
		a = a * c[n / lst] * b[n % lst];
		
		ans ^= (a.a[1][1] - a.a[1][4] + mo) % mo;
	}
	
	printf("%lld", ans);
	
	return 0;
}