[HAOI2011]Problem b

标签：http   prime   href   prepare   oid   space   printf   int   pac   

[HAOI2011]Problem b

给出n个询问，询问\(\sum_{i=a}^b\sum_{j=c}^d(gcd(i,j)==k)\),100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000。

解

法一:Mobius反演

设
\[b\leq d,f(k)=\sum_{i=a}^b\sum_{j=c}^d(gcd(i,j)==k)\]
\[F(k)=\sum_{i=a}^b\sum_{j=c}^d(k|gcd(i,j))=\]
\[([b/k]-[(a-1)/k])([d/k]-[(c-1)/k])\]

由Mobius反演定理有
\[ans=f(k)=\sum_{k|x}F(x)\mu(x/k)=\]
\[\sum_{k|x}([b/x]-[(a-1)/x])([d/x]-[(c-1)/x])\mu(x/k)=\]
\[\sum_{x=1}^{[b/k]}([\frac{b}{xk}]-[\frac{a-1}{kx}])([\frac{d}{xk}]-[\frac{c-1}{xk}])\mu(x)\]

利用式子进行整除分块即可。

参考代码：

#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[50001];
int prime[50001],pt,mb[50001];
il int min(int,int);
il void prepare(int);
int main(){
    int a,b,c,d,k,i,j,ans(0),lsy;
    prepare(50000),scanf("%d",&lsy);
    while(lsy--){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        b/=k,d/=k,(--a)/=k,(--c)/=k,ans&=0;
        if(b>d)swap(b,d),swap(a,c);
        for(i=1;i<=b;i=j+1){
            j=min(b/(b/i),d/(d/i));
            if(a/i)j=min(j,a/(a/i));
            if(c/i)j=min(j,c/(c/i));
            ans+=(b/i-a/i)*(d/i-c/i)*(mb[j]-mb[i-1]);
        }printf("%d\n",ans);
    }
    return 0;
}
il int min(int a,int b){
    return a<b?a:b;
}
il void prepare(int n){
    int i,j;check[1]|=mb[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])prime[++pt]=i,mb[i]=-1;
        for(j=1;j<=pt&&prime[j]<=n/i;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mb[i*prime[j]]=-mb[i];
        }
    }for(i=1;i<=n;++i)mb[i]+=mb[i-1];
}

法二：容斥原理

注意到下标的开始不为1，而我们能够做的是下标从1开始的问题，于是设\((a,b)\)表示a，b范围内，下界为1的满足题意的gcd限制的方案数，不难由容斥原理得

\[ans=(c,d)-(a-1,d)-(b,c-1)+(a-1,c-1)\]

写出单个无下界限制的函数，套用公式即可。

参考代码：

#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[50001];
int prime[50001],pt,mb[50001];
il void read(int&);
il int min(int,int);
il ll mobius(int,int,int);
void prepare(int),pen(ll);
int main(){
    int lsy,a,b,c,d,p;
    read(lsy),prepare(50000);while(lsy--){
        read(a),read(b),read(c),read(d),read(p);
        pen(mobius(b,d,p)-mobius(a-1,d,p)-
            mobius(b,c-1,p)+mobius(a-1,c-1,p)),putchar('\n');
    }
    return 0;
}
il ll mobius(int a,int b,int p){
    int i,j;ll ans(0);a/=p,b/=p;
    if(a>b)swap(a,b);
    for(i=1;i<=a;i=j+1)
        j=min(a/(a/i),b/(b/i)),
            ans+=(ll)(a/i)*(b/i)*(mb[j]-mb[i-1]);
    return ans;
}
void prepare(int n){
    int i,j;check[1]|=mb[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])prime[++pt]=i,mb[i]=-1;
        for(j=1;i*prime[j]<=n&&j<=pt;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mb[i*prime[j]]=-mb[i];
        }
    }for(i=1;i<=n;++i)mb[i]+=mb[i-1];
}
il int min(int a,int b){
    return a<b?a:b;
}
il void read(int &x){
    x&=0;ri char c;while(c=getchar(),c<'0'||c>'9');
    while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
}
void pen(ll x){
    if(x>9)pen(x/10);putchar(x%10+48);
}

[HAOI2011]Problem b

标签：http   prime   href   prepare   oid   space   printf   int   pac   

原文地址：https://www.cnblogs.com/a1b3c7d9/p/10781665.html