洛谷 P2257 YY的GCD

YY的GCD

题目描述

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

给定\(N\), \(M\) ,求\(1 \le x \le N,1 \le y \le M\)且\(gcd(x, y)\)为质数的\((x, y)\)有多少对

kAc这种傻×必然不会了，于是向你来请教……

多组输入

输入输出格式

输入格式：

第一行一个整数\(T\)表述数据组数

接下来T行，每行两个正整数，表示\(N,M\)

输出格式：

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

说明

\(T = 10000\)
\(N, M \le 10000000\)

Solution

思路：按反演的套路把式子列出来以后，通过想办法操作求和，把好算的放在一起，\(\mu\)一个人呆。

以下\(f(x)=\sum_{i=1}^a\sum_{j=1}^b[gcd(a,b)=x]\)，\(g(x)=\lfloor\frac{a}{x}\rfloor \lfloor \frac{b}{x} \rfloor\)，\(p\)代表质数集合

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

\[=\sum_{n \in p}f(n) \]

\[=\sum_{n \in p}\sum_{n|d}^{min(a,b)}\mu(\frac{d}{n})g(d) \]

\[=\sum_{n \in p}\sum_{k=1}^{\min(\lfloor \frac{a}{n} \rfloor, \lfloor \frac{b}{n} \rfloor)}\mu(k)\lfloor\frac{a}{kn}\rfloor \lfloor\frac{b}{kn}\rfloor \]

现在我们不想让\(\mu\)和别人呆在一起，因为两个求和之间也没夹个什么东西，不妨把它们换一换，直接换肯定不行，得找一波操作。

令\(T=kn\)，则有\(k=\frac{T}{n}\)，上式就等于

\[=\sum_{T=1}^{min(a,b)}\sum_{n|T,n\in p}\mu(\frac{T}{n})\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor \]

\[=\sum_{T=1}^{min(a,b)}\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor\sum_{n|T,n\in p}\mu(\frac{T}{n}) \]

后面那个感觉预处理还挺慢的，但根据调和级数不会慢于\(O(nlnn)\)，而且跑一下发现才\(3e7\)多运算。

预处理完后面那个就可以除法分块回答询问了。

---

Code:

#include <cstdio>
#define ll long long
const int N=1e7;
int mu[N+10],f[N+10],pri[N+10],ispri[N+10],cnt,a,b,T;
void init()
{
    mu[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!ispri[i])
        {
            mu[i]=-1;
            pri[++cnt]=i;
        }
        for(int j=1;j<=cnt&&pri[j]*i<=N;j++)
        {
            ispri[pri[j]*i]=1;
            if(i%pri[j]==0) break;
            else mu[pri[j]*i]=-mu[i];
        }
    }
    for(int i=1;i<=cnt;i++)
        for(int j=pri[i];j<=N;j+=pri[i])
            f[j]+=mu[j/pri[i]];
    for(int i=1;i<=N;i++)
        f[i]+=f[i-1];
}
int min(int x,int y){return x<y?x:y;}
int main()
{
    init();
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d",&a,&b);
        ll ans=0;
        for(int l=1,r;l<=min(a,b);l=r+1)
        {
            r=min(a/(a/l),b/(b/l));
            ans+=1ll*(a/l)*(b/l)*(f[r]-f[l-1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}

---

2018.10.19