51nod 1238 最小公倍数之和 V3

51nod 1238 最小公倍数之和 V3

求
\[ \sum_{i=1}^N\sum_{j=1}^N lcm(i,j) \]
\(N\leq 10^{10}\)

先按照套路推一波反演的式子：
\[ Ans=\sum_{g=1}g\sum_{i=1}^{\frac{n}{g}}\sum_{j=1}^{\frac{n}{g}}ij\sum_{d|i,d|j}\mu(d)\\ =\sum_{g=1}g\sum_{d=1}^{\frac{n}{g}}d^2\mu(d)S^2(\frac{n}{dg})\\ =\sum_{T=1}^n\sum_{d|T}d^2\mu(d)\frac{T}{d}S^2(\frac{n}{T})\\ \]
难点在于求下面的函数的前缀和。
\[ G(n)=\sum_{d|T}d^2\mu(d)\frac{T}{d} \]
设：
\[ A(n)=n^2\mu(n)\\ B(n)=n \]
则：
\[ G(n)=A*B \]
其中\(*\)表示狄利克雷卷积。

考虑用杜教筛，也就是构造一个函数\(C(n)\)，使得\(G*C\)有些美妙的性质。

考虑从\(A(n)=n^2\mu(n)\)下手，将\(n^2\)消掉，只留下\(\mu(n)\)。

设\(C(n)=n^2\),
\[ A*C=\sum_{d|n}d^2\mu(d)(\frac{n}{d})^2\\ =n^2\sum_{d|n}\mu(d)\\ =[n==1] \]
所以：
\[ \begin{align} G*C&=(A*C)*B\\ &=\epsilon *B\\ &=\sum_{d=1}[d==1]\frac{n}{d}\\ &=n \end{align} \]
然后就是杜教筛的套路：
\[ \sum_{i=1}^n\sum_{d|n}G(n)(\frac{n}{d})^2=\sum_{i=1}^ni\\ \Rightarrow \sum_{i=1}^ni^2\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}G(j)=\frac{n*(n+1)}{2}\\ \Rightarrow \sum_{i=1}^ni^2S_G(\lfloor\frac{n}{i}\rfloor)=\frac{n*(n+1)}{2}\\ \Rightarrow S_G(n)=\frac{n*(n+1)}{2}-\sum_{i=2}^ni^2S_G(\lfloor\frac{n}{i}\rfloor)\\ \]
代码：

#include<bits/stdc++.h>
#define ll long long
#define maxx 3000005

using namespace std;
inline ll Get() {ll x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}while('0'<=ch&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}return x*f;}

const ll mod=1e9+7;

ll ksm(ll t,ll x) {
    ll ans=1;
    for(;x;x>>=1,t=t*t%mod)
        if(x&1) ans=ans*t%mod;
    return ans;
}

ll n;
int p[maxx];
ll f[maxx];
bool vis[maxx];
const ll inv2=ksm(2,mod-2),inv6=ksm(6,mod-2);
ll cal(ll n) {return n*(n+1)%mod*inv2%mod;}
ll cal2(ll n) {return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;}

void pre(int n) {
    for(int i=2;i<=n;i++) {
        if(!vis[i]) p[++p[0]]=i,f[i]=1-i+mod;
        for(int j=1;j<=p[0]&&1ll*i*p[j]<=n;j++) {
            vis[i*p[j]]=1;
            if(i%p[j]==0) {
                f[i*p[j]]=f[i];
                break;
            }
            f[i*p[j]]=(1-p[j])*f[i]%mod;
        }
    }
    f[1]=1;
    for(int i=1;i<=n;i++) {
        f[i]=((f[i]*i+f[i-1])%mod+mod)%mod;
    }
}

map<ll,ll>st;
ll Sum(ll n) {
    if(n<=3000000) return f[n];
    if(st.find(n)!=st.end()) return st[n];
    ll ans=cal(n%mod);
    ll last=1;
    for(ll i=2;i<=n;i=last+1) {
        ll now=n/(n/i);
        ans=(ans-(cal2(now%mod)-cal2(last%mod)+mod)*Sum(n/i)%mod+mod)%mod;
        last=now;
    }
    return st[n]=ans;
}

ll solve(ll n) {
    ll ans=0;
    ll last=0;
    for(ll i=1;i<=n;i=last+1) {
        ll now=n/(n/i);
        (ans+=(Sum(now)-Sum(last)+mod)*cal(n/i%mod)%mod*cal(n/i%mod))%=mod;
        last=now;
    }
    return ans;
}

int main() {
    pre(3000000);
    n=Get();
    cout<<solve(n);
    return 0;
}

转载于:https://www.cnblogs.com/hchhch233/p/10645777.html