洛谷P2257 YY的GCD 题解

洛谷P2257 YY的GCD 题解

题意

多组询问。

给定 N , M N, M N,M,求 1 ≤ x ≤ N 1 \leq x \leq N 1xN 1 ≤ y ≤ M 1 \leq y \leq M 1yM gcd ⁡ ( x , y ) \gcd(x,y) gcd(x,y) 为质数的 ( x , y ) (x,y) (x,y) 有多少对。

一句话题意

P \mathcal{P} P 为所有素数组成的集合,求
∑ i = 1 n ∑ j = 1 m [ gcd ⁡ ( i , j ) ∈ P ] \sum_{i=1}^{n}\sum_{j=1}^{m}\left[\gcd(i,j)\in\mathcal{P}\right] i=1nj=1m[gcd(i,j)P]
不妨假设 n ≤ m n\le m nm

然后么推推柿子

友情提醒,这道题涉及经典操作:转化枚举量、构造莫比乌斯反演

不知道的可以先去看看别的题,比如这个题 link
∑ i = 1 n ∑ j = 1 m [ gcd ⁡ ( i , j ) ∈ P ] = ∑ k = 1 n ∑ i = 1 n ∑ j = 1 m [ gcd ⁡ ( i , j ) = k ] [ k ∈ P ] = ∑ k = 1 n ∑ i = 1 ⌊ n k ⌋ ∑ j = 1 ⌊ m k ⌋ [ gcd ⁡ ( i , j ) = 1 ] [ k ∈ P ] = ∑ k = 1 n ∑ i = 1 ⌊ n k ⌋ ∑ j = 1 ⌊ m k ⌋ ∑ d ∣ gcd ⁡ ( i , j ) μ ( d ) [ k ∈ P ] = ∑ k = 1 n ∑ d = 1 ⌊ n k ⌋ μ ( d ) ⌊ n k d ⌋ ⌊ m k d ⌋ [ k ∈ P ] = ∑ l = 1 n ⌊ n l ⌋ ⌊ m l ⌋ ∑ k ∣ l ∧ k ∈ P μ ( l k ) \begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{m}\left[\gcd(i,j)\in\mathcal{P}\right] \\= &\sum_{k=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^m\left[\gcd(i,j)=k\right][k \in \mathcal{P}] \\= &\sum_{k=1}^{n}\sum_{i=1}^{\left\lfloor{\frac{n}{k}}\right\rfloor}\sum_{j=1}^{\left\lfloor{\frac{m}{k}}\right\rfloor}[\gcd(i,j)=1][k \in \mathcal{P}] \\=&\sum_{k=1}^{n}\sum_{i=1}^{\left\lfloor{\frac{n}{k}}\right\rfloor}\sum_{j=1}^{\left\lfloor{\frac{m}{k}}\right\rfloor}\sum_{d\mid \gcd(i,j)}\mu(d)[k \in \mathcal{P}] \\=&\sum_{k=1}^{n}\sum_{d=1}^{\left\lfloor{\frac{n}{k}}\right\rfloor}\mu(d)\left\lfloor{\dfrac{n}{kd}}\right\rfloor\left\lfloor{\dfrac{m}{kd}}\right\rfloor[k \in \mathcal{P}] \\=&\sum_{l=1}^{n}\left\lfloor{\dfrac{n}{l}}\right\rfloor\left\lfloor{\dfrac{m}{l}}\right\rfloor\sum_{k \mid l\land k \in \mathcal{P}}\mu\left(\dfrac{l}{k}\right) \end{aligned} =====i=1nj=1m[gcd(i,j)P]k=1ni=1nj=1m[gcd(i,j)=k][kP]k=1ni=1knj=1km[gcd(i,j)=1][kP]k=1ni=1knj=1kmdgcd(i,j)μ(d)[kP]k=1nd=1knμ(d)kdnkdm[kP]l=1nlnlmklkPμ(kl)
注: ∧ \land 表示逻辑与,可认为其等价于 “且”

f ( x ) = ∑ k ∣ x ∧ k ∈ P μ ( x k ) f(x) = \sum\limits_{k \mid x\land k \in \mathcal{P}}\mu\left(\dfrac{x}{k}\right) f(x)=kxkPμ(kx)

然后这个题其实到这里就差不多了

前面的用数论分块

后面那个 μ \mu μ 可以在筛法里面求

时间复杂度 O ( n log ⁡ log ⁡ n + Q n ) O(n\log \log n + Q\sqrt{n}) O(nloglogn+Qn )

#include <bits/stdc++.h>
using namespace std;
// #define int long long
#define INF 0x3f3f3f3f3f3f3f3f
namespace FastIO
{
    #define gc() readchar()
    #define pc(a) putchar(a)
    #define SIZ (int)(1e6+15)
    char buf1[SIZ],*p1,*p2;
    char readchar()
    {
        if(p1==p2)p1=buf1,p2=buf1+fread(buf1,1,SIZ,stdin);
        return p1==p2?EOF:*p1++;
    }
    template<typename T>void read(T &k)
    {
        char ch=gc();T x=0,f=1;
        while(!isdigit(ch)){if(ch=='-')f=-1;ch=gc();}
        while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=gc();}
        k=x*f;
    }
    template<typename T>void write(T k)
    {
        if(k<0){k=-k;pc('-');}
        static T stk[66];T top=0;
        do{stk[top++]=k%10,k/=10;}while(k);
        while(top){pc(stk[--top]+'0');}
    }
}using namespace FastIO;
#define N (int)(1e7+15)
int Q,x,y;
int prime[6000005],mu[N],pcnt,sum[N];
bool ck[N];
void Mobius()
{
    mu[1]=1;
    for(int i=2; i<=N-5; i++)
    {
        if(!ck[i])
        {
            prime[++pcnt]=i;
            mu[i]=-1;
        }
        for(int j=1; j<=pcnt&&i*prime[j]<=N-5; j++)
        {
            int pos=i*prime[j];
            ck[pos]=1;

            if(i%prime[j])
            {
                mu[pos]=-mu[i];
            }else
            {
                mu[pos]=0;
                break;
            }
        }
    }
    for(int i=1; i<=pcnt; i++)
        for(int j=1; prime[i]*j<=N-5; j++)
            sum[j*prime[i]]+=mu[j];
            
    for(int i=1; i<=N-5; i++)
        sum[i]+=sum[i-1];
}
long long solve(int n,int m)
{
    long long res=0;
    for(int l=1,r; l<=min(n,m); l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        res+=(long long)(sum[r]-sum[l-1])*(n/l)*(m/l);
    }
    return res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    read(Q);
    Mobius();
    while(Q--)
    {
        read(x);read(y);
        write(solve(x,y));pc('\n');
    }
    return 0;
}

这份代码很慢,跑了7.74s, 目前最优解就跑了 5.69s(截止20220511)

其实这个 O ( n log ⁡ log ⁡ n ) O(n \log\log n) O(nloglogn) 的枚举并没有必要

考虑在线性筛的过程中更新答案
f ( x ) = ∑ k ∣ x ∧ k ∈ P μ ( x k ) f(x) = \sum\limits_{k \mid x\land k \in \mathcal{P}}\mu\left(\dfrac{x}{k}\right) f(x)=kxkPμ(kx)
x = a b x=ab x=ab ,其中 a a a x x x 的最小质因数

p p p 为当前线性筛枚举的一个用于标记的素数(prime[j]

  • x ∈ P x \in \mathcal{P} xP ,则显然 f ( x ) = μ ( 1 ) = 1 f(x) = \mu(1) = 1 f(x)=μ(1)=1

  • x ∉ P ∧ b   m o d   a = 0 x\notin \mathcal P \land b \bmod a = 0 x/Pbmoda=0

    则有 x = a p b ′ , b = a p − 1 b ′ , p x = a^p b^{\prime},b=a^{p-1}b^{\prime},p x=apb,b=ap1b,p 为大于 1 1 1 的正整数 (即 x x x 有多个最小质因数 )

    • b b b 没有多个最小质因数,则当且仅当 p = a p=a p=a μ ( x p ) = μ ( b ) ≠ 0 \mu\left(\dfrac{x}{p}\right) = \mu(b) \ne 0 μ(px)=μ(b)=0

      因此有 f ( x ) = μ ( b ) f(x) = \mu(b) f(x)=μ(b)

    • b b b 含有多个最小质因数,则 ∀ p ∈ P ′ , μ ( x p ) = 0 , P ′ \forall p \in \mathcal{P}^{\prime},\mu\left(\dfrac{x}{p}\right)=0,\mathcal{P}^{\prime} pP,μ(px)=0,P 为当前筛出的素数集

      不过我们无须单独考虑这种情况,因为此时仍然有 μ ( x p ) = μ ( b ) \mu\left(\dfrac{x}{p}\right) = \mu(b) μ(px)=μ(b)

  • x ∉ P ∧ b   m o d   a ≠ 0 x\notin \mathcal{P} \land b \bmod a \ne 0 x/Pbmoda=0

    则此时 x x x 有唯一且最小的质因数 a a a

    显然有 μ ( x ) = − μ ( b ) \mu(x) = -\mu(b) μ(x)=μ(b) ,因为它仅仅多了一个本质不同的质因子

    则对于 μ ( x p ) = − μ ( b p ) \mu\left(\dfrac{x}{p}\right) = -\mu\left(\dfrac{b}{p}\right) μ(px)=μ(pb) 也是成立的

    因此 f ( b ) f(b) f(b) 中的每一项在 f ( x ) f(x) f(x) 中都能找到对应的

    注意到
    f ( x ) = f ( a b ) = ∑ k ∣ a b ∧ k ∈ P μ ( a b k ) f(x) = f(ab) =\sum\limits_{k \mid ab\land k \in \mathcal{P}}\mu\left(\dfrac{ab}{k}\right) f(x)=f(ab)=kabkPμ(kab)
    由于 gcd ⁡ ( a , b ) = 1 \gcd(a,b)=1 gcd(a,b)=1 (显然),则 f ( x ) f(x) f(x) 仅仅比 f ( b ) f(b) f(b) 多了一项 μ ( a b a ) = μ ( b ) \mu\left(\dfrac{ab}{a}\right) = \mu(b) μ(aab)=μ(b)

    则有 f ( x ) = − f ( b ) + μ ( b ) f(x) = -f(b)+\mu(b) f(x)=f(b)+μ(b)

综上,有
f ( x ) = { μ ( 1 ) , x ∈ P , μ ( b ) , x ∉ P ∧ b   m o d   a = 0 , − f ( b ) + μ ( b ) , x ∉ P ∧ b   m o d   a ≠ 0. f(x) = \begin{cases} \mu(1),&x \in \mathcal{P}, \\\\\mu(b),&x \notin \mathcal{P} \land b \bmod a=0, \\\\-f(b)+\mu(b),&x\notin \mathcal{P} \land b\bmod a \ne 0. \end{cases} f(x)=μ(1),μ(b),f(b)+μ(b),xP,x/Pbmoda=0,x/Pbmoda=0.
于是我们就可以在 O ( n + Q n ) O(n + Q\sqrt{n}) O(n+Qn ) 的时间复杂度内解决这个问题了

然后稍微卡一卡常,就跑的飞快

对了,目前最优解其实是我的代码哦(虽然没啥用

代码如下

#include <bits/stdc++.h>
using namespace std;
// #define int long long
typedef long long ll;
#define INF 0x3f3f3f3f3f3f3f3f
namespace FastIO
{
    #define gc() readchar()
    #define pc(a) putchar(a)
    #define SIZ (int)(1e6+15)
    char buf1[SIZ],*p1,*p2;
    char readchar()
    {
        if(p1==p2)p1=buf1,p2=buf1+fread(buf1,1,SIZ,stdin);
        return p1==p2?EOF:*p1++;
    }
    template<typename T>void read(T &k)
    {
        char ch=gc();T x=0,f=1;
        while(!isdigit(ch)){if(ch=='-')f=-1;ch=gc();}
        while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=gc();}
        k=x*f;
    }
    template<typename T>void write(T k)
    {
        if(k<0){k=-k;pc('-');}
        static T stk[66];T top=0;
        do{stk[top++]=k%10,k/=10;}while(k);
        while(top){pc(stk[--top]+'0');}
    }
}using namespace FastIO;
#define N (int)(1e7+15)
int Q,mx,x[10005],y[10005];;
int prime[1000005],mu[N],pcnt,f[N];
bool ck[N];
void Mobius()
{
    mu[1]=1;
    for(int i=2; i<=mx+5; i++)
    {
        if(!ck[i])
        {
            prime[++pcnt]=i;
            mu[i]=-1;
            f[i]=1;
        }
        for(int j=1; j<=pcnt&&i*prime[j]<=mx+5; j++)
        {
            int pos=i*prime[j];
            ck[pos]=1;
            if(i%prime[j])
            {
                f[pos]=-f[i]+mu[i];
                mu[pos]=-mu[i];
            }else
            {
                f[pos]=mu[i];
                mu[pos]=0;
                break;
            }
        }
        f[i]+=f[i-1];
    }
}
ll solve(int n,int m)
{
    ll res=0;
    for(int l=1,r; l<=min(n,m); l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        res+=(ll)(f[r]-f[l-1])*(n/l)*(m/l);
    }
    return res;
}

signed main()
{
    read(Q);
    for(int i=1; i<=Q; i++)
        read(x[i]),read(y[i]);
    mx=max(*max_element(x+1,x+1+Q),*max_element(y+1,y+1+Q));
    Mobius();
    for(int i=1; i<=Q; i++)
        write(solve(x[i],y[i])),pc('\n');
    return 0;
}

如果有什么错误或者问题的话,欢迎留言

参考文献

[1] https://www.luogu.com.cn/blog/An-Amazing-Blog/solution-p2257

[2] https://siyuan.blog.luogu.org/solution-p2257

转载请说明出处

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值