\[ans=\sum_{i=1}^{a}\sum_{j=1}^b[gcd(a,b)=d] \]

我們假設

\[a\leq b \]

\[\begin{aligned} ans & = \sum_{i=1}^{a}\sum_{j=1}^b[gcd(a,b)=d] \\ & = \sum_{i=1}^{\left\lfloor\dfrac{a}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\dfrac{b}{d}\right\rfloor}[gcd(i,j)=1] \\ & = \sum_{i=1}^{\left\lfloor\dfrac{a}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\dfrac{b}{d}\right\rfloor}\sum_{t|gcd(i,j)}\mu(t) \\ & = \sum_{t=1}^{\left\lfloor\dfrac{a}{d}\right\rfloor}\mu(t)\left\lfloor\dfrac{a}{dt}\right\rfloor\left\lfloor\dfrac{b}{dt}\right\rfloor \end{aligned} \]

數論分塊

#include<cstdio>
#define Starseven main
#define ll long long
namespace lyt {
	void read(int &x){
	char ch=getchar();int re=0,op=1;
	while(ch<'0'||ch>'9'){if(ch=='-') op=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){re=(re<<3)+(re<<1)+ch-'0';ch=getchar();}
	x = re * op;
	return ;
	}
	void read(long long &x){
	char ch=getchar();long long re=0,op=1;
	while(ch<'0'||ch>'9'){if(ch=='-') op=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){re=(re<<3ll)+(re<<1ll)+ch-'0';ch=getchar();}
	x = re * op;
	return ;
	}
	void write(int x){
		if(x<0){putchar('-');x=-x;}
		if(x>9) write(x/10);
		putchar(x%10+'0');
		return ;
	}//記得自己加空格和換行 
	void write(long long x){
		if(x<0){putchar('-');x=-x;}
		if(x>9) write(x/10);
		putchar(x%10+'0');
		return ;
	}//記得自己加空格和換行
	int max(int x,int y){return x<y?y:x;}
	int min(int x,int y){return x<y?x:y;}
	int abs(int x){return x<0?-x:x;}
	long long max(long long x,long long y){return x<y?y:x;}
	long long min(long long x,long long y){return x<y?x:y;}
	long long abs(long long x){return x<0?-x:x;}
	double abs(double x){return x<0?-x:x;}
	void swap(int &a,int &b) {a ^= b ^= a ^= b;}
	void swap(long long &a,long long &b) {a ^= b ^= a ^= b;}
}using namespace lyt;
const int maxn = 1e6;
int a, b, d, prime[maxn + 20], num;
ll mu[maxn + 20];
bool vis[maxn + 20];

void Init() {
	mu[1] = 1;
	for (int i = 2; i <= maxn; i++) {
		if (!vis[i]) {
			prime[++num] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= num && prime[j] * i <= maxn; j++) {
			int x = prime[j] * i;
			vis[x] = true;
			if (i % prime[j] == 0) {
				mu[x] = 0;
				break;
			} 
			mu[x] = -mu[i];
		}
	}
	for (int i = 1; i <= maxn; i++) {
		mu[i] += mu[i - 1];
	}
	return ;
}

int Starseven(void) {
	read(a);
	read(b);
	read(d);
	Init();
	if (a > b) swap(a,b);
	a /= d; b /= d;
	ll ans = 0;
	for (int i = 1, j; i <= a; i = j + 1) {
		j = min(a / (a / i) , b / (b / i));
		ll ta = (ll)a / i, tb = (ll)b / i;
		ans += (mu[j] - mu[i - 1]) * ta * tb;
	}
	write(ans);
	puts("");
	return 0;	
}