题目
题目描述
对于所有
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];
int calc(int a,int b,const int &X){
if(a+b < X)
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];
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;
rp = min(rp,V/X);
int_ a = a0/p, b = b0/p;
int_ l = a/(X-1)+1;
int_ ddg = b/(X-1);
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)
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;
}
|