标签:现在 没有 质数 mes typedef 简化 简洁 oid 公约数
大力水手问禅师:“大师,我觉得我光有力气是不够的。比如我吃菠菜可以让力气更大,但是却没有提升跑步的速度。请问怎样才能跑得更快?我试过吃白菜,没有效果。”
禅师浅笑,答:“方法很简单,不过若想我教你,你先看看这道\(UOJ\) \(Round\)的\(C\)题。”
令 \(p=998244353\)(\(7×17×223+17×17×223+1\),一个质数。
给你整数 \(n,c,d\)。现在有整数\(x_1,…,x_n\)和\(b_1,…,b_n\)满足\(0≤x_1,…,x_n,b_1,…,b_n<p\),且对于 \(1≤i≤n\) 满足:
\[\sum_{j=1}^n gcd(i,j)^c*lcm(i,j)^d*x_j \equiv b_i(mod \ \ p) \]
其中 \(v≡u(mod\ \ p)\) 表示 \(v\) 和 \(u\) 除以 \(p\) 的余数相等。\(gcd(i,j)\) 表示 \(i\) 和 \(j\) 的最大公约数,\(lcm(i,j)\) 表示 \(i\) 和 \(j\) 的最小公倍数。
有 \(q\) 个询问,每次给出 \(b_1,…,b_n\),请你解出 \(x_1,…,x_n\) 的值。
大力莫比乌斯反演……
直接上公式(为了方便书写,我们将\(mod \ \ p\)省略):
\[sum_{j= 1 } ^ n gcd(i, j) ^ c * lcm(i, j) ^ d * x_j = b_i\]
\[\sum_{j = 1} ^ n gcd(i, j) ^ c * \frac{i ^ d * j ^ d}{gcd(i, j) ^ d} * x_j = b_i \]
我们记\(f(n) = n^{c - d},\ \ h(n) = n ^ d\)。
\[\sum_{j=1}^n f(gcd(i,j))*h(j)*x_j = \frac{b_i}{h(i)}\]
这个地方我们发现\(gcd(i,j)^{c-d}\)比较的难以处理,但是假设我们知道了一个函数\(f_r(n)\),满足\(f(n) = \sum_{d | n}f_r(d)\),那么在知道了\(f(n)\)之后,求\(f_r(n)\)是很好求的,只需要一次莫比乌斯反演就行了。
于是原式变成了这样:
\[\sum_{j = 1} ^ n \sum_{d | gcd(i,j)} f_r(d) * h(j) * x_j \equiv \frac{b_i}{h(i)} (mod \ \ p)\]
为了简洁一点,我们记\(B(i) = \frac{b_i}{h(i)}\)。
\[\sum_{j = 1} ^ n \sum_{d | gcd(i,j)} f_r(d) * h(j) * x_j = B(i)\]
移一下项
\[\sum_{d | i} f_r(d) \sum_{d | j} h(j) * x_j = B(i) \]
考虑到\(\sum_{d | j} h(j) * x_j\)这一项是求所有是\(d\)的倍数的\(h(j)*x_j\)的值,这个值只和\(d\)有关。所以我们可以记\(A(d) = \sum_{d | j} h(j) * x_j\)
所以原式继续简化:
\[\sum_{d | i} f_r(d) * A(d) = B(i)\]
这样其实就可以进行第二次莫比乌斯反演了……为;为了看起来更清楚一些,我们再记\(F(d) = f_r(d) * A(d)\)
那么这个莫比乌斯反演就很清楚了:
\[\sum_{d | i} F(d) = B(i) \]
反演一下可得:
\[F(i) = \sum_{d | i}\mu(d) B(\frac{i}{d})\]
然后我们将\(F(n)\)展开可得:\(f_r(d) * A(d) = F(n)\) 即\(A(d) = \frac{F(n)}{f_r(d)}\)
继续讲\(A(n)\)展开可得:\(\sum_{n | j} h(j) * x_j = A(n)\),再来一次莫比乌斯反演,可以得到:\(h(j) * x_j = \sum_{n | j} \mu(\frac{j}{n})* A(j)\)。这样就可以轻松求得\(x_j\)的值,无解多解什么的中间除法的时候判一下就行了。
#pragma GCC optimize(2)
#include<bits/stdc++.h>
using namespace std;
const int Md=998244353;
const int N=3e5+500;
typedef long long ll;
inline int Add(const int &x,const int &y) {return (x+y)>=Md?(x+y-Md):(x+y);}
inline int Sub(const int &x,const int &y) {return (x-y)<0?(x-y+Md):(x-y);}
inline int Mul(const int &x,const int &y) {return (ll)x*y%Md;}
inline int Chg(int x) {int X=x;while(X<0) X+=Md;while(X>=Md) X-=Md;return X;}
int Powe(int x,int y) {
int ans=1;
while(y) {
if(y&1) ans=Mul(ans,x);
x=Mul(x,x);
y>>=1;
}
return ans;
}
inline int Div(const int &x,const int &y) {return Mul(x,Powe(y,Md-2));}
void read(int &x) {
x=0;int f=1;char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
x*=f;
}
void read(ll &x) {
x=0;int f=1;char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
x*=f;
}
int n,q,c,d,cnt;
int b[N],F[N],mu[N],pri[N],ntp[N],f[N],A[N],g[N],Iv1[N],Iv2[N],H[N];
void Init() {
mu[1]=1;
for(int i=2;i<N;i++) {
if(!ntp[i]) mu[i]=-1,pri[++cnt]=i;
for(int j=1;j<=cnt&&i*pri[j]<N;j++) {
ntp[i*pri[j]]=1;
if(i%pri[j]==0) {
mu[i*pri[j]]=0;
break;
}
else mu[i*pri[j]]=-mu[i];
}
}
}
int main() {
Init();
read(n);read(c);read(d);read(q);
for(int i=1;i<=n;i++) {
f[i]=Powe(i,c>=d?(c-d):(c-d+Md-1)%(Md-1));
g[i]=Powe(i,d);
Iv1[i]=Powe(g[i],Md-2);
}
for(int i=1;i<=n;i++) {
for(int j=i*2;j<=n;j+=i) {
f[j]=Sub(f[j],f[i]);
}
}
for(int i=1;i<=n;i++) Iv2[i]=Powe(f[i],Md-2);
while(q--) {
for(int i=1;i<=n;i++) read(b[i]),b[i]=Mul(b[i],Iv1[i]);
for(int i=1;i<=n;i++) F[i]=H[i]=A[i]=0;
for(int i=1;i<=n;i++) {
for(int j=i;j<=n;j+=i) {
F[j]=Add(F[j],Mul(Chg(mu[j/i]),b[i]));
}
}
bool fl=1;
for(int i=1;i<=n;i++) {
if(Iv2[i]) A[i]=Mul(F[i],Iv2[i]);
else if(F[i]) {
fl=0;
puts("-1");
break;
}
}
if(!fl) continue;
for(int i=1;i<=n;i++) {
for(int j=i;j<=n;j+=i) {
H[i]=Add(H[i],Mul(Chg(mu[j/i]),A[j]));
}
}
for(int i=1;i<=n;i++) {
if(g[i]==0) {
if(H[i]==0) b[i]=1;
else {
fl=0;
puts("-1");
break;
}
}
else b[i]=Mul(H[i],Iv1[i]);
}
if(fl) {
for(int i=1;i<=n;i++) {
printf("%d%c",b[i]," \n"[i==n]);
}
}
}
return 0;
}
标签:现在 没有 质数 mes typedef 简化 简洁 oid 公约数
原文地址:https://www.cnblogs.com/Apocrypha/p/10472238.html