【UOJ#62】【UR #5】怎样跑得更快（莫比乌斯反演）

题面

题解

众所周知，\(lcm(i,j)=\frac{ij}{gcd(i,j)}\)，于是原式就变成了：

\[\sum_{j=1}^n gcd(i,j)^{c-d}i^dj^dx_j\equiv b_i
\]

于是我们就可以写成函数的形式：

\[\sum_{j=1}^n f(gcd(i,j))h(i)h(j)x_j\equiv b_i
\]

然后就开始枚举\(gcd\)。

\[\begin{aligned}
b_i&=\sum_{d=1}^n f(d)\sum_{j=1}^n [gcd(i,j)=d]h(i)h(j)x_j\\
&=\sum_{d=1}^n f(d) \sum_{j=1}^{n}[\frac{gcd(i,j)}{d}=1]h(i)h(j)x_j\\
&=\sum_{d|i} f(d) \sum_{d|j}h(i)h(j)\sum_{k|\frac{gcd(i,j)}{d}}\mu(k)x_j
\end{aligned}\]

条件等价于\(kd|gcd(i,j)\)，令\(T=kd\)。

\[\begin{aligned}
b_i&=h(i)\sum_{d|i}f(d)\sum_{d|j}h(j)\sum_{T|gcd(i,j)}\mu(k)x_j\\
&=h(i)\sum_{T|i}\sum_{T|j}\sum_{d|T}f(d)x_j\mu(\frac{T}{d})h(j)\\
&=h(i)\sum_{T|i}\sum_{T|j}g(j)x_j\sum_{d|T}\mu(\frac{T}{d})f(d)
\end{aligned}\]

后半部分可以提前预处理出来，记做\(fr(T)\)。

继续往下就是：

\[\begin{aligned}
b_i&=h(i)\sum_{T|i}\sum_{T|j}h(j)x_j fr(T)\\
&=h(i)\sum_{T|i}fr(T)\sum_{T|j}h(j)x_j
\end{aligned}\]

考虑把后半部分的那个东西也给提前算出来，记做\(g(T)\)。

那么要求的就变成了\(\displaystyle b_i=h(i)\sum_{T|i}fr(T)g(T)\)

令\(gr(T)=fr(T)g(T)\)，于是\(\displaystyle b_i=h(i)\sum_{T|i}gr(T)\)。

根据莫比乌斯反演可以得到：

\[\frac{b_i}{h(i)}=\sum_{T|i}gr(T)
\]

\[gr(i)=\sum_{T|i}\mu(\frac{i}{T})\frac{b_T}{h(T)}
\]

而\(gr\)是可以算出来的，\(fr\)也是可以算出来的，所以\(g\)也是可以算出来的。

而\(\displaystyle g(T)=\sum_{T|j}h(j)x_j\)，那么通过莫比乌斯反演可以算出\(h(j)x_j\)，直接除掉之后就能算出\(x_j\)了。

如果无解就是存在除法的时候出现了非零数除以\(0\)。

#include<iostream>

#include<cstdio>

#include<cstdlib>

#include<cstring>

#include<cmath>

#include<algorithm>

#include<vector>

using namespace std;

#define ll long long

#define MOD 998244353

#define MAX 100100

inline int read()

{

	int x=0;bool t=false;char ch=getchar();

	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();

	if(ch=='-')t=true,ch=getchar();

	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();

	return t?-x:x;

}

int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}

int n,C,D,Q;

int pri[MAX],tot,mu[MAX];

bool zs[MAX];

void Sieve()

{

	mu[1]=1;

	for(int i=2;i<=n;++i)

	{

		if(!zs[i])pri[++tot]=i,mu[i]=-1;

		for(int j=1;i*pri[j]<=n;++j)

		{

			zs[i*pri[j]]=true;

			if(i%pri[j]==0)break;

			mu[i*pri[j]]=-mu[i];

		}

	}

}

int b[MAX],gr[MAX],h[MAX],invh[MAX],f[MAX],g[MAX],fr[MAX],w[MAX],invfr[MAX];

int main()

{

	n=read();C=read();D=read();Q=read();Sieve();

	for(int i=1;i<=n;++i)h[i]=fpow(i,D),invh[i]=fpow(h[i],MOD-2);;

	for(int i=1,p=(C-D+MOD-1)%(MOD-1);i<=n;++i)f[i]=fpow(i,p);

	for(int i=1;i<=n;++i)

		for(int j=i;j<=n;j+=i)

			fr[j]=(0ll+fr[j]+mu[j/i]*f[i]+MOD)%MOD;

	for(int i=1;i<=n;++i)invfr[i]=fpow(fr[i],MOD-2);

	while(Q--)

	{

		for(int i=1;i<=n;++i)b[i]=read();

		bool fl=true;

		for(int i=1;i<=n;++i)w[i]=1ll*b[i]*invh[i]%MOD;

		for(int i=1;i<=n;++i)if(b[i]&&!h[i])fl=false;

		for(int i=1;i<=n;++i)gr[i]=0;

		for(int i=1;i<=n;++i)

			for(int j=i;j<=n;j+=i)

				gr[j]=(0ll+gr[j]+mu[j/i]*w[i]+MOD)%MOD;

		for(int i=1;i<=n;++i)g[i]=1ll*gr[i]*invfr[i]%MOD;

		for(int i=1;i<=n;++i)if(gr[i]&&!invfr[i])fl=false;

		for(int i=1;i<=n;++i)w[i]=0;

		for(int i=1;i<=n;++i)

			for(int j=i;j<=n;j+=i)

				w[i]=(0ll+w[i]+mu[j/i]*g[j]+MOD)%MOD;

		for(int i=1;i<=n;++i)if(w[i]&&!h[i])fl=false;

		for(int i=1;i<=n;++i)w[i]=1ll*w[i]*invh[i]%MOD;

		if(!fl){puts("-1");continue;}

		for(int i=1;i<=n;++i)printf("%d ",w[i]);

		puts("");

	}

	return 0;

}

巴特西

【UOJ#62】【UR #5】怎样跑得更快（莫比乌斯反演）

【UOJ#62】【UR #5】怎样跑得更快（莫比乌斯反演）

题面

题解

最新文章

热门文章