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

﹏ヽ暗。殇╰゛Y 2021-11-13 18:58 169阅读 0赞

[https://scut.online/p/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

[https_scut.online_p_114]: https://scut.online/p/114