51nod 1237 最大公约数之和 V3(杜教筛)

forever97 posted @ 2016年10月28日 12:53 in 数学-杜教筛 with tags 杜教筛 , 304 阅读

 

【题目链接】 https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1237

 

【题目大意】

    求[1,n][1,n]最大公约数之和

 

【题解】

    枚举最大公约数k,得到答案为2*∑(k*phi_sum(n/k))-n*(n+1)/2

    phi_sum可以利用杜教筛实现

 

【代码】

#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long LL;
const int mod=1333331,inv2=500000004;
const LL MOD=1e9+7;
LL a,b,n,miu[5000010],phi[5000010];
int p[500010],cnt=0,i,tot;
bool v[5000010];
struct HASHMAP{
    int h[mod+10],cnt,nxt[100010];
    LL st[100010],S[100010];
    void push(LL k,LL v){
        int key=k%mod;
        for(int i=h[key];i;i=nxt[i]){
            if(S[i]==k)return;
        }++cnt;nxt[cnt]=h[key];h[key]=cnt;
        S[cnt]=k;st[cnt]=v;
    }
    LL ask(LL k){
        int key=k%mod;
        for(int i=h[key];i;i=nxt[i]){
            if(S[i]==k)return st[i];
        }return -1;
    }
}H;
void Get_Prime(){
    for(miu[1]=phi[1]=1,i=2;i<=5000000;i++){
        if(!v[i])p[tot++]=i,miu[i]=-1,phi[i]=i-1;
        for(int j=0;i*p[j]<=5000000&&j<tot;j++){
            v[i*p[j]]=1;
            if(i%p[j]){
                miu[i*p[j]]=-miu[i];
                phi[i*p[j]]=phi[i]*(p[j]-1);
            }else{
                miu[i*p[j]]=0;
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
        }
    }for(int i=2;i<=5000000;++i)phi[i]=(phi[i-1]+phi[i])%MOD;
}
LL phi_sum(LL n){
    if(n<=5000000)return phi[n];
    LL tmp=H.ask(n),la,A=0;
    if(tmp!=-1)return tmp;
    for(LL i=2;i<=n;i=la+1){
        LL now=n/i; la=n/now;
        (A+=(la-i+1)%MOD*phi_sum(n/i)%MOD)%=MOD;
    }A=((n%MOD)*(n%MOD+1)%MOD*inv2%MOD-A+MOD)%MOD;
	H.push(n,A);return A;
}
int main(){
    scanf("%lld",&n);
    Get_Prime();
    LL la,ans=0;
    for(LL i=1;i<=n;i=la+1){
        LL now=n/i;
        la=n/now;
        ans=(ans+(i+la)%MOD*(la-i+1)%MOD*inv2%MOD*phi_sum(now)%MOD)%MOD;
    }ans=ans*2%MOD; LL k=n%MOD;
    ans=(ans-k*(k+1)%MOD*inv2%MOD+MOD)%MOD;
    printf("%lld\n",ans);
    return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter