UOJ#62【UR #5】怎样跑得更快(反演)

UOJ#62【UR #5】怎样跑得更快

题目描述

大力水手问禅师：“大师，我觉得我光有力气是不够的。比如我吃菠菜可以让力气更大，但是却没有提升跑步的速度。请问怎样才能跑得更快？我试过吃白菜，没有效果。”
禅师浅笑，答：“方法很简单，不过若想我教你，你先看看这道UOJ Round的C题。”
令 (p = 998244353)（(7 imes 17 imes 2^{23} + 1)，一个质数）。
给你整数 (n, c, d)。现在有整数 (x_1, dots, x_n) 和 (b_1, dots, b_n) 满足 (0 leq x_1, dots, x_n, b_1, dots, b_n lt p)，且对于 (1 leq i leq n) 满足：

[egin{equation} sum_{j = 1}^{n} gcd(i, j)^c cdot lcm(i, j)^d cdot x_j equiv b_i pmod{p} end{equation} ]

其中 (v equiv u pmod{p}) 表示 (v) 和 (u) 除以 (p) 的余数相等。(gcd(i, j)) 表示 (i) 和 (j) 的最大公约数，(lcm(i, j)) 表示 (i) 和 (j) 的最小公倍数。
有 (q) 个询问，每次给出 (b_1, dots, b_n)，请你解出 (x_1, dots, x_n) 的值。

数据范围

[1 le n le 10^5, ~~~~q le 3 ]

解题思路

好题 + 1

先说一个有启发的 (Theta(n^3+qn^2)) 做法

直接将系数暴力算出直接跑高斯消元是 (Theta(n^3q)) 的

有一个技巧，因为系数矩阵相同，所以对系数矩阵求逆即可

[A * V = Ans\ A^T*A*V = A^T*Ans\ V =A^T*Ans ]

向量乘矩阵是 (Theta(n^2)) 的

正解（参考vfk解法）

转化一波式子，变为

[sum_{j = 1}^{n} gcd(i, j)^c cdot lcm(i, j)^d cdot x_j equiv b_i pmod{p}\ sum_{j = 1}^{n} gcd(i, j)^{c-d}cdot i^d*j^d cdot x_j equiv b_i pmod{p}\ ]

类似于

[sum_{j = 1}^{n} F(gcd(i, j))cdot h(i)*g(j) cdot x_j equiv b_i pmod{p}\ ]

的式子都可以做

同时设 (F(x) = sum_{i|d}f(i))，有

[f(x) = sum_{d|x}F(d)*mu(frac xd)\ f(x) = F(x) - sum_{d|x且d != x} f(d)\ ]

两种做法均可在 (Theta(NlogN)) 的时间复杂度内求出

将 F 换为 f 得到

[sum_{j =1}^n sum_{d}[d|i][d|j]f(d)cdot h(i)*g(j) cdot x_j equiv b_i pmod{p}\ sum_{d|i}f(d)sum_{j =1}^n [d|j]cdot g(j) cdot x_j equiv b_i*h(i)^{-1} pmod{p}\ let ~~H_i = b_icdot h(i)^{-1},~~G_i= sum_{i|j}g(j) cdot x_j\ sum_{d|i}f(d)G_d=H_i\ ]

这个东西，很熟悉，和刚才设 (F(x) = sum_{i|d}f(i)) 一样可以递推过去，然后再对 (G_i) 进行一次递推

总复杂度 (Theta(qNlogN))

代码

#pragma GCC optimize(3, "inline")
#include <iostream>
#include <cstdio>
#include <cstring>
#define ll long long
using namespace std;

template <typename T>
void read(T &x) {
    x = 0; bool f = 0;
    char c = getchar();
    for (;!isdigit(c);c=getchar()) if (c=='-') f=1;
    for (;isdigit(c);c=getchar()) x=x*10+(c^48);
    if (f) x=-x;
}

const int P = 998244353;
int n, c, d, q;

ll fpw(ll x, ll mi) {
	ll res = 1;
	while (mi) {
		if (mi & 1) res = res * x % P;
		x = x * x % P, mi >>= 1;
	}
	return res;
}

const int N = 200500;
ll f[N], h[N], G[N], b[N], X[N], B[N];
void getf(void) {
	for (int i = 1;i <= n; i++) {
		f[i] += fpw(i, (P - 1 + c - d) % (P - 1));
		for (int j = i + i;j <= n; j += i) f[j] -= f[i];
		f[i] = (f[i] % P + P) % P, f[i] = fpw(f[i], P - 2);
	}
}

int main() {
	read(n), read(c), read(d), read(q);
	c %= (P - 1), d %= (P - 1), getf();
	for (int i = 1;i <= n; i++) 
		h[i] = fpw(i, P - 1 - d);
	while (q--) {
		for (int i = 1;i <= n; i++) 
			read(G[i]), G[i] = G[i] * h[i] % P;
		bool fl = 0;
		for (int i = 1;i <= n; i++) {
			for (int j = i + i;j <= n; j += i)
				G[j] -= G[i];
			G[i] = (G[i] % P + P) % P;
			if (G[i] && !f[i]) { fl = 1; break; }
			G[i] = G[i] * f[i] % P;
		}
		if (fl) { puts("-1"); continue; }
		for (int i = n;i >= 1; i--) {
			for (int j = i + i;j <= n; j += i)
				G[i] -= G[j];
			G[i] = (G[i] % P + P) % P;
		}
		for (int i = 1;i <= n; i++) printf ("%lld ", G[i] * h[i] % P);
		puts("");
	}
	return 0;
}

给个最短代码?

#include <bits/stdc++.h>
#define ll long long
#define rep() for(int i=1;i<=n;i++)
#define Rep() for(int j=i+i;j<=n;j+=i)
const int P=998244353;
const int N=200500;
ll f[N],h[N],G[N];
int n,c,d,q;
ll w(ll x,ll m) {
	ll r=1;
	for (;m;m>>=1,x=x*x%P)(m&1)&&(r=r*x%P);
	return r;
}
int main() {
	scanf ("%d%d%d%d",&n,&c,&d,&q);c%=(P-1),d%=(P-1);
	rep(){
		h[i]=w(i,P-1-d),f[i]+=w(i,(P-1+c-d)%(P-1));
		Rep()f[j]-=f[i];f[i]=(f[i]%P+P)%P,f[i]=w(f[i],P-2);
	}
	while (q--) {
		rep() scanf("%lld",G+i),G[i]=G[i]*h[i]%P;ll fl=0;
		rep(){
			Rep() G[j]-=G[i];G[i]=(G[i]%P+P)%P;
			if(G[i]&&!f[i]){fl=1;break;}G[i]=G[i]*f[i]%P;}
		if(fl){puts("-1");continue;}
		for (int i=n;i>=1;G[i]=(G[i]%P+P)%P,i--) Rep()G[i]-=G[j]; 
		rep()printf("%lld ",G[i]*h[i] % P);puts("");
	}
	return 0;
}