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

题目分析

显然不可能高斯消元。

考虑反演。

(b_i=sumlimits_{j=1}^ngcd(i,j)^Ccdot ext{lcm}(i,j)^Dcdot x_j)

(b_i=sumlimits_{j=1}^ngcd(i,j)^Ccdot frac{i^Dcdot j^D}{gcd(i,j)^D}cdot x_j)

(b_i=sumlimits_{j=1}^ngcd(i,j)^{C-D}cdot i^Dcdot j^Dcdot x_j)

实际上形如(b_i=sumlimits_{j=1}^nf(gcd(i,j))cdot g(i)cdot h(j)cdot x_j)都可以做。

我们按照套路化一下式子。

(b_i=sumlimits_{d|i}sumlimits_{d|j}[gcd(i,j)=d]cdot f(d)cdot g(i)cdot h(j)cdot x_j)

([gcd(i,j)=d])换成(sumlimits_{k|frac{gcd(i,j)}{d}}mu(k))

(b_i=sumlimits_{d|i}sumlimits_{d|j}sumlimits_{k|frac{gcd(i,j)}{d}}mu(k)cdot f(d)cdot g(i)cdot h(j)cdot x_j)

(b_i=sumlimits_{d|i}sumlimits_{d|j}sumlimits_{k cdot d|gcd(i,j)}mu(k)cdot f(d)cdot g(i)cdot h(j)cdot x_j)

(b_i=sumlimits_{T|i}sumlimits_{T|j}sumlimits_{d|T}mu(frac{T}{d})cdot f(d)cdot g(i)cdot h(j)cdot x_j)

(frac{b_i}{g(i)}=sumlimits_{T|i}sumlimits_{T|j}sumlimits_{d|T}mu(frac{T}{d})cdot f(d)cdot h(j)cdot x_j)

(fr(T)=sumlimits_{d|T}mu(frac{T}{d})cdot f(d))

(frac{b_i}{g(i)}=sumlimits_{T|i}sumlimits_{T|j}fr(T)cdot h(j)cdot x_j)

(frac{b_i}{g(i)}=sumlimits_{T|i}fr(T)sumlimits_{T|j}h(j)cdot x_j)

(q(T)=sumlimits_{T|j}h(j)cdot x_j)

(frac{b_i}{g(i)}=sumlimits_{T|i}fr(T)cdot q(T))

(fr(i)cdot q(i)=sumlimits_{T|i}mu(frac{i}{T})cdot frac{b_T}{g(T)})

所以就可以求出(q(i))了。

求出(q(i))后,再次反演,

(h(i)cdot x_i=sumlimits_{i|j}mu(frac{j}{i})cdot q(j))

那么就很容易求出(x_i)了。

注意一下无解的情况即可。

#include <bits/stdc++.h>
using namespace std;
inline int Getint(){
	int x=0,f=1;char ch=getchar();
	while(!isdigit(ch))ch!='-'?:f=-1,ch=getchar();
	while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
	return x*f;
}
typedef long long ll;
const int Maxn=100005,mod=998244353,pmod=mod-1;
int n,q,c,d,b[Maxn],ans[Maxn],h[Maxn],fr[Maxn],g[Maxn];
int mu[Maxn],Prime[Maxn];
bool vis[Maxn];
ll Pow(ll x,ll k){
	ll ret=1;
	while(k){
		if(k&1)ret=ret*x%mod;
		k>>=1;x=x*x%mod;
	}
	return ret;
}
void init(){
	mu[1]=1;
	for(int i=2;i<=100000;i++){
		if(!vis[i]){Prime[++Prime[0]]=i;mu[i]=-1;}
		for(int j=1;j<=Prime[0]&&i*Prime[j]<=100000;j++){
			vis[i*Prime[j]]=1;
			if(i%Prime[j]==0){mu[i*Prime[j]]=0;break;}
			mu[i*Prime[j]]=-mu[i];
		}
	}
	int mi=((c-d)%pmod+pmod)%pmod;
	for(int i=1;i<=100000;i++){
		int tmp=Pow(i,mi);
		for(int j=1;i*j<=100000;j++)
			fr[i*j]=(fr[i*j]+(ll)mu[j]*tmp)%mod;
	}
	for(int i=1;i<=100000;i++)g[i]=Pow(i,d);
}
void solve(){
	memset(h,0,sizeof(h));
	memset(ans,0,sizeof(ans));
	int mi=((-d)%pmod+pmod)%pmod;
	for(int i=1;i<=n;i++)b[i]=(ll)b[i]*Pow(i,mi)%mod;
	for(int i=1;i<=n;i++)
		for(int j=1;i*j<=n;j++)
			h[i*j]=(h[i*j]+(ll)mu[j]*b[i])%mod;
	for(int i=1;i<=n;i++){
		if(fr[i]==0&&h[i]!=0){puts("-1");return;}
		h[i]=(ll)h[i]*Pow(fr[i],mod-2)%mod;
	}
	for(int i=1;i<=n;i++)
		for(int j=1;i*j<=n;j++)
			ans[i]=(ans[i]+(ll)mu[j]*h[i*j])%mod;
	for(int i=1;i<=n;i++){
		if(g[i]==0&&ans[i]!=0){puts("-1");return;}
		if(g[i])ans[i]=(ll)ans[i]*Pow(g[i],mod-2)%mod;
		else ans[i]=0;
	}
	for(int i=1;i<=n;i++)cout<<(ans[i]+mod)%mod<<" 
"[i==n];
}
int main(){
	n=Getint();c=Getint();d=Getint();q=Getint();
	init();
	while(q--){
		for(int i=1;i<=n;i++)b[i]=Getint();
		solve();
	}
}
原文地址:https://www.cnblogs.com/Trrui/p/10005672.html