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

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

题面

UOJ

题解

众所周知，\(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;
}

转载于:https://www.cnblogs.com/cjyyb/p/11074280.html