莫比乌斯反演模板

莫比乌斯函数定义：

设 $n = p_{1}^{k_{1}} \cdot p_{2}^{k_{2}} \cdot \dots \cdot p_{m}^{k_{m}}$ ，其中 p 为素数，则定义如下：

$μ (n) = {\begin{matrix} <mstyle displaystyle="false" scriptlevel="0"> 1 </mstyle> & <mstyle displaystyle="false" scriptlevel="0"> n = 1 </mstyle> \\ <mstyle displaystyle="false" scriptlevel="0"> (- 1)^{m} </mstyle> & <mstyle displaystyle="false" scriptlevel="0"> \prod_{i = 1}^{m} k_{i} = 1 </mstyle> \\ <mstyle displaystyle="false" scriptlevel="0"> 0 </mstyle> & <mstyle displaystyle="false" scriptlevel="0"> <mtext> otherwise </mtext> (k_{i} > 1) </mstyle> \end{matrix}$

可知有平方因子的数的莫比乌斯函数值为 0。

由于Mobius函数是积性函数，即 $μ (a \cdot b) = μ (a) \cdot μ (b), (a, b) = 1$ ,所以
可以在 $O (n)$ 的时间内筛出 $[1, n]$ 的函数值，一个数 $x = p_{1}^{k_{1}} \cdot p_{2}^{k_{2}} \cdot \dots \cdot p_{m}^{k_{m}}$ ，
若 $k_{i} = 1, k_{j} > 1, j > i$ 则x只会被 $\frac{x}{p_{j}} m o d p_{j} = 0$ 删除

int mu[MAXN],prime[MAXN],cnt;
bool noprime[MAXN];
void mublus(){
    mu[1]=1;
    for(int i=2;i<MAXN;++i){
        if(!noprime[i]){
            prime[++cnt]=i;
            mu[i]=-1;//质数为-1
        }
        for(int j=1;j<=cnt&&prime[j]*i<MAXN;++j){
            noprime[prime[j]*i]=1;
            if(i%prime[j]==0){
                mu[prime[j]*i]=0;// prime[j]*i有因子prime[j]^2 ,故为0；
                break;
            }
            else mu[prime[j]*i]=-mu[i];// i 比prime[j]*i多了一个质因子。
        }
    }
}

$\sum_{d ∣ n} μ (d) = [n = 1]$

证明：当 n = 1 时有 $μ (1) = 1$ ，显然成立。否则设 $n = p_{1}^{k_{1}} \cdot p_{2}^{k_{2}} \cdot \dots \cdot p_{m}^{k_{m}}$ ， $d = p_{1}^{x_{1}} \cdot p_{2}^{x_{2}} \cdot \dots \cdot p_{m}^{x_{m}}$

根据 $μ (n)$ 的定义，只需考虑 $x_{i} = 0$ 或 $x_{i} = 1$ 的情况。我们设 d 中存在 r 个 $x_{i}$ 为 1，那么有：

$\sum_{d ∣ n} μ (d) = \sum_{r = 0}^{m} (\frac{m}{r}) (- 1)^{r} (n \neq 1)$

根据二项式定理：

$(x + y)^{n} = \sum_{k = 0}^{n} (\frac{n}{k}) x^{n - k} y^{k} (x + y)$

我们令 x = 1, y = -1，即得证：

$\sum_{r = 0}^{m} (\frac{m}{r}) (- 1)^{r} = (1 - 1)^{m} = 0 (n \neq 1)$

莫比乌斯反演

假设对于数论函数 $f (n)$ 和 $F (n)$ ，满足以下关系式： $F (n) = \sum_{d ∣ n} f (d)$

则将其默比乌斯反演公式定义为：
$f (n) = <munder> \sum d ∣ n </munder> μ (d) F (\frac{n}{d})$
或者满足： $F (n) = \sum_{n ∣ d} f (d)$

有
$f (n) = <munder> \sum n ∣ d </munder> μ (\frac{d}{n}) F (d)$

1.求区间[1,a], [1,b]中gcd=k的个数

$<munderover> \sum i = 1 a </munderover> <munderover> \sum j = 1 b </munderover> [g c d (i, j) = k]$
即求
$<munderover> \sum i = 1 ⌊ \frac{a}{k} ⌋ </munderover> <munderover> \sum j = 1 ⌊ \frac{b}{k} ⌋ </munderover> [g c d (i, j) = 1]$

令 $f (k)$ 为满足gcd=k的个数， $F (k)$ 为满足 $g c d = x \cdot k ， x > 0$ 的个数.

则有 $F (n) = \sum_{n ∣ d} f (d)$ （d是n的倍数， $\sum f (d)$ 显然就是F的定义）

所以有反演公式： $f (n) = \sum_{n ∣ d} μ (\frac{d}{n}) F (d)$

显然， $F (n) = ⌊ \frac{a}{n} ⌋ \cdot ⌊ \frac{b}{n} ⌋$

（a/n是第一个区间里n的倍数的个数，两者相乘则为 i，j 都是n的倍数的个数，所以gcd(i,j)是n的倍数，显然就是F的定义）

所以 $f (n) = \sum_{n ∣ d} μ (\frac{d}{n}) F (d) = \sum_{n ∣ d} μ (\frac{d}{n}) \cdot ⌊ \frac{a}{d} ⌋ \cdot ⌊ \frac{b}{d} ⌋$

由于k=1

有 $f (1) = \sum_{1 ∣ d} μ (d) F (d) = \sum_{1 ∣ d} μ (d) \cdot ⌊ \frac{a}{d} ⌋ \cdot ⌊ \frac{b}{d} ⌋$

HDU-1695

求[a,b],[c,d]里gcd=k的个数，a=c=1，gcd(1,2)和gcd(2,1)算作一个。

$a n s = f (1) = \sum_{1 ∣ i}^{\min (b, d)} μ (i) F (i) = \sum_{1 ∣ i} μ (i) \cdot ⌊ \frac{b}{i} ⌋ \cdot ⌊ \frac{d}{i} ⌋$
为何sum的上界是min(b,d)，因为i大于min时，F(i)=0,且F(i)表示两个区间里各取一个数，满足gcd=i的倍数，若i大于其中一个区间的上界，则不能满足gcd=i的倍数。
最后，ans是有重复的，例如[1,5]和[1,10]两个区间，[1,3]和[3,1]算作两次，但是[1,6],[6,1]中[6,1]不可能倍计算，所以只算了一次。所以，在[1,min(b,d)]里的答案都计算了两次。

#include<bits/stdc++.h>
using namespace std;
#define Init(arr,val) memset(arr,val,sizeof(arr))
const int inf=0x3f3f3f3f,mod=1000000007,MAXN=1e5+8;
typedef long long ll;
int mu[MAXN],prime[MAXN],cnt;
bool noprime[MAXN];
void mubius(){
    mu[1]=1;
    for(int i=2;i<MAXN;++i){
        if(!noprime[i]){
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&prime[j]*i<MAXN;++j){
            noprime[prime[j]*i]=1;
            if(i%prime[j]==0){
                mu[prime[j]*i]=0;
                break;
            }
            else mu[prime[j]*i]=-mu[i];
        }
    }
}
int main() {
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    std::cout.tie(0);
    mubius();
    int T,a,b,c,d,k;
    cin>>T;
    for(int Case=1;Case<=T;++Case){
        cin>>a>>b>>c>>d>>k;
        if(k==0){//没有gcd=0
            cout<<"Case "<<Case<<": 0\n";
            continue;
        }
        b/=k,d/=k;//将gcd=k化成gcd=1.
        ll ans1=0,ans2=0;
        if(b>d)swap(b,d);//b，d里b是最小值。
        for(int i=1;i<=b;++i){
        	ans1+=1ll*mu[i]*(b/i)*(d/i);
        	ans2+=1ll*mu[i]*(b/i)*(b/i);
        }
        cout<<"Case "<<Case<<": "<<ans1-(ans2>>1)<<endl;
    }
    return 0;
}

这样做是 O(n)。
例如，b=10000，d=11000，i>6000时， $⌊ \frac{b}{i} ⌋ = ⌊ \frac{d}{i} ⌋ = 1$ ，但是我们还是计算了4000次。

所以，不如先求 $μ$ 的前缀和，再求出对于b，d两数， $⌊ \frac{b}{i} ⌋, ⌊ \frac{d}{i} ⌋$ 均不变的范围，在该范围内，

$a n s = s u m (μ) \cdot ⌊ \frac{b}{i} ⌋ \cdot ⌊ \frac{d}{i} ⌋$

对于一个在区间[1 ，x]内的位置 i，设 $⌊ \frac{x}{i} ⌋ = k$ ，对于向下取整都为k的最后一个位置last，

$⌊ \frac{x}{l a s t} ⌋ = k$ ，且 $⌊ \frac{x}{l a s t + 1} ⌋ = k - 1 ， x m o d (l a s t + 1) = 0$

得： $l a s t = \frac{x}{⌊ \frac{x}{i} ⌋ - 1} - 1$ ,因为x和 $⌊ \frac{x}{i} ⌋ - 1$ 是整除关系，把-1去掉，分子相当于+1，所以 $\frac{x}{⌊ \frac{x}{i} ⌋} - 1$ 不是整除，

比原来少了1，再加上1就和原来一样了，所以为了方便，减1全部去掉。

最后再和另一个区间里的取最小值。

for(int i=1;i<=b;++i){
        	ans1+=1ll*mu[i]*(b/i)*(d/i);
        	ans2+=1ll*mu[i]*(b/i)*(b/i);
	}
//由上面变成了下面，只需处理一下前缀和：
//for(int i=1;i<MAXN;++i)sumu[i]=sumu[i-1]+mu[i];
for(int last,i=1;i<=b;i=last+1){
            last=min(b/(b/i),d/(d/i));
            ans1+=1ll*(sumu[last]-sumu[i-1])*(b/i)*(d/i);
            ans2+=1ll*(sumu[last]-sumu[i-1])*(b/i)*(b/i);
        }

复杂度： $O (\sqrt{n} + \sqrt{m})$

P2257 YY的GCD

题意：给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对。

还是一样，设 $f (d) = \sum_{i = 1}^{N} \sum_{j = 1}^{M} [g c d (i, j) = d]$ 假设 $N \leq M$

$F (n) = \sum_{i = 1}^{N} \sum_{j = 1}^{M} [n ∣ g c d (i, j)] = ⌊ \frac{N}{n} ⌋ ⌊ \frac{M}{n} ⌋$

由反演公式得到：

$f (n) = \sum_{n ∣ d} μ (\frac{d}{n}) F (d) = \sum_{n ∣ d} μ (\frac{d}{n}) ⌊ \frac{N}{d} ⌋ ⌊ \frac{M}{d} ⌋$

$a n s = \sum_{p \in p r i m e} f (p) = \sum_{p \in p r i m e} \sum_{p ∣ d} μ (\frac{d}{p}) ⌊ \frac{N}{d} ⌋ ⌊ \frac{M}{d} ⌋$

$= \sum_{p \in p r i m e} \sum_{k = 1} μ (k) ⌊ \frac{N}{d p} ⌋ ⌊ \frac{M}{k p} ⌋, d : = k p$

令 $T = k p$ ,有：

$a n s = \sum_{p \in p r i m e} \sum_{k = 1} μ (k) ⌊ \frac{N}{T} ⌋ ⌊ \frac{M}{T} ⌋ = \sum_{p \in p r i m e} \sum_{p ∣ T} μ (\frac{T}{p}) ⌊ \frac{N}{T} ⌋ ⌊ \frac{M}{T} ⌋$

重点来了，如何交换两个求和符号的顺序？画图大法好。

$p = p_{1}$ 时， $T_{1} = p_{1}, 2 p_{1}, 3 p_{1}, \dots, n_{1} p_{1}$
$p = p_{2}$ 时， $T_{2} = p_{2}, 2 p_{2}, 3 p_{2}, \dots, n_{2} p_{2}$
$\dots$
$p = p_{n}$ 时， $T_{n} = p_{n}, 2 p_{n}, 3 p_{n}, \dots, n p_{n}$
其中 $T_{i} \leq N$
所以，T从最小的质数 $p_{1}$ 开始枚举，由于任意合数 $x$ ，都能分解为 $x = k_{i} p_{i} = k_{j} p_{j}$ ,所以T遍历 $p_{1}$ 开始到x的

所有数， $\frac{T}{p}$ 枚举的是 $k_{i} \dots$ ，所以有：

$a n s = \sum_{p \in p r i m e} \sum_{p ∣ T} μ (\frac{T}{p}) ⌊ \frac{N}{T} ⌋ ⌊ \frac{M}{T} ⌋$

$= \sum_{T = 1}^{x} \sum_{p ∣ T, p \in p r i m e} μ (\frac{T}{p}) ⌊ \frac{N}{T} ⌋ ⌊ \frac{M}{T} ⌋$

$= \sum_{T = 1}^{x} ⌊ \frac{N}{T} ⌋ ⌊ \frac{M}{T} ⌋ \sum_{p ∣ T, p \in p r i m e} μ (\frac{T}{p})$ （T从2开始也是一样的，T=1时后面的求和等于0）

右边的求和可以通过预处理获得前缀和，对于 $\sum_{p ∣ T, p \in p r i m e} μ (\frac{T}{p})$ ，我们考虑每个质数p对T=kp的贡献。

#include<bits/stdc++.h>
using namespace std;
#define Init(arr,val) memset(arr,val,sizeof(arr))
const int inf=0x3f3f3f3f,mod=10007,MAXN=1e7+8;
typedef long long ll;
int mu[MAXN],prime[MAXN],cnt;
bool noprime[MAXN];
int sumu[MAXN];
void Mobius() {
   mu[1]=1;
   for(int i=2; i<MAXN; ++i) {
       if(!noprime[i]) {
   			prime[++cnt]=i;
   			mu[i]=-1;
       }
       for(int j=1; j<=cnt&&prime[j]*i<MAXN; ++j) {
           noprime[prime[j]*i]=1;
           if(i%prime[j]==0) {
               mu[i*prime[j]]=0;
               break;
           } else
               mu[i*prime[j]]=-mu[i];
       }
   }
   for(int i=1; i<=cnt; ++i)
       for(int j=1; j*prime[i]<MAXN; ++j)
           sumu[j*prime[i]]+=mu[j];
   for(int i=1; i<MAXN; ++i)
       sumu[i]+=sumu[i-1];
}
int main() {
   std::ios::sync_with_stdio(0);
   std::cin.tie(0);
   std::cout.tie(0);
   Mobius();
   int t;
   cin>>t;
   while(t--){
       int x,y,last;
       cin>>x>>y;
       ll ans=0;
       if(x>y)swap(x,y);
       for(int i=1;i<=x;i=last+1){
           last=min(x/(x/i),y/(y/i));
           ans+=1ll*(x/last)*(y/last)*(sumu[last]-sumu[i-1]);
       }
       cout<<ans<<endl;
   }
   return 0;
}