莫比乌斯反演

莫比乌斯函数定义：

设 (n = p_1 ^ {k_1} cdot p_2 ^ {k_2} cdotcdotscdot p_m ^ {k_m})，其中 p 为素数，则定义如下：

[mu(n) = egin{cases} 1 & n = 1 \ (-1) ^ m & prodlimits_{i = 1} ^ {m} k_i = 1 \ 0 & extrm{otherwise}(k_i gt 1) end{cases} ]

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

1. 由于Mobius函数是积性函数，即 (mu(a cdot b)=mu(a)cdot mu(b), (a,b)=1),所以
   可以在(O(n))的时间内筛出([1,n])的函数值，一个数 (x= p_1 ^ {k_1} cdot p_2 ^ {k_2} cdotdotscdot p_m ^ {k_m})，
   若(k_i=1 , k_j>1, j>i) 则x只会被(frac x {p_j} mod 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多了一个质因子。
        }
    }
}

2. (sum_{dmid n}mu(d) = [n = 1])

证明：当 n = 1 时有 (mu(1) = 1)，显然成立。否则设 (n = p_1 ^ {k_1} cdot p_2 ^ {k_2}cdotcdotscdot p_m ^ {k_m})​​ ，(d = p_1 ^ {x_1} cdot p_2 ^ {x_2} cdotcdotscdot p_m ^ {x_m})

​​ 根据(mu(n)) 的定义，只需考虑 (x_i = 0) 或 (x_i = 1)的情况。我们设 d 中存在 r 个 (x_i)为 1，那么有：

(sum_{dmid n}mu(d) = sum_{r = 0}^minom m r(-1)^r(n eq 1))

根据二项式定理：

((x+y)^{n}=sum _{k=0}^{n}inom n kx^{n-k}y^{k} (x+y))
​​
我们令 x = 1, y = -1，即得证：

(sum_{r = 0}^minom m r(-1)^r = (1 - 1)^m = 0(n eq 1))

莫比乌斯反演

假设对于数论函数 (f(n)) 和(F(n))，满足以下关系式：(F(n)=sum _{{d|n}}f(d))

则将其默比乌斯反演公式定义为：

[f(n)=sum _{d|n}mu (d) F({frac {n}{d}}) ]

或者满足：(F(n)=sum_{n|d}f(d))

有

[f(n)=sum_{n|d}mu(frac d n)F(d) ]

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

[sum_{i=1}^{a}sum_{j=1}^{b}[gcd(i,j)=k] ]

即求

[sum_{i=1}^{lfloorfrac a k floor}sum_{j=1}^{lfloorfrac b k floor}[gcd(i,j)=1] ]

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

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

所以有反演公式：(f(n)=sum_{n|d}mu(frac d n)F(d))

显然，(F(n)=lfloorfrac a n floorcdotlfloorfrac b n floor)

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

所以(f(n)=sum_{n|d}mu(frac d n)F(d)=sum_{n|d}mu(frac d n)cdotlfloorfrac a d floorcdotlfloorfrac b d floor)

由于k=1

有(f(1)=sum_{1|d}mu(d)F(d)=sum_{1|d}mu(d)cdotlfloorfrac a d floorcdotlfloorfrac b d floor)

HDU-1695

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

(ans=f(1)=sum_{1|i}^{min(b,d)}mu(i)F(i)=sum_{1|i}mu(i)cdotlfloorfrac b i floorcdotlfloorfrac d i floor)
为何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
";
            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时，(lfloorfrac b i floor=lfloorfrac d i floor=1)，但是我们还是计算了4000次。

所以，不如先求(mu) 的前缀和，再求出对于b，d两数，(lfloorfrac b i floor, lfloorfrac d i floor)均不变的范围，在该范围内，

(ans=sum(mu)cdotlfloorfrac b i floorcdotlfloorfrac d i floor)

对于一个在区间[1 ，x]内的位置 i，设 (lfloorfrac x i floor=k) ，对于向下取整都为k的最后一个位置last，

(lfloorfrac x {last} floor=k)，且 (lfloorfrac {x} {last+1} floor=k-1，xmod ({last+1})=0)

得：(last=frac{x}{lfloorfrac x i floor -1}-1),因为x和(lfloorfrac x i floor -1)是整除关系，把-1去掉，分子相当于+1，所以(frac{x}{lfloorfrac x i floor }-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[gcd(i,j)=d]) 假设(Nle M)

(F(n)=sum_{i=1}^N sum_{j=1}^M[n|gcd(i,j)]=lfloorfrac N n floorlfloorfrac M n floor)

由反演公式得到：

(f(n)=sum_{n|d}mu(frac d n)F(d)=sum_{n|d}mu(frac d n)lfloorfrac N d floorlfloorfrac M d floor)

(ans=sum_{pin prime}f(p)=sum_{pin prime}sum_{p|d}mu(frac d p)lfloorfrac N d floorlfloorfrac M d floor)

(=sum_{pin prime}sum_{k=1}mu(k)lfloorfrac N {dp} floorlfloorfrac M {kp} floor, d:=kp)

令(T=kp),有：

(ans=sum_{pin prime}sum_{k=1}mu(k)lfloorfrac N {T} floorlfloorfrac M {T} floor=sum_{pin prime}sum_{p|T}mu(frac T p)lfloorfrac N {T} floorlfloorfrac M {T} floor)

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

(p=p_1)时，(T_1=p_1,2p_1,3p_1,dots,n_1p_1)
(p=p_2)时，(T_2=p_2,2p_2,3p_2,dots,n_2p_2)
(dots)
(p=p_n)时，(T_n=p_n,2p_n,3p_n,dots,np_n)
其中(T_ile N)
所以，T从最小的质数(p_1)开始枚举，由于任意合数 (x)，都能分解为(x=k_ip_i=k_jp_j),所以T遍历(p_1)开始到x的

所有数，(frac T p)枚举的是(k_i dots)，所以有：

(ans=sum_{pin prime}sum_{p|T}mu(frac T p)lfloorfrac N {T} floorlfloorfrac M {T} floor)

(=sum_{T=1}^{x}sum_{p|T, pin prime}mu(frac T p)lfloorfrac N {T} floorlfloorfrac M {T} floor)

(=sum_{T=1}^{x}lfloorfrac N {T} floorlfloorfrac M {T} floorsum_{p|T, pin prime}mu(frac T p))（T从2开始也是一样的，T=1时后面的求和等于0）

右边的求和可以通过预处理获得前缀和，对于(sum_{p|T, pin prime}mu(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;
}