LCM(i,j)求和 （莫比乌斯函数）（NKOJ3958/BZOJ2693）

题意：

求$\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)$
n,m=1e7，多组输入

---

$O(N)$版本：

遇到lcm首先是换成gcd，$Ans=\sum_{i=1}^n\sum_{j=1}^m\dfrac{i*j}{gcd(i,j)}$，再换成因子角度，设：$gcd(i,j)=d,min(n,m)=N$，有：

$$\sum_{d=1}^N\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]\dfrac{i*j}{d}$$ 而 $gcd==d$按照套路可以进行莫比乌斯的倍数反演： $$g(d)=\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]i*j（去掉\sum_{d=1}^N和/d才是积性函数，才容易反演）\\构造：f(d)=\sum_{d|D}^Ng(D)=\sum_{i=1}^n\sum_{j=1}^m[d|gcd(i,j)]i*j\\显然：f(d)=d^2*\dfrac{(1+[\frac{n}{d}])([\frac{n}{d}])}{2}*\dfrac{(1+[\frac{m}{d}])([\frac{m}{d}])}{2}\\反演：g(d)=\sum_{d|D}^Nu(\frac{D}{d})f(D)\\那么：g(1)=\sum_{i=1}^Nu(i)f(i)\\代入：Ans=\sum_{d=1}^N\frac{1}{d}\sum_{d|D}^Nu(\frac{D}{d})f(D)\\令i=D/d得：Ans=\sum_{d=1}^N\frac{1}{d}\sum_{i=1}^{[\frac{N}{d}]}u(i)f(id)\\最后带入f得：Ans=\sum_{d=1}^Nd\sum_{i=1}^{[\frac{N}{d}]}u(i)*i^2*\dfrac{(1+[\frac{n}{id}])([\frac{n}{id}])}{2}*\dfrac{(1+[\frac{m}{id}])([\frac{m}{id}])}{2}$$

观察发现$[\frac{N}{d}]$和$[\frac{n}{id}]$可以整除分块，那么时间复杂度为$O(\sqrt n*\sqrt n)=O(n)$，一般的题目是可以过了，但是案例数多的就不行了

---

$O(\sqrt N)$版本：

$$Ans=\sum_{d=1}^Nd\sum_{i=1}^{[\frac{N}{d}]}u(i)*i^2*\dfrac{(1+[\frac{n}{id}])([\frac{n}{id}])}{2}*\dfrac{(1+[\frac{m}{id}])([\frac{m}{id}])}{2}\\=\sum_{d=1}^Nd\sum_{i|d}u(i)*i^2*\dfrac{(1+[\frac{n}{id}])([\frac{n}{id}])}{2}*\dfrac{(1+[\frac{m}{id}])([\frac{m}{id}])}{2} （这步在写杜教筛（上）中有证明）$$

改变枚举变量，令$D=i*d得：$

$$Ans=\sum_{D=1}^ND*\dfrac{(1+[\frac{n}{D}])([\frac{n}{D}])}{2}*\dfrac{(1+[\frac{m}{D}])([\frac{m}{D}])}{2}\sum_{i|D}u(i)*i$$

想办法快速求出$\sum_{i|D}u(i)*i：$

$$\because u(i)*i为积性\\\therefore\sum_{i|D}u(i)*i为积性\\设：F(D)=\sum_{i|D}u(i)*i\\显然：F(p^q)=1-p\\那么：F(p_1^{q_1}*p_2^{q_2}*…p_k^{q_k})=(1-p_1)*...(1-p_k)\\\therefore D*F(D)可以直接放进线性筛预，再处理出前缀和\\\therefore只需要对[\frac{n}{D}]分块即可，时间复杂度O(\sqrt N)$$

代码：

#include<bits/stdc++.h>
#define LL long long
using namespace std;
int read(){int t;
    scanf("%d",&t);return t;
}

const int N=1e7+5;
const int mod=1e8+9;//注意模数
LL pri[N>>1],now;
bool vis[N];
LL sum[N];
init(){
    now=0;vis[1]=1;
    for(LL i=2;i<N;i++){//预处理D*F(D)
        if(!vis[i])pri[++now]=i,sum[i]=(i*(1-i)%mod+mod)%mod;
        for(LL j=1,mul=2*i;j<=now&&mul<N;j++,mul=pri[j]*i){
            vis[mul]=1;
            if(i%pri[j]==0){sum[mul]=sum[i]*pri[j]%mod;break;}
            sum[mul]=sum[i]*sum[pri[j]]%mod;
        }
    }
    sum[1]=1;
    for(int i=2;i<N;i++)sum[i]=(sum[i]+sum[i-1])%mod;//变成前缀数组
}

LL cal(LL x,LL y){
    return ((1+x)*x/2%mod)*((1+y)*y/2%mod)%mod;
}

int main(){
    init();
    int t=read();
    while(t--){
        LL n=read(),m=read();
        LL NN=min(n,m);
        LL ans=0;
        for(LL l=1,r;l<=NN;l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans=(ans+(sum[r]-sum[l-1]+mod)%mod*cal(n/l,m/l)%mod)%mod;
        }
        printf("%lld\n",ans);
    }
}