[HAOI2011]Problem b

description

求

\[\sum_{i=1}^{n} \sum_{j=1}^{m}[gcd(i, j)==k] \]

solution

\[\begin{align} \sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i, j)==k] &= \sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[gcd(i, j)==1]\\ &=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d|gcd(i, j)}\mu(d)\\ &=\sum_{d=1}\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}[d|i]\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[d|j]\\ &=\sum_{d=1}\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor \end{align} \]

\[\sum_{i=a}^{n}\sum_{j=b}^{m}[gcd(i, j)==k]=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i, j)==k]-\sum_{i=1}^{a-1}\sum_{j=1}^{m}[gcd(i, j)==k]-\sum_{i=1}^{n}\sum_{j=1}^{b-1}[gcd(i, j)==k]+\sum_{i=1}^{a-1}\sum_{j=1}^{b-1}[gcd(i, j)==k] \]

code

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<map>
#include<bitset>
#include<queue>
#include<set>
#define R register
#define mod 998244353
#define debug puts("wn")
#define next MabLcdG
#define pos(i, j) (i*m+j)
using namespace std;
typedef long long ll; 
typedef long double ld;
typedef double dl;
typedef unsigned long long ull;
template<typename T>void read(T &x);
template<typename T>void write(T x);
template<typename T>void writesp(T x);
template<typename T>void writeln(T x);

const ll N=5e4+5, lim=5e4;

ll mu[N], prime[N], tot;
bool flag[N];
inline void calc_mu(){
	flag[1]=true; mu[1]=1;
	for(R ll i=2; i<=lim; i++){
		if(!flag[i]){
			prime[++tot]=i;
			mu[i]=-1;
		}
		for(R ll j=1; j<=tot && i*prime[j]<=lim; j++){
			flag[i*prime[j]]=true;
			if(i%prime[j]==0){
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(R ll i=1; i<=lim; i++) mu[i]+=mu[i-1];
}

ll T;

inline ll solve(ll n, ll m){
	ll res=0;
	for(R ll l=1, r; l<=n && l<=m; l=r+1){
		r=min(n/(n/l), m/(m/l));
		res+=(mu[r]-mu[l-1])*(n/l)*(m/l);
	} 
	return res;
}
int main(){
	calc_mu();
	read(T);
	ll a, b, c, d, k;
	while(T--){
		read(a); read(b); read(c); read(d); read(k);
		writeln(solve(b/k, d/k)-solve((a-1)/k, d/k)-solve(b/k, (c-1)/k)+solve((a-1)/k, (c-1)/k));
	}
}

template<typename T>void read(T &x){
	char wn=getchar(); ll t=1; x=0;
	while(wn<'0' || wn>'9'){
		if(wn=='-') t=-1;
		wn=getchar();
	}
	while(wn>='0' && wn<='9'){
		x=x*10+wn-'0'; wn=getchar();
	}
	x*=t;
}

template<typename T>void write(T x){
	if(x<0){putchar('-'); x=-x;}
	if(x<=9){putchar(x+'0'); return;}
	write(x/10); putchar(x%10+'0');
}

template<typename T>void writeln(T x){
	write(x); putchar('\n');
}

template<typename T>void writesp(T x){
	write(x); putchar(' ');
}