ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

[POI2007]ZAP-Queries

2019-04-28 08:52:59  阅读:328  来源: 互联网

标签:int POI2007 ans mu read il Queries sum ZAP


[POI2007]ZAP-Queries

求n组\(\sum_{i=1}^a\sum_{j=1}^b(gcd(i,j)==d)\), 1≤n≤50 000,1≤d≤a,b≤50 000。

不难看出是约数组合计数问题,而解决该问题常用思路有容斥原理和Mobius反演。

法一:容斥原理

尽可能特殊化问题,令\(a/=d,b/=d\),原问题等价于在这个范围里寻找互质的数的对数,而容斥原理关键在于构造不合理情况,即有别的约数,于是设\(D(a,b,i)\)表示a,b范围里有公约数i的个数,不难得知\(D(a,b,i)=[a/i][b/i]\),于是由容斥原理有
\[ans=\sum_{i=1}^{min(a,b)}\mu(i)D(a,b,i)=\sum_{i=1}^{min(a,b)}\mu(i)[a/i][b/i]\]
预处理\(\mu\)函数前缀和,进行整出分块即可。

法二:Mobius反演


\[f(d)=\sum_{i=1}^a\sum_{j=1}^b(gcd(i,j)==d)\]
\[F(d)=\sum_{i=1}^a\sum_{j=1}^b(d|gcd(i,j))=[a/d][b/d]\]
由Mobius反演定理我们有
\[f(d)=\sum_{d|x}F(x)\mu(x/d)=\sum_{d|x}[a/x][b/x]\mu(x/d)\]
于是
\[ans=f(d)=\sum_{d|x}[a/x][b/x]\mu(x/d)=\]
\[\sum_{x=1}^{min([a/d],[b/d)]}[a/x][b/x]\mu(x/d)\]

对式子进行整除分块即可。

参考代码

#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[5137],pt,mb[50001];
il void read(int&);
il int min(int,int);
void pen(int),prepare(int);
int main(){
    int lsy,a,b,n,i,j,ans;read(lsy);
    prepare(50000);while(lsy--){
        read(a),read(b),read(n),ans&=0;
        a/=n,b/=n;if(a>b)swap(a,b);
        for(i=1;i<=a;i=j+1)
            j=min(a/(a/i),b/(b/i)),
                ans+=(a/i)*(b/i)*(mb[j]-mb[i-1]);
        pen(ans),putchar('\n');
    }
    return 0;
}
void prepare(int n){
    int i,j;check[1]|=mb[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])mb[i]=-1,prime[++pt]=i;
        for(j=1;j<=pt&&i*prime[j]<=n;++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;
}
void pen(int x){
    if(x>9)pen(x/10);putchar(x%10+48);
}
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();
}

标签:int,POI2007,ans,mu,read,il,Queries,sum,ZAP
来源: https://www.cnblogs.com/a1b3c7d9/p/10781632.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有