ZJT's Blog

大変に気分がいい~

[LOJ 2476] [2018 集训队互测 Day 3] 蒜头的奖杯

2018-7-3 3:152018-7-3 8:9
集训队互测LOJ莫比乌斯反演

题目链接

这似乎是SDOI那几道题的强化版?

考虑枚举$i$,并对$(i,j),(i,k)$进行反演(本文用$(x,y)$表示$x,y$的$\gcd$): $$ \begin{aligned} &\sum_{i,j,k}A_iB_jC_kD_{(i,j)}E_{(i,k)}F_{(j,k)}\\ =&\sum_{i,j,k} A_iB_jC_kF_{(j,k)}\sum_{p|(i,j)}(\mu D)_p\sum_{q|(i,k)}(\mu E)_q\\ =&\sum_{p,q}(\mu D)_p(\mu E)_q\sum_{\mathbb{lcm}(p,q)|i}A_i\sum_{p|j}B_j\sum_{q|k}C_kF_{(j,k)}\\ &(记d=(p,q),x=p/d,y=q/d,n'=\lfloor n/d\rfloor,\\ &A'_d=\sum_{d|x,x\leq n}A_x,\\ &D'=\mu D,E'=\mu E)\\ =&\sum_d\sum_{xy\leq n',(x,y)=1}D'_{dx}E'_{dy}A'_{dxy}\sum_{dx|j}\sum_{dy|k}B_jC_kF_{(j,k)} \end{aligned} $$

假设$d,x,y$已经确定,且$x\leq y$,考虑计算后面的求和: $$ \begin{aligned} &\sum_{dx|j}\sum_{dy|k}B_jC_kF_{(j,k)}\\ =&\sum_{j=1}^{n'}\sum_{k=1}^{n'}[x|j,y|k]B_{dj}C_{dk}F_{d(j,k)}\\ =&\sum_{j=1}^{n'}\sum_{k=1}^{n'}[x|j,y|k]B'_jC'_kF'_{(j,k)}\\ =&\sum_{y|k}C'_k\sum_{d|k}(\mu F')_d\sum_{d|j}[x|j]B'_j \end{aligned} $$

我们暴力枚举所有的$d,x$。上面这坨东西只涉及到对因子求和以及对倍数求和,通过枚举质因子,我们可以每次用$O(n/d\log \log n/d)$的时间计算这坨东西对于每个$y$的值。这样复杂度就是$O(\sum_dn/d * \sqrt{n/d}\log \log n/d)=O(n^{1.5}\log \log n)$。中间那个根号是因为每次要在$\sqrt{n/d}$的范围内枚举$x$的取值。$x>y$的情况同理,求和式子里把$j$提到前面,再每次暴力枚举$y$就行了。

#include <bits/stdc++.h>
#define MAXN 100010
#define ULL unsigned long long

ULL read(){
    char c;
    while((c=getchar())<'0' || c>'9');
    ULL res=c-'0';
    while((c=getchar())>='0' && c<='9') res=res*10+c-'0';
    return res;
}

int n,n1;
int p[MAXN],np;
ULL A[MAXN],B[MAXN],C[MAXN],D[MAXN],E[MAXN],F[MAXN];
ULL B1[MAXN],C1[MAXN],F1[MAXN],S[MAXN];

void init(){
    static bool flag[MAXN];
    for(int i=2;i<=n;i++){
        if(!flag[i]) p[++np]=i;
        for(int j=1;j<=np && i*p[j]<=n;j++){
            flag[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
}

//s_x = \sum_{d|x} f_d
void s1(ULL *a,int n){
    for(int i=1;i<=np && p[i]<=n;i++)
        for(int j=p[i],k=1;j<=n;j+=p[i],k++)
            a[j]+=a[k];
}

//s_x = \sum_{x|d} f_d
void s2(ULL *a,int n){
    for(int i=1;i<=np && p[i]<=n;i++)
        for(int j=n/p[i],k=j*p[i];j>=1;j--,k-=p[i])
            a[j]+=a[k];
}

//s_x = \sum_{d|x} \mu(x/d)*f_d
void s3(ULL *a,int n){
    for(int i=1;i<=np && p[i]<=n;i++)
        for(int j=n/p[i],k=j*p[i];j>=1;j--,k-=p[i])
            a[k]-=a[j];
}

int gcd(int x,int y){
    while(y){ int t=x%y; x=y; y=t; }
    return x;
}

void gao_x(int x,ULL *B1,ULL *C1){
    memset(S,0,sizeof(S[0])*(n1+1));
    for(int i=x;i<=n1;i+=x) S[i]=B1[i];
    s2(S,n1);
    for(int i=1;i<=n1;i++) S[i]*=F1[i];
    s1(S,n1);
    for(int i=1;i<=n1;i++) S[i]*=C1[i];
    s2(S,n1);
}

ULL solve_d(int d){
    n1=n/d;
    for(int i=1;i<=n1;i++) B1[i]=B[i*d],C1[i]=C[i*d],F1[i]=F[i*d];
    s3(F1,n1);
    ULL res=0;
    for(int x=1;x*x<=n1;x++){
        gao_x(x,B1,C1);
        for(int y=x;y*x<=n1;y++)
            if(gcd(x,y)==1)
                res+=S[y]*A[d*x*y]*D[d*x]*E[d*y];
    }
    for(int y=1;y*y<=n1;y++){
        gao_x(y,C1,B1);
        for(int x=y+1;y*x<=n1;x++)
            if(gcd(x,y)==1)
                res+=S[x]*A[d*x*y]*D[d*x]*E[d*y];
    }
    return res;
}

int main(){
#ifdef DEBUG
    freopen("2476.in","r",stdin);
#endif
    n=read();
    for(int i=1;i<=n;i++) A[i]=read();
    for(int i=1;i<=n;i++) B[i]=read();
    for(int i=1;i<=n;i++) C[i]=read();
    for(int i=1;i<=n;i++) D[i]=read();
    for(int i=1;i<=n;i++) E[i]=read();
    for(int i=1;i<=n;i++) F[i]=read();
    init();
    s2(A,n);
    s3(D,n); s3(E,n);
    ULL ans=0;
    for(int d=1;d<=n;d++)
        ans+=solve_d(d);
    printf("%llu\n",ans);
    return 0;
}

查看详细内容

[自选题 #148] Simple Summation Problem

2017-12-26 18:392017-12-26 18:39
集训队作业自选题莫比乌斯反演

传送门

反演: $$ \begin{aligned} \sum_i F(i)&=\sum_i\sum_{d|i}F(d)[\frac id=1]\\ &=\sum_i\sum_{d|i}F(d)\sum_{p|\frac id}\mu(p)\\ &=\sum_i\sum_{d|i}F(d)\sum_{d|p|i}\mu(\frac pd)\\ &=\sum_p(\sum_{p|i}1)(\sum_{d|p}F(d)\mu(\frac pd))\\ &=\sum_p\lfloor\frac np\rfloor(\sum_{d|p}F(d)\mu(\frac pd)) \end{aligned} $$

把后面那坨东西拿出来:$G(n)=\sum_{d|n}F(d)\mu(\frac nd)$,由于积性函数的狄利克雷卷积仍然是积性函数,我们考虑$G(p^d)$($p$为质数),发现$d=1$时这个东西必定是$0$。

所以,只有每个质因子次数都大于$1$的数才会产生贡献,而些数的个数是$O(\sqrt n)$级别的(证明大概是每个这样的数都能表示成$x^2y^3$的形式,然后拿积分去近似这个复杂度,就是$\int_0^{\sqrt n}{(\frac n{x^2})}^{1/3}=O(\sqrt n)$)。所以把这些数全部搜出来,然后挨个算贡献就行了。

#include <cstdio>
#include <cstring>
#define MAXN 10000010
#define LL long long

LL n;
bool flag[MAXN];
int prime[MAXN],nump;
LL ans=0;
int debug=0;

void init(){
    for(LL i=2;i*i<=n;i++){
        if(!flag[i]) prime[++nump]=i;
        for(int j=1;j<=nump && i*i*prime[j]*prime[j]<=n;j++){
            flag[i*prime[j]]=1;
            if(i%prime[j]==0) break;
        }
    }
}

void dfs(int id,LL x,LL f){
    LL p=prime[id];
    LL f1=1,f2=p,f3=p*p;
    if(id>nump || x>n/f3){
        ans+=n/x*f;
        return;
    }
    dfs(id+1,x,f);
    LL d=n/x/f3;
    for(int i=2;d;i++){
        d/=p;
        LL temp=(i%p?f2:f3)-((i-1)%p?f1:f2);
        dfs(id+1,x*f3,f*temp);
        f1=f2; f2=f3;
        f3*=p;
    }
}

int main(){
#ifdef DEBUG
    freopen("148.in","r",stdin);
#endif
    scanf("%lld",&n);
    init();
    dfs(1,1,1);
    printf("%lld\n",ans);
}

查看详细内容