• SCUT


    https://scut.online/p/114

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

    枚举g:
    (A(n)=nsumlimits_{g|n}frac{1}{g^2} sumlimits_{i=1}^{n} i [gcd(i,n)==g])

    最内层除以g:
    (A(n)=nsumlimits_{g|n}frac{1}{g} sumlimits_{i=1}^{frac{n}{g}} i [gcd(i,frac{n}{g})==1])

    考虑里面的:
    (B(n)=sumlimits_{i=1}^{n} i [gcd(i,n)==1])

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

    代回 (A(n)) 里面:
    (A(n)=sumlimits_{g|n}frac{n}{g} B(frac{n}{g}))
    (=sumlimits_{g|n} g B(g))
    (=frac{1}{2}(1+sumlimits_{g|n} g^2 varphi(g)))

    代回 (S(n)) 里面:
    $S(n)=sumlimits_{i=1}^{n} A(i) ( )= sumlimits_{i=1}^{n} frac{1}{2}(1+sumlimits_{g|i} g^2 varphi(g))( )= frac{n}{2}+frac{1}{2}sumlimits_{i=1}^{n}sumlimits_{g|i} g^2 varphi(g)$

    记:
    (C(n)=sumlimits_{i=1}^{n}sumlimits_{g|i} g^2 varphi(g))

    显然:
    (C(n)=sumlimits_{d=1}^{n}d^2 varphi(d) lfloorfrac{n}{d} floor)

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

    转化为快速求 $sumlimits_{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
    ",Ans(n));
        }
        return 0;
    }
    
  • 相关阅读:
    解决Excel vba case过程中遇到的代码难题
    Excel VBA——数据类型
    Outlook与Exchange的概念
    使用Aouth2进行身份验证
    OutLook开发人员文档学习
    SharePoint网站预热(Warm Up),提升响应速度
    SharePoint 2013下,使用ajax调用ashx报Http 302错误
    【SharePoint】图说搜索服务Search Service详细配置(一)
    NLog 日志框架搭建讲解(亲测有效,代码齐全)
    RabbitMQ 教程(四)RabbitMQ并发处理
  • 原文地址:https://www.cnblogs.com/Yinku/p/11001111.html
Copyright © 2020-2023  润新知