hdu5341 Gcd and Lcm

最新推荐文章于 2018-09-05 17:41:24 发布

lshacm

最新推荐文章于 2018-09-05 17:41:24 发布

阅读量539

点赞数

分类专栏： acm_数论

本文链接：https://blog.csdn.net/u013514182/article/details/47208009

版权

acm_数论专栏收录该内容

21 篇文章 0 订阅

订阅专栏

题意：求 $\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^n\sum_{l=1}^nlcm[gcd(i,j),gcd(k,l)]$ ，n<=1e7,数据不超过100组
解法：直接求n^4一定会超时，因此需要逐步化简，直到可以求解为止。
1.首先求F[n,d]= $\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)==d$ ,
令F0[n]= $\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)==1$ ，这个式子可以化简为 $\sum_{d=1}^nu(d)*(n/d)*(n/d)$ ，u(d)是莫比乌斯函数，可以由线性筛得到，预处理它的前n项和，F0[n]可以在O(sqrt(n))复杂度内分块得到。F[n,d]=F0[n/d]，所以可以在O(sqrt(n))内求出F[n,d]

2.合并相同gcd项，上式= $\sum_{d1=1}^n\sum_{d2=1}^nF[n,d1]*F[n,d2]*lcm[d1,d2]$
= $\sum_{d1=1}^n\sum_{d2=1}^nF0[n/d1]*F0[n/d2]*lcm[d1,d2]$
= $\sum_{d1=1}^n\sum_{d2=1}^nF0[n/d1]*F0[n/d2]*d1*d2/(d1,d2)$
令(d1,d2)=d,d1=d*a,d2=d*b
= $\sum_{d=1}^n\sum_{a=1}^{n/d}\sum_{b=1}^{n/d}F0[n/d/a]*F0[n/d/b]*a*b*d*(gcd(a,b)==1)$
令g[i]= $\sum_{j=1}^{i}j*F0[i/j]$
= $\sum_{d=1}^nd*\sum_{p=1}^{n/d}u[p]*p*p*(g[n/d/p])^2$ (这一步参考贾志鹏论文中推导lcm之和的过程)
令T=d*p
= $\sum_{T=1}^n(g[n/T])^2*(T*\sum_{p|T}u[p]*p)$
令S[T]= $T*\sum_{p|T}u[p]*p$ ，为积性函数，可以预处理出S[T]前n项和
由于n/T最多有sqrt(n)个不同取值，所以可以用分块求ans
g[i]= $\sum_{j=1}^{i}j*F0[i/j]$ ，i/j最多有sqrt(i)个不同取值，所以可以用分块求g[i]
又F0[n]也可以用分块求解，并且对于每次的n，最多有sqrt(n)个不同的f取值，所以相当于对于每个询问n，先把需要求的f值预处理出来，这样就少一层sqrt(n)，总复杂度为O(n)

注：代码需要向右拖动才可以看到完整的，不要漏掉了sumu0的求和。代码丑陋了，不好意思~希望朋友们不要介意~

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<vector>
#include<cmath>
#include<string.h>
#define ll long long
using namespace std;

const ll maxn = 10000000+10;
const ll mod = (1ll<<32);
ll t,n,ans;
ll f0[maxn],g0[maxn];
ll sumhe[maxn],sumu[maxn],sumu0[maxn];
bool not_prime[maxn];
ll pri[maxn],cnt;
ll ans0,ans1;

void mobius(){//筛法求预处理莫比乌斯函数和几个要用到的数组
    sumu[1]=1,sumu0[1]=1,sumu[0]=0,sumhe[0]=0,sumu0[0]=0,cnt=0;
    for(ll i=2;i<=maxn-5;i++){
        if(!not_prime[i]){ 
            pri[cnt++]=i;
            sumu[i]=((-1)%mod+mod)%mod;
            sumu0[i]=((1-i)%mod+mod)%mod;
        }
        for(ll j=0;j<cnt&&i*pri[j]<=maxn-5;j++){
            not_prime[i*pri[j]]=1;
            if(i%pri[j]==0) { sumu[i*pri[j]]=0,sumu0[i*pri[j]]=sumu0[i]%mod;break; }
            else sumu[i*pri[j]]=((-sumu[i])%mod+mod)%mod,sumu0[i*pri[j]]=(sumu0[pri[j]]*sumu0[i])%mod;
        }
    }
    for(ll i=1;i<=maxn-5;i++) sumu[i]=(sumu[i-1]+sumu[i])%mod,f0[i]=g0[i]=mod+1,sumhe[i]=(sumhe[i-1]+i)%mod,sumu0[i]=(sumu0[i-1]+i*sumu0[i]%mod)%mod;
}

ll calf(ll n){
    if(f0[n]!=mod+1) return f0[n];
    ll ans=0,tmp;
    for(ll i=1;i<=n;i=tmp+1){
        tmp=n/(n/i);
        ans=(ans+((sumu[tmp]-sumu[i-1])%mod+mod)%mod*(n/i)%mod*(n/i)%mod)%mod;
    }
    return f0[n]=ans;
}

ll calg(ll n){//sqrt(n)计算g0[x]
    if(g0[n]!=mod+1) return g0[n];
    ll ans=0,tmp;
    for(ll i=1;i<=n;i=tmp+1){
        tmp=n/(n/i);
        ans=(ans+((sumhe[tmp]-sumhe[i-1])%mod+mod)%mod*calf(n/i)%mod)%mod;
    }
    return g0[n]=ans;
}

int main(){
    mobius();
    //freopen("a.txt","r",stdin);
    scanf("%lld",&t);
    while(t--){
        scanf("%lld",&n);
        ans=0;
        ll tmp=0;
        for(ll i=1;i<=n;i=tmp+1){
            tmp=n/(n/i);
            ll x=calg(n/i);
            ans=(ans+((sumu0[tmp]-sumu0[i-1])%mod+mod)%mod*x%mod*x%mod)%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}