【LOJ3058】「HNOI2019」白兔之舞

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

【题目链接】

【思路要点】

  • 首先,求出 L L L 的任意原根 g g g ,以及其 k k k 次单位根 w w w
  • 记初始给定的矩阵为 m a t mat mat ,那么答案 a n s t ans_t anst 应当满足
    a n s t = ∑ i = 0 L ( L i ) ( m a t i ) x , y [ i % k = t ] ans_t=\sum_{i=0}^{L}\binom{L}{i}(mat^i)_{x,y}[i\%k=t] anst=i=0L(iL)(mati)x,y[i%k=t]
  • [ i % k = t ] [i\%k=t] [i%k=t] 展开为 1 k ∑ j = 0 k − 1 w j ( i − t ) \frac{1}{k}\sum_{j=0}^{k-1}w^{j(i-t)} k1j=0k1wj(it) ,则有
    k × a n s t = ∑ i = 0 L ( L i ) ( m a t i ) x , y ∑ 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×anst=i=0L(iL)(mati)x,yj=0k1wj(it)
    k × a n s t = ∑ j = 0 k − 1 w − j t ∑ i = 0 L ( L i ) w i j ( m a t i ) x , y 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×anst=j=0k1wjti=0L(iL)wij(mati)x,y
    k × a n s t = ∑ j = 0 k − 1 w − j t ( ( w j m a t + I ) L ) x , y k\times ans_t=\sum_{j=0}^{k-1}w^{-jt}((w^jmat+I)^L)_{x,y} k×anst=j=0k1wjt((wjmat+I)L)x,y
  • 那么 w − ( t 2 ) × k × a n s t = ∑ j = 0 k − 1 w − ( j + t 2 ) × w ( j 2 ) ( ( w j m a t + 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} w(2t)×k×anst=j=0k1w(2j+t)×w(2j)((wjmat+I)L)x,y
  • 使用任意模数 F F T FFT FFT 计算之即可。
  • 时间复杂度 O ( L + N 3 K L o g L + K L o g K ) O(\sqrt{L}+N^3KLogL+KLogK) O(L +N3KLogL+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;
}
<think>我们正在查询与LOJ6279相关的编程题目或解决方案。LOJ(LibreOJ)是一个在线的评测系统,主要收录算法竞赛题目。根据题号6279,我们需要确定该题目的具体内容。由于我无法直接访问网络,我将基于已知信息进行推理。在算法竞赛中,LOJ的题目编号通常与特定的题目对应。题号6279可能对应一个特定的问题。根据常见的题目类型,我们可以推测它可能涉及数据结构或算法。然而,为了准确回答,我们可以尝试回忆或推断常见的LOJ题目。例如,LOJ6000系列题目多为数据结构题,如线段树、树状数组等。实际上,通过搜索记忆,我回想起LOJ6279是“数列分块入门3”(BlockArray3)题目。这是“分块”系列题目中的一题。题目描述:给出一个长为n的数列,以及n个操作,操作涉及区间加法,询问区间内小于某个值x的前驱(即最大的小于x的数)。因此,用户需要的是关于LOJ6279题目的描述和解决方案。解决方案:由于题目要求分块处理,我们将数列分成若干块,每块大小为sqrt(n)。对于每个块,我们维护一个有序的数组(用于二分查找)以及一个加法标记(记录整个块被加了多少)。操作:1.区间加法:对于整块,直接更新加法标记;对于不完整的块,暴力更新元素值,然后重新排序该块(或者重新构建有序数组)。2.查询区间内小于x的前驱:遍历区间内的所有块。对于整块,在块的有序数组中二分查找小于(x-加法标记)的最大值;对于不完整的块,暴力遍历每个元素,判断是否小于x(注意加上加法标记)并更新最大值。代码实现(C++)示例:```cpp#include<iostream>#include<vector>#include<cmath>#include<algorithm>#include<climits>usingnamespacestd;constintMAXN=100010;constintBLOCK_SIZE=320;//约sqrt(MAXN)intn;vector<int>a;vector<int>block_add;//每块的加法标记vector<vector<int>>sorted_blocks;//每块的有序数组voidinit(){intblock_num=(n+BLOCK_SIZE-1)/BLOCK_SIZE;block_add.assign(block_num,0);sorted_blocks.resize(block_num);for(inti=0;i<n;i++){sorted_blocks[i/BLOCK_SIZE].push_back(a[i]);}for(inti=0;i<block_num;i++){sort(sorted_blocks[i].begin(),sorted_blocks[i].end());}}voidupdate_block(intblock_id){sorted_blocks[block_id].clear();intstart=block_id*BLOCK_SIZE;intend=min((block_id+1)*BLOCK_SIZE,n);for(inti=start;i<end;i++){sorted_blocks[block_id].push_back(a[i]);}sort(sorted_blocks[block_id].begin(),sorted_blocks[block_id].end());}voidadd(intl,intr,intc){intblock_l=l/BLOCK_SIZE,block_r=r/BLOCK_SIZE;if(block_l==block_r){for(inti=l;i<=r;i++){a[i]+=c;}update_block(block_l);}else{for(inti=l;i<(block_l+1)*BLOCK_SIZE;i++){a[i]+=c;}update_block(block_l);for(inti=block_l+1;i<block_r;i++){block_add[i]+=c;}for(inti=block_r*BLOCK_SIZE;i<=r;i++){a[i]+=c;}update_block(block_r);}}intquery(intl,intr,intc){intblock_l=l/BLOCK_SIZE,block_r=r/BLOCK_SIZE;intans=-1;if(block_l==block_r){for(inti=l;i<=r;i++){intval=a[i]+block_add[block_l];if(val<c){if(val>ans)ans=val;}}}else{for(inti=l;i<(block_l+1)*BLOCK_SIZE;i++){intval=a[i]+block_add[block_l];if(val<c){if(val>ans)ans=val;}}for(inti=block_l+1;i<block_r;i++){//在块i的有序数组中二分查找小于c-block_add[i]的最大值inttarget=c-block_add[i];autoit=lower_bound(sorted_blocks[i].begin(),sorted_blocks[i].end(),target);if(it!=sorted_blocks[i].begin()){it--;intval=*it+block_add[i];if(val>ans)ans=val;}}for(inti=block_r*BLOCK_SIZE;i<=r;i++){intval=a[i]+block_add[block_r];if(val<c){if(val>ans)ans=val;}}}returnans;}intmain(){cin>>n;a.resize(n);for(inti=0;i<n;i++){cin>>a[i];}init();for(inti=0;i<n;i++){intop,l,r,c;cin>>op>>l>>r>>c;l--;r--;//如果题目输入是1-indexed,则转换为0-indexedif(op==0){add(l,r,c);}else{cout<<query(l,r,c)<<endl;}}return0;}```注意:上述代码是示例代码,可能需要根据题目具体输入要求进行调整。§§相关问题§§1.LOJ6279题目的具体描述是什么?2.分块算法的时间复杂度是多少?3.除了分块,还有哪些方法可以解决LOJ6279?4.如何优化分块算法的常数?5.分块算法适用于哪些场景?</think>### LOJ 6279 题目分析与解决方案 **题目描述**:LOJ 6279 是 LibreOJ 上的 **"数列分块入门 3"**,属于分块算法的经典练习题。核心要求如下: - **操作 1**:区间 $[l, r]$ 内每个数加上值 $c$ - **操作 2**:查询区间 $[l, r]$ 内小于 $c$ 的**最大元素值**(即前驱),不存在时返回 $-1$ #### 分块算法设计 分块的核心思想是将数列划分为 $\sqrt{n}$ 大小的块,平衡暴力操作与批量处理。 ##### 关键数据结构 1. **块内有序数组**: - 每块维护一个排序后的数组,用于二分查找前驱。 - 示例:块 $B_i$ 的排序数组为 $S_i$。 2. **加法懒标记**: - 记录整块被加的总值 $add_i$(不影响块内有序性)。 ##### 操作实现 1. **区间加法($[l, r] + c$)**: - **整块**:直接更新懒标记 $add_i \gets add_i + c$,时间复杂度 $O(1)$。 - **碎块**:暴力修改元素值,并**重构有序数组**(排序),时间复杂度 $O(\sqrt{n} \log \sqrt{n})$。 2. **查询前驱($[l, r]$ 内 $< c$ 的最大值)**: - **整块**:在有序数组 $S_i$ 中二分查找 $c - add_i$ 的前驱,时间复杂度 $O(\log \sqrt{n})$。 - **碎块**:暴力遍历元素,检查 $a_j + add_i < c$ 并更新最大值,时间复杂度 $O(\sqrt{n})$。 #### 代码实现(C++) ```cpp #include <iostream> #include <vector> #include <cmath> #include <algorithm> #include <climits> using namespace std; const int MAXN = 100010; int n, block_size; vector<int> a; vector<int> add_tag; // 块的加法标记 vector<vector<int>> blocks; // 每块的有序数组 void init() { block_size = sqrt(n); int block_num = (n + block_size - 1) / block_size; add_tag.assign(block_num, 0); blocks.resize(block_num); for (int i = 0; i < n; i++) blocks[i / block_size].push_back(a[i]); for (int i = 0; i < block_num; i++) sort(blocks[i].begin(), blocks[i].end()); } void rebuild_block(int block_id) { blocks[block_id].clear(); int start = block_id * block_size; int end = min(start + block_size, n); for (int i = start; i < end; i++) blocks[block_id].push_back(a[i]); sort(blocks[block_id].begin(), blocks[block_id].end()); } void range_add(int l, int r, int c) { int block_l = l / block_size, block_r = r / block_size; if (block_l == block_r) { // 同一块内 for (int i = l; i <= r; i++) a[i] += c; rebuild_block(block_l); } else { // 左侧碎块 for (int i = l; i < (block_l + 1) * block_size; i++) a[i] += c; rebuild_block(block_l); // 中间整块 for (int i = block_l + 1; i < block_r; i++) add_tag[i] += c; // 右侧碎块 for (int i = block_r * block_size; i <= r; i++) a[i] += c; rebuild_block(block_r); } } int query_predecessor(int l, int r, int c) { int ans = -1; int block_l = l / block_size, block_r = r / block_size; // 左侧碎块暴力查询 for (int i = l; i < min(r + 1, (block_l + 1) * block_size); i++) { int val = a[i] + add_tag[block_l]; if (val < c && val > ans) ans = val; } // 中间整块二分查找 for (int i = block_l + 1; i < block_r; i++) { int target = c - add_tag[i]; auto it = lower_bound(blocks[i].begin(), blocks[i].end(), target); if (it != blocks[i].begin()) { it--; int val = *it + add_tag[i]; if (val < c && val > ans) ans = val; } } // 右侧碎块暴力查询 if (block_l != block_r) { for (int i = block_r * block_size; i <= r; i++) { int val = a[i] + add_tag[block_r]; if (val < c && val > ans) ans = val; } } return ans; } int main() { cin >> n; a.resize(n); for (int i = 0; i < n; i++) cin >> a[i]; init(); for (int i = 0; i < n; i++) { int op, l, r, c; cin >> op >> l >> r >> c; l--; r--; // 转为0-indexed if (op == 0) range_add(l, r, c); else cout << query_predecessor(l, r, c) << endl; } return 0; } ``` #### 算法分析 - **时间复杂度**: - 单次修改/查询:$O(\sqrt{n} \log \sqrt{n})$(碎块排序主导)。 - 总操作 $m$ 次:$O(m \sqrt{n} \log n)$。 - **空间复杂度**:$O(n)$。 #### 优化技巧 1. **减少排序次数**: - 碎块修改时只重构受影响块的有序数组。 2. **块大小调整**: - 实测调整块大小为 $n^{0.6}$ 可能更快(需测试)。 #### 应用场景 分块算法适用于**强制在线**的区间问题(如 LOJ 的数列分块系列题),在 $O(\sqrt{n})$ 复杂度下平衡修改与查询[^1]。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值