[洛谷P2257] YY的GCD (莫比烏斯反演)

YY的GCD

題目描述

神犇YY虐完數論後給傻×kAc出了一題

給定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)爲質數的(x, y)有多少對

kAc這種傻×必然不會了，因而向你來請教……

多組輸入

輸入輸出格式

輸入格式：

第一行一個整數T 表述數據組數

接下來T行，每行兩個正整數，表示N, M

輸出格式：

T行，每行一個整數表示第i組數據的結果

輸入輸出樣例

輸入樣例#1： 複製

2
10 10
100 100

輸出樣例#1： 複製

30
2791

說明

\(T = 10000\)

\(N, M <= 10000000\)

Solution

根據shenlao告訴個人一個推題目賊爽的公式
\[\sum_{d|gcd(x,y)}\mu(d)=[gcd(x,y)=1]\]
原題意,求
\[Ans=\sum_{p\in prime}\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=p]\]
而後就開始上癮痛苦的推結論
\[Ans=\sum_{p\in prime}\sum_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{p}\rfloor}[gcd(i,j)=1]\]
\[Ans=\sum_{p\in prime}\sum_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{p}\rfloor}\sum_{d|gcd(i,j)}\mu(d)\]
更換枚舉項,由枚舉gcd(i,j)的約數改成直接枚舉d,這中間對於i有\(\frac{n}{p}\)個數,其中d的倍數的個數爲\(\frac {n}{d\times p}\),對於j同理,相乘以後就保證了\(d|gcd(i,j)\)這個條件
\[Ans=\sum_{p\in prime}\sum _{d=1}^{min({\lfloor\frac{n}{p}\rfloor},{\lfloor\frac{m}{p}\rfloor})}\mu(d){\lfloor\frac{n}{d\times p}\rfloor}{\lfloor\frac{m}{d\times p}\rfloor}\]
令\(d\times p=C\)
\[Ans=\sum_{C=1}^{min(n,m)}\sum_{p|C,p\in prime}\mu(\frac {C}{p}){\lfloor\frac{n}{C}\rfloor}{\lfloor\frac{m}{C}\rfloor}\]
\[Ans=\sum_{C=1}^{min(n,m)}{\lfloor\frac{n}{C}\rfloor}{\lfloor\frac{m}{C}\rfloor}\sum_{p|C,p\in prime}\mu(\frac {C}{p})\]
至此公式推導結束,\(O(n)\)就能夠過了,但是因爲多組詢問,因此咱們預處理一下.而後整除分塊\(O(\sqrt n)\)作,不會的能夠看本蒟蒻的博客整除分塊

Code

#include<bits/stdc++.h>
#define rg register
#define il inline
#define lol long long
#define NN 10000000
#define Min(a,b) (a)<(b)?(a):(b)

using namespace std;
const int N=1e7+10;

void in(int &ans) {
    ans=0; char i=getchar();
    while(i<'0' || i>'9') i=getchar();
    while(i>='0' && i<='9')  ans=(ans<<1)+(ans<<3)+i-'0',i=getchar();
}

int n,m,T,tot;
int mu[N],prime[N];
lol f[N];
bool vis[N];

void print(lol x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) print(x/10);
    putchar(x%10+'0');
}

il void init() {
    mu[1]=1;
    for(rg int i=2;i<=NN;i++) {
        if(!vis[i]) prime[++tot]=i,mu[i]=-1;
        for(rg int j=1;j<=tot && i*prime[j]<=NN;j++) {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0) break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(rg int i=1;i<=tot;i++)//注意枚舉質數在前
        for(rg int j=1;j*prime[i]<=NN;j++)//常數在後
            f[prime[i]*j]+=1ll*mu[j];//對後面的數有貢獻的是常數而不是質數
    for(rg int i=1;i<=NN;i++) f[i]+=f[i-1];
}

int main()
{
    in(T); init();
    while(T--) {
        in(n),in(m);
        lol ans=0; if(n>m) swap(n,m);
        for(rg int l=1,r;l<=n;l=r+1) {
            r=Min(n/(n/l),m/(m/l));
            ans+=1ll*(n/l)*(m/l)*(f[r]-f[l-1]);
        }
        print(ans);putchar('\n');
    }
    return 0;
}

博主蒟蒻,隨意轉載.但必須附上原文連接

http://www.cnblogs.com/real-l/