【BZOJ2693】jzptab （莫比乌斯反演）

Description

给你$n$,$m$，求 $sum^n_{i=1} sum^m_{j=1} lcm(x,y)$

答案对$100000009$取模。

多组数据。

Input

第一行有一个正整数$t$表示数据组数

接下来$t$行每行有两个正整数$n$,$m$

Output

$t$行，第$i$行为第$i$组询问的答案。

HINT

对于$100\%$的数据：$tleq 10000,n,mleq {10}^7$

$100000009$不是一个质数。

 题解：

第一次打莫演，手推公式。（我为什么会做这道题）

egin{aligned}
ans&=sum^n_{i=1}sum^m_{j=1}frac{ij}{gcd(i,j)}\
&=sum^n_dsum^n_{i=1}sum^m_{j=1}frac{ij}{d}[gcd(i,j)=d]\
&=sum^n_dsum^{lfloorfrac{n}{d} floor}_{i=1}sum^{lfloorfrac{m}{d} floor}_{j=1}frac{d^2ij}{d}[gcd(i,j)=1]\
&=sum^n_dsum^{lfloorfrac{n}{d} floor}_{i=1}sum^{lfloorfrac{m}{d} floor}_{j=1}dijsum_{kmid gcd(i,j)}mu(k)&(sum_{dmid n}mu(d)=[n=1])\
&=sum^n_kmu(k)sum^n_dsum^{lfloorfrac{n}{d} floor}_{i=1}[kmid i]sum^{lfloorfrac{m}{d} floor}_{j=1}[kmid j]dij\
&=sum^n_kmu(k)sum^n_dsum^{lfloorfrac{n}{kd} floor}_{i=1}sum^{lfloorfrac{m}{kd} floor}_{j=1}dk^2ij\
&=sum^n_kmu(k)sum^n_dsum^{lfloorfrac{n}{kd} floor}_{i=1}isum^{lfloorfrac{m}{kd} floor}_{j=1}jcdot dk^2\
&设 T=kd\
ans&=sum^n_{T=1}sum^{lfloorfrac{n}{T} floor}_{i=1}isum^{lfloorfrac{m}{T} floor}_{j=1}jsum_{dmid T}mu(frac{T}{d})frac{T^2}{d}
end{aligned}

第一个部分 $ sum^{lfloorfrac{n}{T} floor}_{i=1}isum^{lfloorfrac{m}{T} floor}_{j=1}j $ 用等差数列求和 $O(1)$ 求出

第二个部分 $ sum_{dmid T}mu(frac{T}{d})frac{T^2}{d} $ 线性筛 $O(n)$ 预处理处

设$ g(x)=sum_{dmid T}mu(frac{T}{d})frac{T^2}{d} $，考虑怎么求$g(x)$

如果说x为质数，那么根据公式$g(x)=x-x^2$

如果$x$不为质数，我们设$x=i imes p$，其中$p$为质数，那么有两种情况

$p mid i$，由于$i$和$p$互质而$g(x)$为积性函数，$g(x)=g(i imes p)=g(i) imes g(p)$

$pmid i$，这个时候就有点不是很好搞了……

我们可以把i表示为$t imes p^k$（$t$与$p$互质）

那么我们就尝试一下从乘了一个$p$会有什么影响这个方面来考虑一下

考虑$g(p^k)$的值，显然根据$mu$的定义，只有$mu(1)$和$mu(p)$能够提供贡献（其他的$p$的指数都$>1$，所以都是$0$）

那么我们就可以得到$g(p^k)=f(1)p^k+f(p)p^{k-1}$

然后写出$g(p^{k+1})$的表达式，会发现是$f(1)p^{k+1}+f(p)p^k$

也就是说$g(p^{k+1})=g(p^k)p$

那么就可以得到$g(x)=g(i^p)=g(t imes p^k imes p)=g(t) imes g(p^k) imes p=g(x) imes p$

然后就可以顺利筛出来啦

最外层循环用数论分块，总时间 $O(sqrt{n})$

CODE:

#include<iostream>
#include<cstdio>
using namespace std;

#define mod 100000009LL
#define N 10000005
int t,n,m,cnt,ans;
long long pri[N],g[N],sum[N];
bool vis[N];

void init(){
    sum[1]=g[1]=1;
    for(int i=2;i<N;i++){
        sum[i]=1LL*i*(i+1)/2%mod;
        if(!vis[i]){
            g[i]=(i-1LL*i*i%mod+mod)%mod;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*pri[j]<N;j++){
            vis[i*pri[j]]=true;
            if(i%pri[j])
                g[i*pri[j]]=g[i]*g[pri[j]]%mod;
            else
                g[i*pri[j]]=g[i]*pri[j]%mod;
        }
    }
    for(int i=1;i<N;i++)(g[i]+=g[i-1])%=mod;
}

int main(){
    scanf("%d",&t);
    init();
    while(t--){
        scanf("%d%d",&n,&m);
        if(n>m)swap(n,m);
        ans=0;
        for(int i=1,pos=0;i<=n;i=pos+1){
            pos=min(n/(n/i),m/(m/i));
            ans+=sum[n/i]*sum[m/i]%mod*(g[pos]-g[i-1]+mod)%mod;
            ans%=mod;
        }
        printf("%d
",ans);
    }
}