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
    
     #define ll long long#define maxx 3000005using 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
     
      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<