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

(Description)

给你(n,m)，求(sum_{i=1}^nsum_{j=1}^mlcm(i,j))

答案对(100000009)取模。

多组数据。

(Input)

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

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

(Output)

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

(Sample Input)

1
4 5

(Sample Output)

122

(HINT)

对于(100\%)的数据：(t≤10000,n,m≤107)

(100000009)不是一个质数。

---

这道题一看就是一个数学题，一个莫比乌斯反演的入门题，如果没有了解过莫比乌斯反演的点这里

---

这道题给的很直接，于是我们就直接开始通过莫比乌斯反演导式子

$$sum_{i=1}^nsum_{j=1}mlcm(i,j)$$

由小学数学可得，(frac{ij}{gcd(i,j)}=lcm(i,j))

$$=sum_{i=1}^nsum_{j=1}mfrac{ij}{gcd(i,j)}$$

我们发现(frac{ij}{gcd(i,j)})这个式子看着难受，于是设(d=gcd(i,j))，枚举(d)，但是因为枚举的(d)不一定是(gcd(i,j))，于是最后判断一下

$$=sum_{d=1}^nsum_{i=1}nsum_{j=1}^mfrac{ij}{d}[d=gcd(i,j)]$$

将(d)除出来得到，因为其中(i、j)除于(d)后变成了(frac{i}{d}、frac{j}{d})，所以要使(frac{ij}{d})变化后值不变，就得到(frac{frac{i}{d}frac{j}{d}*d^2}{d})，又因为(i,j)在循环中已经除于了 (d)，所以(frac{ij}{d})变成了(frac{ijd^2}{d}=dij)

$$=sum_{d=1}^nsum_{i=1}{left lfloor frac{n}{d} ight floor}sum_{j=1}^{left lfloor frac{m}{d} ight floor}dij[gcd(i,j)=1]$$

又有莫比乌斯函数的性质：(sum_{dmid n}mu(d)=[n=1])，可以考虑设(d)为(k)，(n)为(gcd(i,j))，可得

$$sum_{kmid gcd(i,j)}mu(k)=[gcd(i,j)=1]$$

带入原式中得

$$=sum_{d=1}^nsum_{i=1}{left lfloor frac{n}{d} ight floor}sum_{j=1}^{left lfloor frac{m}{d} ight floor}dijsum_{kmid gcd(i,j)}mu(k)$$

因为如果要满足(kmid gcd(i,j))，则只需满足(kmid i且kmid j)，就可以把(k)提到前面去，则原式可以变成

$$=sum_{k=1}^{nmu(k)sum_{d=1}}nsum_{i=1}^{leftlfloorfrac{n}{d} ight floor}[kmid i]sum_{j=1}^{leftlfloorfrac{m}{d} ight floor}[kmid j] dij$$

将式子除于(k)得到

$$=sum_{k=1}^{nmu(k)sum_{d=1}}nsum_{i=1}^{{leftlfloorfrac{n}{dk} ight floor}sum_{j=1}}{leftlfloorfrac{m}{dk} ight floor}dk^2ij$$

设(dk=T)，并提前枚举(T)，则有

$$=sum_{T=1}^nsum_{i=1}{leftlfloorfrac{n}{T} ight floor}isum_{j=1}^{leftlfloorfrac{m}{T} ight floor}jsum_{dmid T}frac{T^2}{d}mu(frac{T}{d})$$

这就是最终的式子，现在我们考虑对于(sum_{i=1}^{leftlfloorfrac{n}{T} ight floor}isum_{j=1}^{leftlfloorfrac{m}{T} ight floor})用等差数列求和(O(1))求出，(sum_{dmid T}frac{T^2}{d}mu(frac{T}{d}))用线性筛(O(n))求出

现在对于一个函数(g(x)=sum_{dmid T}frac{T^2}{d}mu(frac{T}{d}))

1. 若(x)为质数，则通过找规律可得(g(x)=x-x^2)
2. 若(x=i*prime)且(i\%prime=0)，则(g(x)=g(i)*prime)
3. 若(x=i*prime)且(i%prime!=0) 即 (gcd(i,prime)=1)，则(g(x)=g(i)*g(prime))
   (ps:2,3)条性质都是积性函数的性质

最外层的循环就用整除分块 (O(sqrt n))处理

#include<bits/stdc++.h>
using namespace std;
const int N=10000010,mod=100000009;
typedef long long ll;
bool vis[N];
ll sum[N];
ll prime[N];
ll g[N];
int cnt=0;
void init()
{
    sum[1]=g[1]=1;
    for(int i=2;i<N;i++)
    {
        sum[i]=1LL*(i+1)*i/2%mod;
        if(!vis[i])
        {
            prime[++cnt]=i;
            g[i]=(i-1LL*i*i%mod+mod)%mod;
        }
        for(int j=1;j<=cnt&&i*prime[j]<N;j++)
        {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)g[i*prime[j]]=g[i]*prime[j]%mod;
            else g[i*prime[j]]=g[i]*g[prime[j]]%mod;
        }
    }
    for(int i=1;i<N;i++)g[i]=(g[i]+g[i-1])%mod;
}
ll ans=0;
void solve(int a,int b)
{
    for(int l=1,r=0;l<=a;l=r+1)
    {
        r=min(a/(a/l),b/(b/l));
        ans=(ans+sum[a/l]*sum[b/l]%mod*(g[r]-g[l-1]+mod)%mod)%mod;
    }
}
int main()
{
    int t;
    scanf("%d",&t);
    init();
    while(t--)
    {
        ans=0ll;
        int n,m;
        scanf("%d %d",&n,&m);
        if(n>m)swap(n,m);
        solve(n,m);
        printf("%lld
",ans);   
    }
    return 0;
}
/*
1
4 5
*/