【NKOJ3958/BZOJ2693】【莫比乌斯反演】jzptab

题目要求：

$$\sum_{i=1}^n \sum_{j=1}^m lcm(i,j)$$

看到lcm，首先想到转换为做过的gcd:

$$\sum_{i=1}^n \sum_{j=1}^m \frac{i*j}{gcd(i,j)}$$

相同的套路，设：

$$d=gcd(i,j)$$
$$N=min(n,m)$$
$$f(n,m,d)=\sum_{i=1}^n \sum_{j=1}^m i*j [gcd(i,j)=d]$$
$$原式=\sum_{d=1}^{N}\frac{f(n,m,d)}{d}$$
$$将f(d)中i,j的d因子提出来：\sum_{d=1}^{N}d*f( \lfloor\frac n d \rfloor,\lfloor\frac m d\rfloor,1)$$
$$开始表演，设：F(n,m,d)= \sum_{i=1}^n \sum_{j=1}^m i*j[d|gcd(i,j)$$
$$显然F(n,m,d)=\sum_{d|e}^Nf(n,m,e)$$
$$反演一下：f(n,m,d)=\sum_{d|e}^N\mu(\frac e d)*F(n,m,e)$$
$$设sum(n,m)=\sum_{i=1}^n\sum_{j=1}^m i*j=\frac{n(n+1)} 2 *\frac{m(m+1)} 2$$
$$显然F(n,m,d)=d^2*sum(\lfloor \frac n d \rfloor,\lfloor \frac m d \rfloor)$$
$$f(n,m,1)=\sum_{i=1}^{N}\mu(i)*i^2*sum(\lfloor \frac n i \rfloor,\lfloor \frac m i\rfloor)$$
$$Ans=\sum_{d=1}^{N}d\sum_{i=1}^{ \lfloor \frac N d\rfloor}\mu(i)*i^2*sum(\lfloor \frac n {di}\rfloor,\lfloor \frac m {di}\rfloor)$$

到目前为止，我们可以利用分块优化，在第一个和式利用$\lfloor \frac N d\rfloor$进行分块，在第二个和式利用$\lfloor \frac n {di} \rfloor$进行分块，在$O(\sqrt n *\sqrt n)=O(n)$的复杂度下求出答案

然而题目还有T组数据的要求OTZ，卡${O(\sqrt n)}$过

于是回忆起董神的教诲，尝试改变枚举顺序，把比较烦的di直接转成枚举量

$$Ans=\sum_{d=1}^{min(n,m)}d*sum(\lfloor \frac n d \rfloor,\lfloor \frac m d \rfloor)\sum_{i|d}^d\mu(i)*i$$

前面求sum()肯定是$\log n$复杂度，于是我们需要在后面那一坨上面进行优化

设$g(d)=\sum_{i|d}^d\mu(i)*i$

显然$\mu(i)*i$为积性函数，根据莫比乌斯反演性质，$g(d)$也为积性函数

将$d$拆为$p_1^{k_1}*p_2^{k_2}*p_3^{k_3}*....*p_i^{k_i}$ ($p_i$为质因子)

对于每个$p_i^{k_i},g(p_i^{k_i})=1-pi$（显然）

$g(d)=\Pi_{i=1}^{d的质因子数}g(p_i^{k_i})=\Pi_{i=1}^k (1-p_i)$

很像$phi$，放在线性筛中处理

总复杂度$O(T \sqrt n)$

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=10000005,Inf=1e9;
const ll LInf=1e15,Mod=100000009;
inline ll _R(){
    ll d=0;bool b=1;char t=getchar();
    while(t<'0'||t>'9'){if(t=='-')b=0;t=getchar();}
    for(;t>='0'&&t<='9';t=getchar())d=(d<<3)+(d<<1)+t-'0';
    return b?d:-d;
}

ll T,n,m,t,Ans;
ll Cnt,Prime[N],Vis[N],Sum[N];
void Euler(){
     Sum[1]=1;
     for(ll i=2;i<=N-5;i++){
        if(!Vis[i]){
            Sum[i]=(i-i*i%Mod+Mod)%Mod;
            Prime[++Cnt]=i;
        }
        for(ll j=1;Prime[j]*i<=N-5&&j<=Cnt;j++){
            ll k=Prime[j]*i;
            Vis[k]=1;
            if(i%Prime[j])Sum[k]=Sum[Prime[j]]*Sum[i]%Mod;
            else {
                Sum[k]=Sum[i]*Prime[j]%Mod;
                break;
            }
        }
     }
     for(int i=1;i<=N-5;i++)Sum[i]=(Sum[i-1]+Sum[i])%Mod;
}
ll Cal(ll x,ll y){
    return (x*(x+1)/2%Mod)*(y*(y+1)/2%Mod)%Mod;
}
int main(){
    T=_R();
    Euler();
    while(T--){
        n=_R(),m=_R();
        Ans=0;
        for(ll i=1;i<=min(n,m);i++){
            t=min(n/(n/i),m/(m/i));
            Ans=(Ans+((Sum[t]-Sum[i-1]+Mod)%Mod)*Cal(n/i,m/i)%Mod)%Mod;
            i=t;
        }
        printf("%lld\n",Ans);
    }
}