首页 > 其他 > 详细

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

时间:2019-06-23 21:24:09      阅读:202      评论:0      收藏:0      [点我收藏+]

【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;
}

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

原文:https://www.cnblogs.com/cjyyb/p/11074280.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!