SPOJ - GCDMAT 莫比乌斯反演

Problem

---

题目描述

给定一定大小的矩阵，矩阵中的每个元素的贡献为其行列坐标的gcd。

举个例子，一个3×2的矩阵有这些元素：

$$\left[\begin{aligned}gcd(1,1)\quad gcd(1,2)\\gcd(2,1)\quad gcd(2,2)\\gcd(3,1)\quad gcd(3,2)\\ \end{aligned}\right]$$

你需要回答这个矩阵中以(i1,j1)为左上角，(i2,j2)为右下角的矩形范围内的元素之和。

输入输出格式

输入格式
第一行一个整数T，表示数据组数。(T<=500)

接下来一行是两个由空格隔开的整数，n和m。(1<=n,m<=50000)

接下来T行每行包含一组询问i1,j1,i2,j2。保证i1<=i2,j1<=j2。

输出格式
对于每个询问输出一行一个整数表示答案对10^9+7取模后的结果。

Solution

---

我们要求$\sum_{i=i_1}^{i_2}\sum_{j=j_1}^{j_2}gcd(i,j)$，只需求出$\sum_{i=1}^x\sum_{j=1}^ygcd(i,j)$。

首先我们可以枚举可能的gcd，即$\sum_{d=1}^nd\sum_{i=1}^x\sum_{j=1}^y(gcd(i,j)==d)$

即相当于$\sum_{d=1}^nd\sum_{i=1}^{\frac x d}\sum_{j=1}^{\frac y d}(gcd(i,j)==1)$

---

证明一
不妨设$f(x)=\sum_{i=1}^a\sum_{j=1}^b(gcd(i,j)==x)$,$g(x)=\sum_{i=1}^a\sum_{j=1}^b(gcd(i,j)==x*k)=\frac a x*\frac b x$
显然$g(x)=\sum_{d|x} f(d)$，由莫比乌斯反演，我们可以知道$f(x)=\sum_{d|x}\mu(d)g(\frac x d)$

---

那么原式化为$\sum_{d=1}^nd\sum_{e|d}\mu(e)\frac x {de}\frac y {de}$

然后你会发现式子十分的复杂。那么我们就重设参数$r=de$
则有$\sum_{r=1}^n\sum_{d|r}d*\mu(\frac r d) \frac x r \frac y r$

---

证明二
我们知道对于$x=p_1p_2p_3…p_k$，有$\mu(x)=(-1)^k$，否则$\mu(x)=0$
其次，对于$\phi(x)=x(1-\frac 1 {p_1})(1-\frac 1 {p_2})…(1-\frac 1 {p_k})$，我们把除x之外的所有数都乘起来，得

$$\phi(x)=x(1-\sum_{i=1}^k\frac 1 {p_i}+\sum_{i,j=1,i\not=j}^k\frac 1 {p_ip_j}…)$$
而我们可以发现其式子中的 $p_i,p_ip_j,p_ip_jp_r…$就是x的所有约数，而其相应的符号则与上面的莫比乌斯函数之定义相符，也就是说 $\phi(x)=\sum_{d|x} \frac x d \mu(d)=\sum_{d|x} d\mu(\frac x d)$

---

由此，原式变为$\sum_{r=1}^n\phi(r) \frac x r \frac y r$

最后是同样的套路，我们对$\frac x r \frac y r$进行分块处理，对欧拉函数求个前缀和，就可以做到$\sqrt n$的复杂度。

Code

---

#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=50000,mod=1000000007;
int z,n,m,tot,i1,j1,i2,j2,ans;
int phi[maxn+10],sum[maxn+10],pri[maxn+10],vis[maxn+10];
inline int min(int x,int y){return x<y?x:y;}
inline void pls(int &x){if(x>=mod) x-=mod;}
void init()
{
    sum[1]=phi[1]=1;
    for(int i=2;i<=maxn;i++)
    {
        if(!vis[i]) pri[++tot]=i,phi[i]=i-1;
        for(int j=1;j<=tot&&i*pri[j]<=maxn;j++)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) {phi[i*pri[j]]=phi[i]*pri[j];break;}
            else phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
        sum[i]=phi[i]+sum[i-1];
    }
}
int solve(int x,int y)
{
    int res=0,lim=min(x,y);
    for(int i=1,j;i<=lim;i=j+1)
    {
        j=min(x/(x/i),y/(y/i));
        pls(res+=(ll)(sum[j]-sum[i-1])*(x/i)%mod*(y/i)%mod);
    }
    return res;
}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("in.txt","r",stdin);
    #endif
    scanf("%d%d%d",&z,&n,&m);
    init();
    while(z--)
    {
        scanf("%d%d%d%d",&i1,&j1,&i2,&j2);
        ans=((ll)solve(i2,j2)-solve(i1-1,j2)-solve(i2,j1-1)+solve(i1-1,j1-1)+2*mod)%mod;
        printf("%d\n",ans);
    }
    return 0;
}