【LOJ3058】「HNOI2019」白兔之舞

模数FFT与矩阵快速幂

最新推荐文章于 2019-12-28 17:15:46 发布

cz_xuyixuan

最新推荐文章于 2019-12-28 17:15:46 发布

阅读量592

点赞数

CC 4.0 BY-SA版权

分类专栏：【类型】做题记录【OJ】LOJ 【算法】FFT与NTT 【算法】单位根反演【算法】矩阵乘法

本文链接：https://blog.youkuaiyun.com/qq_39972971/article/details/89512679

【类型】做题记录同时被 3 个专栏收录

767 篇文章

订阅专栏

【资料】好题

145 篇文章

订阅专栏

【OJ】LOJ

123 篇文章

订阅专栏

本文介绍了一种解决特定矩阵幂运算问题的方法，利用任意模数FFT进行快速计算。通过求解矩阵的kkk次单位根，结合矩阵快速幂技巧，实现了O(L+N^3KLogL+KLogK)的时间复杂度。代码示例展示了如何初始化FFT，执行前缀和变换，以及如何通过多项式乘法加速计算。

【题目链接】

点击打开链接

【思路要点】

首先，求出 $L$ 的任意原根 $g$ ，以及其 $k$ 次单位根 $w$ 。
记初始给定的矩阵为 $m a t$ ，那么答案 $ans_t$ 应当满足
$ans_t=\sum_{i=0}^{L}\binom{L}{i}(mat^i)_{x,y}[i\%k=t]$
将 $[i\%k=t]$ 展开为 $\frac{1}{k}\sum_{j=0}^{k-1}w^{j(i-t)}$ ，则有
$k\times ans_t=\sum_{i=0}^{L}\binom{L}{i}(mat^i)_{x,y}\sum_{j=0}^{k-1}w^{j(i-t)}$
$k\times ans_t=\sum_{j=0}^{k-1}w^{-jt}\sum_{i=0}^{L}\binom{L}{i}w^{ij}(mat^i)_{x,y}$
$k\times ans_t=\sum_{j=0}^{k-1}w^{-jt}((w^jmat+I)^L)_{x,y}$
那么 $w^{-\binom{t}{2}}\times k\times ans_t=\sum_{j=0}^{k-1}w^{-\binom{j+t}{2}}\times w^{\binom{j}{2}}((w^jmat+I)^L)_{x,y}$
使用任意模数 $F F T$ 计算之即可。
时间复杂度 $O(\sqrt{L}+N^3KLogL+KLogK)$ 。

【代码】

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 262144;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
template <typename T> void chkmax(T &x, T y) {x = max(x, y); }
template <typename T> void chkmin(T &x, T y) {x = min(x, y); } 
template <typename T> void read(T &x) {
	x = 0; int f = 1;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
	for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
	x *= f;
}
template <typename T> void write(T x) {
	if (x < 0) x = -x, putchar('-');
	if (x > 9) write(x / 10);
	putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {
	write(x);
	puts("");
}
namespace AnyModuloFFT {
	const int MAXN = 262144;
	const long double pi = acosl(-1);
	struct point {long double x, y; };
	point operator + (point a, point b) {return (point) {a.x + b.x, a.y + b.y}; }
	point operator - (point a, point b) {return (point) {a.x - b.x, a.y - b.y}; }
	point operator * (point a, point b) {return (point) {a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x}; }
	point operator / (point a, long double x) {return (point) {a.x / x, a.y / x}; }
	int N, Log, home[MAXN];
	point tmp[MAXN];
	void FFTinit() {
		for (int i = 0; i < N; i++) {
			int tmp = i, ans = 0;
			for (int j = 1; j <= Log; j++) {
				ans <<= 1;
				ans += tmp & 1;
				tmp >>= 1;
			}
			home[i] = ans;
		}
	}
	void FFT(point *a, int mode) {
		for (int i = 0; i < N; i++)
			if (home[i] < i) swap(a[i], a[home[i]]);
		for (int len = 2; len <= N; len <<= 1) {
			point delta = (point) {cosl(2 * pi / len * mode), sinl(2 * pi / len * mode)};
			for (int i = 0; i < N; i += len) {
				point now = (point) {1, 0};
				for (int j = i, k = i + len / 2; k < i + len; j++, k++) {
					point tmp = a[j];
					point tnp = a[k] * now;
					a[j] = tmp + tnp;
					a[k] = tmp - tnp;
					now = now * delta;
				}
			}
		}
		if (mode == -1) {
			for (int i = 0; i < N; i++)
				a[i] = a[i] / N;
		}
	}
	void times(int *a, int *b, int *c, int P, int limit) {
		N = 1, Log = 0;
		while (N <= 2 * limit) {
			N <<= 1;
			Log++;
		}
		static point ax[MAXN], ay[MAXN];
		static point bx[MAXN], by[MAXN];
		for (int i = 0; i <= limit; i++) {
			ax[i] = (point) {a[i] & 32767, 0};
			ay[i] = (point) {a[i] >> 15, 0};
			bx[i] = (point) {b[i] & 32767, 0};
			by[i] = (point) {b[i] >> 15, 0};
		}
		for (int i = limit + 1; i < N; i++) {
			ax[i] = (point) {0, 0};
			ay[i] = (point) {0, 0};
			bx[i] = (point) {0, 0};
			by[i] = (point) {0, 0};
		}
		FFTinit();
		FFT(ax, 1), FFT(ay, 1), FFT(bx, 1), FFT(by, 1);
		static point x[MAXN], y[MAXN], z[MAXN];
		for (int i = 0; i < N; i++) {
			x[i] = ax[i] * bx[i];
			y[i] = ax[i] * by[i] + ay[i] * bx[i];
			z[i] = ay[i] * by[i];
		}
		FFT(x, -1), FFT(y, -1), FFT(z, -1);
		auto num = [&] (point x) {
			return (long long) (x.x + 0.5) % P;
		};
		for (int i = 0; i < N; i++) {
			int res = num(z[i]);
			res = (32768ll * res + num(y[i])) % P;
			res = (32768ll * res + num(x[i])) % P;
			c[i] = res;
		}
	}
}
int n, k, l, x, y, P, g, w;
int mat[3][3], f[MAXN], h[MAXN], res[MAXN];
int power(int x, int y) {
	if (y == 0) return 1;
	int tmp = power(x, y / 2);
	if (y % 2 == 0) return 1ll * tmp * tmp % P;
	else return 1ll * tmp * tmp % P * x % P;
}
bool PrimitiveRoot(int g) {
	int phi = P - 1;
	for (int i = 2; i * i <= phi; i++)
		if (phi % i == 0) {
			while (phi % i == 0) phi /= i;
			if (power(g, (P - 1) / i) == 1) return false;
		}
	if (phi != 1 && power(g, (P - 1) / phi) == 1) return false;
	return true;
}
int pmat() {
	static int cur[3][3], res[3][3], tmp[3][3];
	for (int i = 0; i <= n - 1; i++)
	for (int j = 0; j <= n - 1; j++) {
		res[i][j] = i == j;
		cur[i][j] = mat[i][j] + (i == j);
		if (cur[i][j] >= P) cur[i][j] -= P;
	}
	int lft = l;
	for (int p = 1; lft != 0; p <<= 1) {
		if (lft & p) {
			lft ^= p;
			for (int i = 0; i <= n - 1; i++)
			for (int j = 0; j <= n - 1; j++) {
				ll tres = 0;
				for (int k = 0; k <= n - 1; k++)
					tres += 1ll * res[i][k] * cur[k][j];
				tmp[i][j] = tres % P;
			}
			memcpy(res, tmp, sizeof(res));
		}
		for (int i = 0; i <= n - 1; i++)
		for (int j = 0; j <= n - 1; j++) {
			ll tres = 0;
			for (int k = 0; k <= n - 1; k++)
				tres += 1ll * cur[i][k] * cur[k][j];
			tmp[i][j] = tres % P;
		}
		memcpy(cur, tmp, sizeof(cur));
	}
	return res[x][y];
}
int main() {
	read(n), read(k), read(l), read(x), read(y), read(P), x--, y--;
	g = 2; while (!PrimitiveRoot(g)) g++; w = power(g, (P - 1) / k);
	for (int i = 0; i <= n - 1; i++)
	for (int j = 0; j <= n - 1; j++)
		read(mat[i][j]);
	for (int i = 0, j = 2 * k - 2; i <= 2 * k - 2; i++, j--) {
		f[i] = power(w, P - 1 - j * (j - 1ll) / 2 % (P - 1));
		if (i <= k - 1) h[i] = 1ll * pmat() * power(w, i * (i - 1ll) / 2 % (P - 1)) % P;
		for (int a = 0; a <= n - 1; a++)
		for (int b = 0; b <= n - 1; b++)
			mat[a][b] = 1ll * mat[a][b] * w % P;
	}
	AnyModuloFFT :: times(f, h, res, P, 2 * k - 2);
	int invk = power(k, P - 2);
	for (int i = 0, j = 2 * k - 2; i <= k - 1; i++, j--)
		writeln(1ll * res[j] * power(w, i * (i - 1ll) / 2 % (P - 1)) % P * invk % P);
	return 0;
}