SCUT - 114 - 作业之数学篇 - 杜教筛

https://scut.online/p/114

\(A(n)=\sum\limits_{i=1}^{n} \frac{lcm(i,n)}{gcd(i,n)}\)
\(=\sum\limits_{i=1}^{n} \frac{in}{gcd^2(i,n)}\)

枚举g：
\(A(n)=n\sum\limits_{g|n}\frac{1}{g^2} \sum\limits_{i=1}^{n} i [gcd(i,n)==g]\)

最内层除以g：
\(A(n)=n\sum\limits_{g|n}\frac{1}{g} \sum\limits_{i=1}^{\frac{n}{g}} i [gcd(i,\frac{n}{g})==1]\)

考虑里面的：
\(B(n)=\sum\limits_{i=1}^{n} i [gcd(i,n)==1]\)

显然：
\(B(n)=\frac{1}{2}n([n==1]+\varphi(n))\)

代回 \(A(n)\) 里面：
\(A(n)=\sum\limits_{g|n}\frac{n}{g} B(\frac{n}{g})\)
\(=\sum\limits_{g|n} g B(g)\)
\(=\frac{1}{2}(1+\sum\limits_{g|n} g^2 \varphi(g))\)

代回 \(S(n)\) 里面：
$S(n)=\sum\limits_{i=1}^{n} A(i) $
\(= \sum\limits_{i=1}^{n} \frac{1}{2}(1+\sum\limits_{g|i} g^2 \varphi(g))\)
\(= \frac{n}{2}+\frac{1}{2}\sum\limits_{i=1}^{n}\sum\limits_{g|i} g^2 \varphi(g)\)

记：
\(C(n)=\sum\limits_{i=1}^{n}\sum\limits_{g|i} g^2 \varphi(g)\)

显然：
\(C(n)=\sum\limits_{d=1}^{n}d^2 \varphi(d) \lfloor\frac{n}{d}\rfloor\)

后面是一个分块，前面是个杜教筛。然后我们召唤鹏哥就可以通过这道题了。

转化为快速求 $\sum\limits_{i=1}^{n} i^2 \varphi(i) $，卷一个id平方然后查鹏哥的cheatsheet过了。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

const int mod=1e9+7;
const int inv2=(mod+1)>>1;
const int MAXN=4e6;

int pri[MAXN+1];
int &pritop=pri[0];
ll B[MAXN+1];
ll PrefixB[MAXN+1];
int pk[MAXN+1];

void sieve(int n=MAXN) {
    pk[1]=1;
    B[1]=1;
    for(int i=2; i<=n; i++) {
        if(!pri[i]) {
            pri[++pritop]=i;
            pk[i]=i;
            ll tmp=1ll*i*i;
            //.
            if(tmp>=mod)
                tmp%=mod;
            B[i]=tmp*(i-1);
            if(B[i]>=mod)
                B[i]%=mod;
            //.
        }
        for(int j=1; j<=pritop; j++) {
            int &p=pri[j];
            int t=i*p;
            //.
            if(t>n)
                break;
            pri[t]=1;
            if(i%p) {
                pk[t]=p;
                B[t]=B[i]*B[p];
                if(B[t]>=mod)
                    B[t]%=mod;
                //.
            } else {
                pk[t]=pk[i]*p;
                if(pk[t]==t) {
                    ll tmp=1ll*t*t;
                    if(tmp>=mod)
                        tmp%=mod;
                    B[t]=tmp*(t-i);

                    if(B[t]>=mod)
                        B[t]%=mod;
                    //.
                } else {
                    B[t]=B[t/pk[t]]*B[pk[t]];
                    if(B[t]>=mod)
                        B[t]%=mod;
                }
                break;
            }
        }
    }
    for(int i=1; i<=n; i++) {
        PrefixB[i]=PrefixB[i-1]+B[i];
        if(PrefixB[i]>=mod)
            PrefixB[i]-=mod;
    }
}

inline int phi(int n) {
    int res=n;
    for(int i=1; i<=pritop&&pri[i]*pri[i]<=n; i++) {
        if(n%pri[i]==0) {
            res-=res/pri[i];
            while(n%pri[i]==0)
                n/=pri[i];
        }
        if(n==1)
            return res;
    }
    res-=res/n;
    return res;
    //.
}

inline ll qpow(ll x,int n) {
    ll res=1;
    while(n) {
        if(n&1) {
            res*=x;
            if(res>=mod)
                res%=mod;
        }
        x*=x;
        if(x>=mod)
            x%=mod;
        n>>=1;
    }
    return res;
}

int inv4=qpow(4,mod-2);
int inv6=qpow(6,mod-2);

inline ll s2(ll n) {
    ll tmp=n*(n+1);
    if(tmp>=mod)
        tmp%=mod;
    tmp*=(n*2+1);
    if(tmp>=mod)
        tmp%=mod;
    tmp*=inv6;
    if(tmp>=mod)
        tmp%=mod;
    return tmp;
}

inline ll s3(ll n) {
    ll tmp=n*(n+1);
    if(tmp>=mod)
        tmp%=mod;
    tmp*=tmp;
    if(tmp>=mod)
        tmp%=mod;
    tmp*=inv4;
    if(tmp>=mod)
        tmp%=mod;
    return tmp;
}

unordered_map<int,ll> PB;

inline ll S(int n) {
    if(n<=MAXN)
        return PrefixB[n];
    if(PB.count(n))
        return PB[n];
    ll ret=s3(n);
    for(int l=2,r; l<=n; l=r+1) {
        int t=n/l;
        r=n/t;
        ll tmp=s2(r)-s2(l-1);
        if(tmp<0)
            tmp+=mod;
        tmp*=S(t);
        if(tmp>=mod)
            tmp%=mod;
        ret-=tmp;
        if(ret<0)
            ret+=mod;
    }
    return PB[n]=ret;
}

inline ll Prefix(int l,int r) {
    l--;
    ll PL=S(l);
    ll PR=S(r);
    ll res=PR-PL;
    if(res<0)
        res+=mod;
    return res;
}

inline ll P(int n) {
    ll res=0;
    for(int l=1,r; l<=n; l=r+1) {
        int t=n/l;
        r=n/t;
        ll tmp=1ll*t*Prefix(l,r);
        if(tmp>=mod)
            tmp%=mod;
        res+=tmp;
        if(res>=mod)
            res-=mod;
    }
    return res;
}

inline ll Ans(int n) {
    ll res=(1ll*n*inv2)%mod;
    res+=(P(n)*inv2)%mod;
    return res%mod;
}

int main() {
#ifdef Yinku
    freopen("Yinku.in","r",stdin);
#endif // Yinku
    sieve();
    int n;
    while(~scanf("%d",&n)) {
        printf("%lld\n",Ans(n));
    }
    return 0;
}

转载于:https://www.cnblogs.com/Yinku/p/11001111.html