YY的GCD

题目描述

神犇YY虐完数论后给傻×kAc出了一题html

给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对c++

kAc这种傻×必然不会了，因而向你来请教……spa

多组输入code

输入输出格式

输入格式：htm

第一行一个整数T 表述数据组数blog

接下来T行，每行两个正整数，表示N, Mget

输出格式：博客

T行，每行一个整数表示第i组数据的结果it

输入输出样例

输入样例#1： 复制io

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/