题目

题目描述
对于所有 x ∈ [ A , B ] ,    y ∈ [ C , D ] x\in[A,B],\;y\in[C,D] x∈[A,B],y∈[C,D] 的二元组 x , y x,y x,y，记 d = gcd ⁡ ( x , y ) d=\gcd(x,y) d=gcd(x,y)，如果 x + y d ⩽ V = 999 {x+y\over d}\leqslant V=999 dx+y​⩽V=999，那么将计数器加上 x + y d {x+y\over d} dx+y​ 。

请求出最后计数器的值。

数据范围与提示
A , B , C , D ⩽ 1 0 12 A,B,C,D\leqslant 10^{12} A,B,C,D⩽1012 。但是事实证明这是误导信息。

思路

先讲我的思路。因为大家一看正解，就觉得我的做法是*。

为了方便，容斥一下，只考虑 x ⩽ A x\leqslant A x⩽A 和 y ⩽ B y\leqslant B y⩽B 的约束。一看这个 1 0 12 10^{12} 1012，肯定联想 A \sqrt A A ​ 啊，整除分块？试一试吧。

考虑枚举 d d d，那么 x d , y d {x\over d},{y\over d} dx​,dy​ 有一个范围。但是还要保证二者互质，所以利用莫比乌斯函数。可以得到式子是
∑ d ∑ p μ ( p ) ∑ x ⩽ ⌊ A d p ⌋ ∑ y ⩽ ⌊ B d p ⌋ [ x + y ⩽ ⌊ V p ⌋ ] ( x p + y p ) \sum_{d}\sum_{p} \mu(p) \sum_{x\leqslant\lfloor{A\over dp}\rfloor} \sum_{y\leqslant\lfloor{B\over dp}\rfloor} \left[x+y\leqslant\left\lfloor{V\over p}\right\rfloor\right] (xp+yp) d∑​p∑​μ(p)x⩽⌊dpA​⌋∑​y⩽⌊dpB​⌋∑​[x+y⩽⌊pV​⌋](xp+yp)

那么显然，当 p p p 一定时，后面 x , y x,y x,y 的范围只由 d d d 决定，并且是整除分块的形式。枚举 p p p 则复杂度为 O ( V A ) \mathcal O(V\sqrt{A}) O(VA ​) 。

测试了一下，发现 A = 1 0 12 A=10^{12} A=1012 的运行速度极其慢。于是考虑剪枝。当 ⌊ A d p ⌋ ⩾ ⌊ V p ⌋ \lfloor{A\over dp}\rfloor\geqslant\lfloor{V\over p}\rfloor ⌊dpA​⌋⩾⌊pV​⌋ 时，这个范围就失去了意义，因为 x + y ⩽ ⌊ V p ⌋ x+y\leqslant\lfloor{V\over p}\rfloor x+y⩽⌊pV​⌋ 已经保证了这一点。所以就把这样的 d d d 拿出来统一计算。

那么，我们只需要计算 ⌊ A d p ⌋ < ⌊ V p ⌋ \lfloor{A\over dp}\rfloor<\lfloor{V\over p}\rfloor ⌊dpA​⌋<⌊pV​⌋ 的部分。由于我们是数论分块，本就会让 ⌊ A d p ⌋ \lfloor{A\over dp}\rfloor ⌊dpA​⌋ 的值每次不同，所以复杂度是 O ( ⌊ V p ⌋ ) \mathcal O(\lfloor{V\over p}\rfloor) O(⌊pV​⌋) 。注意我们要对 B B B 也进行同样的操作，才能保证复杂度。

诶，好像复杂度是 O ( V ln ⁡ V ) \mathcal O(V\ln V) O(VlnV) 的，搞定了？并且最终复杂度与 A A A 无关？

正解

枚举两个互质的数 x d , y d ∈ [ 0 , V ] {x\over d},{y\over d}\in[0,V] dx​,dy​∈[0,V]，根据 x , y x,y x,y 的范围可以求出 d d d 的范围，也就得到了这种方案的数量。时间复杂度 O ( V 2 ) \mathcal O(V^2) O(V2)，注意用二维数组存储 gcd ⁡ \gcd gcd，否则可能变为 O ( V 2 ln ⁡ V ) \mathcal O(V^2\ln V) O(V2lnV) 。

代码

代码实现中，我还对 V p , A p , B p \frac Vp,\frac Ap,\frac Bp pV​,pA​,pB​ 作了整除分块。虽然它实际上没太大用……

#include <cstdio>
#include <iostream>
#include <cstring>
#include <vector>
#include <algorithm>
using namespace std;
typedef long long int_;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}
inline void writeint(int x){
	if(x > 9) writeint(x/10);
	putchar((x-x/10*10)^48);
}

const int Mod = 1e9+7;
const int inv2 = (Mod+1)>>1;
const int V = 999;
int mu[V+1]; bool isPrime[V+1];
void sieve(){
	memset(isPrime+2,1,V-1);
	rep(i,1,V) mu[i] = 1;
	rep(i,2,V) if(isPrime[i]){
		mu[i] = -1;
		for(int j=2; j<=V/i; ++j){
			if(j%i == 0) mu[i*j] = 0;
			isPrime[i*j] = false;
			mu[i*j] = -mu[i*j];
		}
	}
}

int jb[V+1]; // a*(a+1)*(2*a+1)/6+a*(a+1)/2
int calc(int a,int b,const int &X){
	if(a+b < X) // never explode
		return ((b*a*(a+1)>>1)
			+(a*b*(b+1)>>1))%Mod;
	int res = (jb[a]+(X-a-1)*a*(a+1))%Mod;
	res += (jb[b]+(X-b-1)*b*(b+1))%Mod;
	res -= jb[X-a-1]; res -= jb[X-b-1];
	res = (res%Mod+Mod)%Mod;
	return int(int_(res)*inv2%Mod);
}
int wxk[V+1]; // prefix sum of mu[i]*i
int solve(int_ a0,int_ b0){
	if(a0 > b0) swap(a0,b0);
	int res = 0, end_p = min(a0,int_(V>>1));
	for(int p=1,rp; p<=end_p; p=rp+1){
		rp = min(a0/(a0/p),b0/(b0/p));
		const int X = V/p; // real boundary
		rp = min(rp,V/X); // keep the same
		int_ a = a0/p, b = b0/p;
		int_ l = a/(X-1)+1; // or X-1 <= a/l
		int_ ddg = b/(X-1); // or X-1 <= b/l
		int now = calc(X-1,X-1,X)*(l-1)%Mod;
		for(int_ r; l<=a; l=r){
			r = min(a/(a/l),b/(b/l))+1;
			if(b/l >= X-1) // useless
				r = min(ddg,a/(a/l))+1;
			int ta = min(a/l,int_(X-1));
			int tb = min(b/l,int_(X-1));
			now = (now+calc(ta,tb,X)*(r-l))%Mod;
		}
		int t = wxk[rp]-wxk[p-1];
		res = (res+int_(now)*t)%Mod;
	}
	return (res%Mod+Mod)%Mod;
}

int main(){
	rep(i,1,V) jb[i] = jb[i-1]+i*(i+1);
	sieve(); int_ a, b, c, d;
	rep(i,1,V) wxk[i] = wxk[i-1]+mu[i]*i;
	scanf("%lld %lld %lld %lld",&a,&b,&c,&d);
	int res = solve(b,d)+solve(a-1,c-1);
	res -= solve(a-1,d)+solve(b,c-1);
	printf("%d\n",(res%Mod+Mod)%Mod);
	return 0;
}