数论练习题

数论

---

Prob 1:

给你 $N,M$，求 $\sum_{a=1}^N\sum_{b=1}^M gcd(a,b)$

思路一：$\sum_{d\mid n}\varphi(d)=n$

推导如下：

$$\sum_{a=1}^N\sum_{b=1}^Mgcd(a,b)$$
$$=\sum_{a=1}^N\sum_{b=1}^M\sum_{d\mid gcd(a,b)}\varphi(d)$$
$$=\sum_d\varphi(d)\cdot\sum_{a=1}^N\sum_{b=1}^M[d\mid gcd(a,b)]$$
$$=\sum_d\varphi(d)\cdot\sum_{a=1}^N\sum_{b=1}^M[d\mid a\ \wedge\ d\mid b]$$
$$=\sum_d\varphi(d)\cdot\sum_{a=1}^N[d\mid a]\cdot\sum_{b=1}^M[d\mid b]$$
$$=\sum_d\varphi(d)\cdot\lfloor\frac{N}{d}\rfloor\cdot\lfloor\frac{M}{d}\rfloor$$

化简到这里，我们不禁想起了另一道经典题目：
给定n，求 $\sum_{i=1}^N\lfloor\frac{N}{i}\rfloor$
解决方法：通过观察（打表）、猜想……我们可以发现，对于很多$i$而言，$\lfloor\frac{N}{i}\rfloor$的值是相等的。不妨设 $\forall i\in[l,r],\lfloor\frac{N}{i}\rfloor=d$，那么$\sum_{i=l}^r\lfloor\frac{N}{i}\rfloor=d*\sum_{i=l}^r$，然而，我们可以证明$[l,r]$这样的区间总共不超过$\sqrt{N}$个，问题可以在$O(\sqrt{N})$内解决。

我们再回到这道题来。
我们会发现，$\lfloor\frac{N}{d}\rfloor$和$\lfloor\frac{M}{d}\rfloor$，分别有$\sqrt{N},\sqrt{M}$个区间，那么$\lfloor\frac{N}{d}\rfloor\cdot\lfloor\frac{M}{d}\rfloor$取值不超过$\sqrt{N}+\sqrt{M}$个区间，我们可以预处理出来$\varphi(d)$的前缀和，然后枚举区间，问题可以在$O(\sqrt{N}+\sqrt{M})$解决。
粘上代码：

#include<bits/stdc++.h>
using namespace std;
const int N=100000,M=N+5;
#define LL long long 
int n,m;
bool cjx[M];
int prime[M],tot;
int phi[M];
LL sum[M];
void pre(){
    tot=0;
    phi[1]=1;
    for(int i=2;i<=N;i++){
        if(!cjx[i]) {
            prime[tot++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<tot;j++){
            if(i*prime[j]>N) break;
            cjx[i*prime[j]]=true;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            } else {
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            }
        }
    }
    for(int i=1;i<=N;i++) sum[i]=sum[i-1]+phi[i];
}
LL solve(int n,int m){
    LL re=0;
    for(int i=1,j;i<=min(n,m);i=j+1){
        j=min(n/(n/i),m/(m/i));
        re+=(sum[j]-sum[i-1])*(n/i)*(m/i);
    }
    return re;
}
int main(){
    pre();
    int T;  scanf("%d",&T);
    while(T--){
        scanf("%d%d",&n,&m);
        cout<<solve(n,m)<<endl;
    }
    return 0;
}

思路二：$\sum_{d\mid n}\mu(d)=[n=1]$

我们换一个思路，考虑$d$可能成为哪些$i,j$的最大公约数

$$\sum_{i=1}^N\sum_{j=1}^Mgcd(i,j)$$
$$=\sum_{i=1}^N\sum_{j=1}^M\sum_d[d=gcd(i,j)]d$$
$$=\sum_dd\cdot\sum_{i=1}^N\sum_{j=1}^M[d=gcd(i,j)]$$
$$=\sum_dd\cdot\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}[gcd(i,j)=1]$$
$$=\sum_dd\cdot\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}\sum_{d'\mid gcd(i,j)}\mu(d')$$
$$=\sum_dd\cdot\sum_{d'}\mu(d')\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}[d'\mid i]\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}[d'\mid j]$$
$$=\sum_dd\cdot\sum_{d'}\mu(d')\cdot\lfloor\frac{\lfloor\frac{N}{d}\rfloor}{d'}\rfloor\cdot\lfloor\frac{\lfloor\frac{M}{d}\rfloor}{d'}\rfloor$$
$$=\sum_dd\cdot\sum_{d'}\mu(d')\cdot\lfloor\frac{N}{d\cdot d'}\rfloor\cdot\lfloor\frac{M}{d\cdot d'}\rfloor$$
我们设$D=d\cdot d'$，那么$d=\frac{D}{d'}$
$$ans=\sum_D\lfloor\frac{N}{D}\rfloor\cdot\lfloor\frac{M}{D}\rfloor\cdot\sum_{d\mid D}\mu(d)\cdot \frac{D}{d}$$
我们发现后面的$H(x)=\sum_{d\mid x}\mu(d)\cdot \frac{x}{d}$是$f(x)=\mu(x)$和$g(x)=x$的狄利克雷卷积，所以$H(x)$是积性函数，我们可以使用筛法$O(N)$预处理出来，剩下过程和思路一相同了。
到了这里，不知道大家有没有发现，和思路一对比之下，我们发现，$H(x)=\varphi(x)$，没错！
这是巧合吗？不！
$n=\sum_{d\mid n}\varphi(d)$
我们设$f(x)=x,g(x)=\varphi(x)$
则$f(x)=\sum_{d\mid x}g(d)$，
我们进行莫比乌斯反演得到：
$g(x)=\sum_{d\mid x}\mu(d)\cdot f(\frac{x}{d})$
也就是
$\varphi(x)=\sum_{d\mid x}\mu(d)\cdot \frac{x}{d}$
我们结合思路一与思路二无意中发现了莫比乌斯反演，哈哈哈。。。

---

Prob 2：

给你 $N,M$，求$\sum_{a=1}^N\sum_{b=1}^Mlcm(a,b)$

主要公式：$\sum_{d\mid n}\mu(d)=[n=1]$

推导如下：

$$\sum_{a=1}^N\sum_{b=1}^Mlcm(a,b)$$
$$=\sum_{a=1}^N\sum_{b=1}^M\frac{a\cdot b}{gcd(a,b)}$$
$$=\sum_{a=1}^N\ \ \sum_{b=1,d=gcd(a,b)}^M\frac{a\cdot b}{d}$$
$$=\sum_d\sum_{a=1}^N\sum_{b=1}^M\frac{a\cdot b}{d}\cdot [d=gcd(a,b)]$$
$$=\sum_dd\cdot\sum_{a=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{M}{d}\rfloor}a\cdot b\cdot [gcd(a,b)=1]$$
不妨设$f(n,m)=\sum_{a=1}^n\sum_{b=1}^ma\cdot b\cdot [gcd(a,b)=1]$
下面化简$f(n,m)=$
$$\sum_{a=1}^n\sum_{b=1}^ma\cdot b\cdot [gcd(a,b)=1]$$
$$=\sum_{a=1}^n\sum_{b=1}^ma\cdot b\cdot\sum_{d\mid gcd(a,b)}\mu(d)$$
$$=\sum_d\mu(d)\sum_{d\mid a}^n\sum_{d\mid b}^ma\cdot b$$
$$=\sum_d\mu(d)\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}a\cdot b\cdot d^2$$
$$=\sum_d\mu(d)\cdot d^2\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}a\cdot b$$
$$=\sum_d\mu(d)\cdot d^2\cdot \frac{\lfloor\frac{n}{d}\rfloor\cdot(\lfloor\frac{n}{d}\rfloor+1)}{2}\cdot\frac{\lfloor\frac{m}{d}\rfloor\cdot(\lfloor\frac{m}{d}\rfloor+1)}{2}$$
下面我们把$f(\lfloor\frac{N}{d}\rfloor,\lfloor\frac{M}{d}\rfloor)$代入原式：$ans=$
$$\frac{1}{4}\sum_dd\cdot\sum_{d'}\mu(d')\cdot d'^2\cdot \lfloor\frac{\lfloor\frac{N}{d}\rfloor}{d'}\rfloor\cdot (\lfloor\frac{\lfloor\frac{N}{d}\rfloor}{d'}\rfloor+1)\cdot \lfloor\frac{\lfloor\frac{M}{d}\rfloor}{d'}\rfloor\cdot (\lfloor\frac{\lfloor\frac{M}{d}\rfloor}{d'}\rfloor+1)$$
$$=\frac{1}{4}\sum_dd\cdot\sum_{d'}\mu(d')\cdot d'^2\cdot \lfloor\frac{N}{d\cdot d'}\rfloor\cdot (\lfloor\frac{N}{d\cdot d'}\rfloor+1)\cdot \lfloor\frac{ M}{d\cdot d'}\rfloor\cdot (\lfloor\frac{ M}{d\cdot d'}\rfloor+1)$$
令$D=d\cdot d'$
$$=\frac{1}{4}\sum_D \lfloor\frac{N}{D}\rfloor\cdot (\lfloor\frac{N}{D}\rfloor+1)\cdot \lfloor\frac{ M}{D}\rfloor\cdot (\lfloor\frac{ M}{D}\rfloor+1)\cdot D\sum_{d\mid D}\mu(d)\cdot d$$
由于积性函数的莫比乌斯变换函数还是积性函数，所以$D\sum_{d\mid D}\mu(d)\cdot d$是积性函数，可以用筛法$O(N)$，预处理出来。我们再观察前面$\frac{1}{4}\sum_D \lfloor\frac{N}{D}\rfloor\cdot (\lfloor\frac{N}{D}\rfloor+1)\cdot \lfloor\frac{ M}{D}\rfloor\cdot (\lfloor\frac{ M}{D}\rfloor+1)$这一堆东西，很明显还是分块来搞。
这样一来，问题最终$O(\sqrt{N}+\sqrt{M})$完美解决。
粘上代码：

#include<bits/stdc++.h>
using namespace std;
const int N=100000;
const int M=100005;
typedef long long LL;
bool cjx[M];
int prime[M],tot;
int mu[M];
LL g[M],sum[M];
void pre(){
    tot=0;
    mu[1]=1;    g[1]=1;
    for(int i=2;i<=N;i++){
        if(!cjx[i]){
            prime[tot++]=i;
            mu[i]=-1;
            g[i]=-i*(i-1);
        }
        for(int j=0;j<tot;j++){
            if(i*prime[j]>N) break;
            cjx[i*prime[j]]=true;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                g[i*prime[j]]=g[i]*prime[j];
                break;
            }   else {
                mu[i*prime[j]]=-mu[i];
                g[i*prime[j]]=g[i]*g[prime[j]];
            }
        }
    }
    for(int i=1;i<=N;i++) sum[i]=sum[i-1]+g[i];
}
LL solve(int n,int m){
    LL re=0;
    for(int i=1,j;i<=min(n,m);i=j+1){
        j=min(n/(n/i),m/(m/i));
        re+=(LL)(sum[j]-sum[i-1])*(n/i)*(n/i+1)*(m/i)*(m/i+1);
    }
    return re>>2;
}
int main(){
    pre();
    int T;  scanf("%d",&T);
    while(T--){
        int n,m;
        scanf("%d%d",&n,&m);
        cout<<solve(n,m)<<endl;
    }
    return 0;
}

---

Prob 3

给你 $n$，求 $\sum_{i=1}^ngcd(i,n)$

解析：
首先转化思路，我们思考$d$可以成为哪些$i$与$n$的最大公约数
$\sum_{i=1}^ngcd(i,n)=\sum_{d\mid n}d\cdot \varphi(\frac{n}{d})$
设$f(x)=x,g(x)=\varphi(x),H(x)=\sum_{d\mid n}d\cdot \varphi(\frac{n}{d})$
则$H(x)$是$f(x)$和$g(x)$的狄利克雷卷积，所以$H(x)$是积性函数。
我们将$n$分解质因数
$H(n)=H(p_1^{\alpha1}\cdot p_2^{\alpha2}\cdot\ldots\cdot p_n^{\alpha n})=H(p_1^{\alpha1})\cdot H(p_2^{\alpha2})\cdot\ldots\cdot H(p_n^{\alpha n})$
所以我们只需要思考怎么计算$H(p^k)$即可。
$H(p^k)=\sum_{d\mid p^k}d\cdot \varphi(\frac{p^k}{d})$
$=\sum_{i=0}^kp^i\cdot \varphi(p^{k-i})$
而我们知道$\varphi(p^k)=p^k-p^{k-1}$
那么我们通过数列化简、整理最终得到
$H(p^k)=(k+1)\cdot p^k-k\cdot p^{k-1}$
问题解决了，我的代码比较丑，大家意识流即可……

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long LL;
LL pow(LL a,LL b){
    LL re=1;
    for(;b;b>>=1,a*=a)
        if(b&1) re*=a;
    return re;
}
int main(){
    LL n;   cin>>n;
    LL ans=1;
    for(int i=2;i*i<=n;i++){
        LL a=0;
        while(n%i==0){
            a++;    n/=i;
        }
        if(a==0) continue;
        LL b=pow(i,a-1);
        ans*=(a+1)*b*i-a*b;
    }
    ans*=2*n-1;
    cout<<ans<<endl;
    return 0;
}