LuoguP3455 [POI2007]ZAP - Queries
本文同步发表于博主的Cmd Markdown Blog.
题目描述
给定(n,m,k),求以下式子的值:
[sum_{i=1}^n sum_{j=1}^m [gcd(i,j) = k]
]
Input & Output's examples
Input's examples
2
4 5 2
6 4 3
Output's examples
3
2
分析
比较显然的Mobius反演题目，没学过莫反的别尝试暴力了
我们来一起推一下式子(注意：以下的分数默认向下取整)。
[sum_{i=1}^n sum_{j=1}^m [gcd(i,j) = k]
]
同时除一下(k):
[sum_{i=1}^{frac nk} sum_{j=1}^{frac mk} [gcd(i,j) = 1]
]
由Mobius函数的定义式：
[sum_{d | n} mu(d) = [n = 1]
]
我们把(gcd(i,j))代入可得。
[sum_{d|gcd(i,j)} mu(d) = [gcd(i,j) = 1]
]
所以，原式可以简化为：
[sum_{i=1}^{frac nk} sum_{j=1}^{frac mk} sum_{d | gcd(i,j)} mu(d)
]
我们变换枚举顺序，首先枚举(d)，将它的求和符号提前，后面内容则变为布尔表达式
[sum_{d=1}^{min { frac nk,frac mk }} sum_{i=1}^{frac nk} sum_{j=1}^{frac mk} [d | gcd(i,j)] mu(d)
]
因为(d | gcd(i,j))，那么我们只需要枚举(d)的倍数即可，同时除一下(d):
[sum_{d=1}^{min { frac nk,frac mk }} sum_{i=1}^{frac {n}{kd}} sum_{j=1}^{frac {m}{kd}} [gcd(i,j) in Z] mu(d)
]
显然(gcd(i,j) in Z)恒成立，那么直接去掉。
[sum_{d=1}^{min { frac nk,frac mk }} sum_{i=1}^{frac {n}{kd}} sum_{j=1}^{frac {m}{kd}} mu(d)
]
(i,j)已经消失，那么我们直接根据乘法分配率：
[sum_{d=1}^{min { frac nk,frac mk }} lfloor frac {n}{kd}
floor lfloor frac {m}{kd}
floor mu(d)
]
因为(mu(d))是不完全积性函数，我们可以(O(n))求和。
但是本题要求多组数据,于是我们可以数论分块一下，然后预处理(mu)前缀和。
代码
/* Headers */
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cctype>
#include<algorithm>
#include<vector>
#include<queue>
#include<stack>
#include<climits>
#include<iostream>
#include<map>
#define FOR(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define ROF(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define FORL(i,a,b,c) for(long long i=(a);i<=(b);i+=(c))
#define ROFL(i,a,b,c) for(long long i=(a);i>=(b);i-=(c))
#define FORR(i,a,b,c) for(register int i=(a);i<=(b);i+=(c))
#define ROFR(i,a,b,c) for(register int i=(a);i>=(b);i-=(c))
#define lowbit(x) x&(-x)
#define LeftChild(x) x<<1
#define RightChild(x) (x<<1)+1
#define RevEdge(x) x^1
#define FILE_IN(x) freopen(x,"r",stdin);
#define FILE_OUT(x) freopen(x,"w",stdout);
#define CLOSE_IN() fclose(stdin);
#define CLOSE_OUT() fclose(stdout);
#define IOS(x) std::ios::sync_with_stdio(x)
#define Dividing() printf("-----------------------------------
");
namespace FastIO{
const int BUFSIZE = 1 << 20;
char ibuf[BUFSIZE],*is = ibuf,*its = ibuf;
char obuf[BUFSIZE],*os = obuf,*ot = obuf + BUFSIZE;
inline char getch(){
if(is == its)
its = (is = ibuf)+fread(ibuf,1,BUFSIZE,stdin);
return (is == its)?EOF:*is++;
}
inline int getint(){
int res = 0,neg = 0,ch = getch();
while(!(isdigit(ch) || ch == '-') && ch != EOF)
ch = getch();
if(ch == '-'){
neg = 1;ch = getch();
}
while(isdigit(ch)){
res = (res << 3) + (res << 1)+ (ch - '0');
ch = getch();
}
return neg?-res:res;
}
inline void flush(){
fwrite(obuf,1,os-obuf,stdout);
os = obuf;
}
inline void putch(char ch){
*os++ = ch;
if(os == ot) flush();
}
inline void putint(int res){
static char q[10];
if(res==0) putch('0');
else if(res < 0){putch('-');res = -res;}
int top = 0;
while(res){
q[top++] = res % 10 + '0';
res /= 10;
}
while(top--) putch(q[top]);
}
inline void space(bool x){
if(!x) putch('
');
else putch(' ');
}
}
inline void read(int &x){
int rt = FastIO::getint();
x = rt;
}
inline void print(int x,bool enter){
FastIO::putint(x);
FastIO::flush();
FastIO::space(enter);
}
/* definitions */
const int MAXN = 5e4 + 10;
int n,m,d,u[MAXN];
int T,prime[MAXN],cnt;
bool isprime[MAXN];
/* functions */
inline long long solve(){
long long ans = 0;
scanf("%d%d%d",&n,&m,&d);
if(n > m) std::swap(n,m);
n /= d, m /= d;
for(int i=1,j;i <= n;i = j + 1){
j = std::min(n / (n / i), m / (m / i));
ans += (long long) (u[j] - u[i - 1]) * (n / i) * (m / i);
}
return ans;
}
inline void init(){
u[1] = 1;
FOR(i,2,MAXN,1){
if(!isprime[i]) u[i] = -1, prime[++cnt] = i;
for(int j=1;j <= cnt && (long long) i * prime[j] <= MAXN;++j){
isprime[i * prime[j]] = true;
if(i % prime[j] == 0) break;
u[i * prime[j]] = -u[i];
}
u[i] += u[i-1];
}
}
int main(int argc,char *argv[]){
init();
scanf("%d",&T);
while(T--) printf("%lld
",solve());
return 0;
}
THE END
LuoguP3455 [POI2007]ZAP - Queries