YY的GCD（洛谷-P2257）

最新推荐文章于 2021-04-06 17:22:23 发布

Alex_McAvoy

最新推荐文章于 2021-04-06 17:22:23 发布

阅读量687

点赞数

分类专栏： # 洛谷 # 数论——莫比乌斯反演

本文链接：https://blog.csdn.net/u011815404/article/details/81432265

版权

洛谷同时被 2 个专栏收录

161 篇文章

订阅专栏

数论——莫比乌斯反演

14 篇文章

订阅专栏

题目描述

神犇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

思路：

问题实质是求 $res=\sum_{i=1}^n\sum_{j=1}^m[gcd(x,y)=prime]$ 的值

设 f(k) 为满足 GCD(i,j)=k 的个数，即： $f(k)=\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=k]$

g(k) 为满足 GCD(i,j)=k 的倍数的个数，即： $g(k)=\sum_{k|d}f(d)=\left \lfloor \frac{n}{k} \right \rfloor \left \lfloor \frac{m}{k} \right \rfloor$

可以看出，g(k) 与 f(d) 间满足莫比乌斯反演的形式： $g(n)=\sum_{n|d}f(d) \Leftrightarrow f(n)=\sum_{n|d} u(\frac{d}{n})g(d)$

那么，对 $res=\sum_{i=1}^n\sum_{j=1}^m[gcd(x,y)=prime]$ 进行化简

那么有： $res=\sum_{p \in prime}\sum_{i=1}^n\sum_{j=1}^m[gcd(x,y)=p]$

将 f(k) 代入，得： $res=\sum_{p \in prime}f(p)$

根据莫比乌斯反演： $res=\sum_{p\in prime}\sum_{p|d} u(\frac{d}{p})g(d)=\sum_{p|d} u(\frac{d}{p})\left \lfloor \frac{n}{d} \right \rfloor\left \lfloor \frac{m}{d} \right \rfloor$

将枚举项 $\left \lfloor \frac{d}{p} \right \rfloor$ 更换为 d

那么有： $res=\sum_{p \in prime} \sum_{d=1}^{minn}u(d)g(d*p)=\sum_{p \in prime} \sum_{d=1}^{minn}u(d)\left \lfloor \frac{n}{d*p} \right \rfloor\left \lfloor \frac{m}{d*p} \right \rfloor$

设 T=d*p，那么有： $res=\sum_{T=1}^{min(n,m)}\sum_{t|T,t\in prime}u(\left \lfloor \frac{T}{t} \right \rfloor)\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor$

即： $res=\sum_{T=1}^{min(n,m)}\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor (\sum_{t|T,t\in prime}u(\left \lfloor \frac{T}{t} \right \rfloor))$

此时，如果是单组查询，可以直接做，时间复杂度为 O(n)，若为多组数据的话，此时需要加一个整除分块，因此要预处理所有 T=d*p 对应的 $\sum_{p|T}u(\left \lfloor \frac{T}{p} \right \rfloor)$ 的值，即设： $F(n)=\sum_{d|n}u (\frac{n}{d})$ ，sum(x) 是 F(x) 的前缀和

源代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<string>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<utility>
#include<stack>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<bitset>
#define PI acos(-1.0)
#define INF 0x3f3f3f3f
#define LL long long
#define Pair pair<int,int>
LL quickPow(LL a,LL b){ LL res=1; while(b){if(b&1)res*=a; a*=a; b>>=1;} return res; }
LL quickModPow(LL a,LL b,LL mod){ LL res=1; a=a%mod; while(b){if(b&1)res=(a*res)%mod; a=(a*a)%mod; b>>=1;} return res; }
LL getInv(LL a,LL mod){ return quickModPow(a,mod-2,mod); } // (a/b)%MOD=(a%MOD * getInv(b)%MOD)%MOD
const double EPS = 1E-10;
const int MOD = 1E9+7;
const int N = 10000000+5;
const int dx[] = {-1,1,0,0,-1,-1,1,1};
const int dy[] = {0,0,-1,1,-1,1,-1,1};
using namespace std;
int mu[N];
int prime[N],cnt;
bool bprime[N];
int F[N];
LL sum[N];

void getMu(int n){//线性筛求莫比乌斯函数

    cnt=0;
    mu[1]=1;//根据定义，μ(1)=1
    memset(bprime,false,sizeof(bprime));

    for(int i=2;i<=n;i++){//求2~N-1的莫比乌斯函数
        if(!bprime[i]){
            prime[++cnt]=i;//存储质数
            mu[i]=-1;//i为质数时，μ(1)=-1
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){//枚举i之前的素数个数
            bprime[i*prime[j]]=true;//不是质数
            if(i%prime[j])//i不是prime[j]的整数倍时，i*prime[j]就不会包含相同质因子
                mu[i*prime[j]]=-mu[i];//mu[k]=mu[i]*mu[prime[j]]，因为prime[j]是质数，mu值为-1
            else{
                mu[i*prime[j]]=0;
                break;//留到后面再筛
            }
        }
    }

    for(int j=1;j<=cnt;j++)
        for(int i=1;i*prime[j]<=n;i++)
            F[i*prime[j]]+=mu[i];
    for(int i=1;i<=n;i++)
        sum[i]=sum[i-1]+(LL)F[i];
}

int main() {
    getMu(10000000);

    int t;
    scanf("%d",&t);
    while(t--){
        int n,m;
        scanf("%d%d",&n,&m);

        int minn=min(n,m);
        LL res=0;
        for(int left=1,right;left<=minn;left=right+1){//整除分块
            right=min(n/(n/left),m/(m/left));
            res+=1ll*(n/left)*(m/left)*(sum[right]-sum[left-1]);
        }
        printf("%lld\n",res);
    }
    return 0;
}