简单反演+神仙优化题。
首先我们知道
σ 0 ( x y ) = ∑ i ∣ x ∑ j ∣ y [ ( i , j ) = 1 ] \sigma_0(xy)=\sum\limits_{i|x}\sum\limits_{j|y}[(i,j)=1] σ0(xy)=i∣x∑j∣y∑[(i,j)=1]
σ 0 ( x y z ) = ∑ i ∣ x ∑ j ∣ y ∑ k ∣ z [ ( i , j ) = 1 ] [ ( i , k ) = 1 ] [ ( j , k ) = 1 ] \sigma_0(xyz)=\sum\limits_{i|x}\sum\limits_{j|y}\sum\limits_{k|z}[(i,j)=1][(i,k)=1][(j,k)=1] σ0(xyz)=i∣x∑j∣y∑k∣z∑[(i,j)=1][(i,k)=1][(j,k)=1]
( σ 0 \sigma_0 σ0是除数函数,相当于题目描述中的 d d d)
只要将 x , y , z x,y,z x,y,z素因子分解后观察即可得到上面的结论。它的好处是出现了我们喜欢的 [ a = 1 ] [a=1] [a=1]
就可以莫比乌斯反演了
∑ i = 1 a ∑ j = 1 b ∑ k = 1 c σ 0 ( i j k ) \sum\limits_{i=1}^a\sum\limits_{j=1}^b\sum\limits_{k=1}^c\sigma_0(ijk) i=1∑aj=1∑bk=1∑cσ0(ijk)
= ∑ i = 1 a ∑ j = 1 b ∑ k = 1 c ∑ u ∣ i ∑ v ∣ j ∑ w ∣ k [ ( u , v ) = 1 ] [ ( u , w ) = 1 ] [ ( v , w ) = 1 ] =\sum\limits_{i=1}^a\sum\limits_{j=1}^b\sum\limits_{k=1}^c\sum\limits_{u|i}\sum\limits_{v|j}\sum\limits_{w|k}[(u,v)=1][(u,w)=1][(v,w)=1] =i=1∑aj=1∑bk=1∑cu∣i∑v∣j∑w∣k∑[(u,v)=1][(u,w)=1][(v,w)=1]
= ∑ u = 1 a ∑ v = 1 b ∑ w = 1 c ⌊ a u ⌋ ⌊ b v ⌋ ⌊ c w ⌋ [ ( u , v ) = 1 ] [ ( u , w ) = 1 ] [ ( v , w ) = 1 ] =\sum\limits_{u=1}^a\sum\limits_{v=1}^b\sum\limits_{w=1}^c\lfloor\frac{a}{u}\rfloor\lfloor\frac{b}{v}\rfloor\lfloor\frac{c}{w}\rfloor[(u,v)=1][(u,w)=1][(v,w)=1] =u=1∑av=1∑bw=1∑c⌊ua⌋⌊vb⌋⌊wc⌋[(u,v)=1][(u,w)=1][(v,w)=1]
= ∑ i = 1 a ∑ j = 1 b ∑ k = 1 c ⌊ a i ⌋ ⌊ b j ⌋ ⌊ c k ⌋ ∑ u ∣ ( i , j ) μ ( u ) ∑ v ∣ ( i , k ) μ ( v ) ∑ w ∣ ( j , k ) μ ( w ) =\sum\limits_{i=1}^a\sum\limits_{j=1}^b\sum\limits_{k=1}^c\lfloor\frac{a}{i}\rfloor\lfloor\frac{b}{j}\rfloor\lfloor\frac{c}{k}\rfloor\sum\limits_{u|(i,j)}\mu(u)\sum\limits_{v|(i,k)}\mu(v)\sum\limits_{w|(j,k)}\mu(w) =i=1∑aj=1∑bk=1∑c⌊ia⌋⌊jb⌋⌊kc⌋u∣(i,j)∑μ(u)v∣(i,k)∑μ(v)w∣(j,k)∑μ(w)
= ∑ u = 1 min ( a , b ) μ ( u ) ∑ v = 1 min ( a , c ) μ ( v ) ∑ w = 1 min ( b , c ) μ ( w ) ∑ [ u , v ] ∣ i a ⌊ a i ⌋ ∑ [ u , w ] ∣ j b ⌊ b j ⌋ ∑ [ v , w ] ∣ k c ⌊ c k ⌋ =\sum\limits_{u=1}^{\min(a,b)}\mu(u)\sum\limits_{v=1}^{\min(a,c)}\mu(v)\sum\limits_{w=1}^{\min(b,c)}\mu(w)\sum\limits_{[u,v]|i}^a\lfloor\frac{a}{i}\rfloor\sum\limits_{[u,w]|j}^b\lfloor\frac{b}{j}\rfloor\sum\limits_{[v,w]|k}^c\lfloor\frac{c}{k}\rfloor =u=1∑min(a,b)μ(u)v=1∑min(a,c)μ(v)w=1∑min(b,c)μ(w)[u,v]∣i∑a⌊ia⌋[u,w]∣j∑b⌊jb⌋[v,w]∣k∑c⌊kc⌋
设 f y ( x ) = ∑ x ∣ i y ⌊ y i ⌋ f_y(x)=\sum\limits_{x|i}^y\lfloor\frac{y}{i}\rfloor fy(x)=x∣i∑y⌊iy⌋,可以 O ( n log n ) O(n\log n) O(nlogn)预处理出 f a , f b , f c f_a,f_b,f_c fa,fb,fc所有值
a n s = ∑ u = 1 min ( a , b ) μ ( u ) ∑ v = 1 min ( a , c ) μ ( v ) ∑ w = 1 min ( b , c ) μ ( w ) f a ( [ u , v ] ) f b ( [ u , w ] ) f c ( [ v , w ] ) ans=\sum\limits_{u=1}^{\min(a,b)}\mu(u)\sum\limits_{v=1}^{\min(a,c)}\mu(v)\sum\limits_{w=1}^{\min(b,c)}\mu(w)f_a([u,v])f_b([u,w])f_c([v,w]) ans=u=1∑min(a,b)μ(u)v=1∑min(a,c)μ(v)w=1∑min(b,c)μ(w)fa([u,v])fb([u,w])fc([v,w])
然而痛苦的地方就在于这样推并不能降复杂度,目前为止我们的复杂度依然是 O ( n 3 ) O(n^3) O(n3)。事实上可以把它推到 O ( n 2 log n ) O(n^2\log n) O(n2logn)的,那样做可以过cf原题但是并不利于我们这题的后续优化。
但是也许我们可以不 O ( n 3 ) O(n^3) O(n3)地枚举 ( u , v , w ) (u,v,w) (u,v,w)三元组,因为事实上这样枚举的过程中有很多情况的贡献是 0 0 0:
-
如果 μ ( u ) , μ ( v ) , μ ( w ) \mu(u),\mu(v),\mu(w) μ(u),μ(v),μ(w)中有一个是 0 0 0,那么贡献是 0 0 0。这样的情况还挺多的。
-
如果 x > y x>y x>y,那么 f y ( x ) = 0 f_y(x)=0 fy(x)=0。也就是说 [ u , v ] > a [u,v]>a [u,v]>a或者 [ u , w ] > b [u,w]>b [u,w]>b或者 [ v , w ] > c [v,w]>c [v,w]>c都会导致贡献是 0 0 0。
所以我们有一种做法就是枚举最大公约数 g g g,枚举 u = i g , v = j g u=ig,v=jg u=ig,v=jg并注意满足条件 ( i , j ) = 1 , μ ( u ) ≠ 0 , μ ( v ) ≠ 0 , [ u , v ] = i j g ≤ max ( a , b , c ) (i,j)=1,\mu(u)\neq0,\mu(v)\neq0,[u,v]=ijg\leq\max(a,b,c) (i,j)=1,μ(u)̸=0,μ(v)̸=0,[u,v]=ijg≤max(a,b,c),连边 ( u , v ) (u,v) (u,v),建出一张图。这张图的三元环刚好对应到 ( u , v , w ) (u,v,w) (u,v,w)三元组,三元环计数的同时统计答案即可。
至于这样做的复杂度为什么是对的,我们已经发现有相当多的点对 ( u , v ) (u,v) (u,v)并不满足上述条件(意味着它们产生的贡献是 0 0 0),除去这些之后剩下的点对连出的边不会太多。正如shadowice大佬已经说了他实测下来这个边数只有七十多万。那么只要你的三元环计数速度合格就好。
关于三元环计数,可以参考iki9的这篇博客,我们可以在 O ( m m ) O(m\sqrt m) O(mm)的时间内完成三元环计数。
但是注意到上述连边没有考虑到自环。也就是说如果三元组 ( u , v , w ) (u,v,w) (u,v,w)中存在两个数相等或者三个数都相等,我们并没有将答案计算在内。因此我们需要单独算。这个很方便:三个数都相等的情况直接枚举这个数是几,把答案加起来就好;恰好有两个数相等的情况就在我们枚举 ( u , v ) (u,v) (u,v)连边的过程中解决,只要顺手把 ( u , u , v ) (u,u,v) (u,u,v)和 ( u , v , v ) (u,v,v) (u,v,v)的答案都加起来即可。
在做的过程中时刻注意剪枝,只要出现 μ \mu μ值为 0 0 0就跳过。另外可以顺手加一些循环展开、register之类的卡卡常。还有大家应该都注意到的cache miss问题,我们选择使用vector而不是邻接表存图以访问连续的内存来加速。
还有一个小trick:我们统计的是 i j k ijk ijk的约数个数之和,一个数的约数个数不会很多, i , j , k i,j,k i,j,k又都小于等于 1 0 5 10^5 105,估计最终的答案不会太大。事实上答案是不会超过 long long \text{long long} long long的,因此我们在统计的过程中完全不需要取模,只要在最后输出的时候模一下即可。
给一份稍微可读一点的代码
#include <cctype>
#include <cstdio>
#include <climits>
#include <algorithm>
#include <vector>
template <typename T> inline void read(T& x) {
int f = 0, c = getchar(); x = 0;
while (!isdigit(c)) f |= c == '-', c = getchar();
while (isdigit(c)) x = x * 10 + c - 48, c = getchar();
if (f) x = -x;
}
template <typename T> void write(T x) {
if (x < 0) putchar('-'), x = -x;
if (x > 9) write(x / 10);
putchar(x % 10 + 48);
}
template <typename T> inline void writeln(T x) { write(x); puts(""); }
template <typename T> inline bool chkmin(T &x, const T &y) { return y < x ? (x = y, 1) : 0; }
template <typename T> inline bool chkmax(T &x, const T &y) { return x > y ? (x = y, 1) : 0; }
#if __cplusplus >= 201103L
#define reg
#else
#define reg register
#endif
typedef long long LL;
const int maxn = 1e5 + 207;
const LL mod = 1e9 + 7;
int mu[maxn], pri[maxn];
bool ff[maxn];
LL fa[maxn], fb[maxn], fc[maxn];
int a, b, c, mx, mn;
struct Edge {
int u, v, w;
Edge(int x, int y, int z = 0) : u(x), v(y), w(z) {}
Edge() : u(0), v(0), w(0) {}
};
struct Node {
int to, w;
Node(int x, int y = 0) : to(x), w(y) {}
Node() : to(0), w(0) {}
};
std::vector<Node> G[maxn];
int deg[maxn];
Edge es[maxn << 4];
int vis[maxn];
int etot;
int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }
inline void sieve(int n) {
mu[1] = 1;
for (reg int i = 2; i <= n; ++i) {
if (!ff[i]) { pri[++pri[0]] = i; mu[i] = -1; }
for (reg int j = 1; j <= pri[0]; ++j) {
int x = i * pri[j];
if (x > n) break;
ff[x] = 1;
if (i % pri[j]) mu[x] = -mu[i];
else { mu[x] = 0; break; }
}
}
}
int main() {
int T; read(T);
sieve(1e5);
while (T--) {
read(a); read(b); read(c);
mx = std::max(std::max(a, b), c);
mn = std::min(std::min(a, b), c);
// 预处理 f 函数值
for (reg int i = 1; i <= mx; i += 4) {
fa[i] = 0; fa[i + 1] = 0; fa[i + 2] = 0; fa[i + 3] = 0;
fb[i] = 0; fb[i + 1] = 0; fb[i + 2] = 0; fb[i + 3] = 0;
fc[i] = 0; fc[i + 1] = 0; fc[i + 2] = 0; fc[i + 3] = 0;
}
for (reg int i = 1; i <= a; ++i)
for (reg int j = i; j <= a; j += i)
fa[i] += a / j;
for (reg int i = 1; i <= b; ++i)
for (reg int j = i; j <= b; j += i)
fb[i] += b / j;
for (reg int i = 1; i <= c; ++i)
for (reg int j = i; j <= c; j += i)
fc[i] += c / j;
// u = v = w
LL ans = 0;
for (reg int i = 1; i <= mn; ++i)
if (mu[i])
ans += mu[i] * mu[i] * mu[i] * fa[i] * fb[i] * fc[i];
// u, v, w 有两个相等
etot = 0;
for (reg int i = 1; i <= mx; i += 4) {
deg[i] = 0; deg[i + 1] = 0; deg[i + 2] = 0; deg[i + 3] = 0;
}
for (reg int g = 1; g <= mx; ++g)
for (reg int i = 1; i * g <= mx; ++i)
if (mu[i * g])
for (reg int j = i + 1; 1ll * i * j * g <= mx; ++j)
if (mu[j * g] && gcd(i, j) == 1) {
int u = i * g, v = j * g, lcm = i * j * g;
int uuv = mu[u] * mu[u] * mu[v], uvv = mu[u] * mu[v] * mu[v];
ans += uuv * (fa[u] * fb[lcm] * fc[lcm] + fa[lcm] * fb[u] * fc[lcm] + fa[lcm] * fb[lcm] * fc[u]);
ans += uvv * (fa[v] * fb[lcm] * fc[lcm] + fa[lcm] * fb[v] * fc[lcm] + fa[lcm] * fb[lcm] * fc[v]);
++deg[u]; ++deg[v];
es[++etot] = Edge(u, v, lcm);
}
// u, v, w 两两不同
for (reg int i = 1; i <= mx; ++i) G[i].clear();
for (reg int i = 1; i <= etot; ++i) {
int u = es[i].u, v = es[i].v, w = es[i].w;
if (deg[u] > deg[v] || (deg[u] == deg[v] && u < v))
G[u].push_back(Node(v, w));
else G[v].push_back(Node(u, w));
}
for (reg int x = 1; x <= mx; ++x) if (mu[x]) {
unsigned sz = G[x].size();
for (reg unsigned i = 0; i < sz; ++i)
vis[G[x][i].to] = G[x][i].w;
for (reg unsigned i = 0; i < sz; ++i) {
int y = G[x][i].to, wxy = G[x][i].w;
if (!mu[y]) continue;
for (reg unsigned j = 0, sszz = G[y].size(); j < sszz; ++j) {
int z = G[y][j].to;
if (vis[z]) {
int wyz = G[y][j].w, wxz = vis[z], mmuu = mu[x] * mu[y] * mu[z];
if (!mmuu) continue;
ans += mmuu * (fa[wxy] * fb[wyz] * fc[wxz]
+ fa[wxy] * fb[wxz] * fc[wyz]
+ fa[wyz] * fb[wxy] * fc[wxz]
+ fa[wyz] * fb[wxz] * fc[wxy]
+ fa[wxz] * fb[wxy] * fc[wyz]
+ fa[wxz] * fb[wyz] * fc[wxy]);
}
}
}
for (reg unsigned i = 0; i < sz; ++i)
vis[G[x][i].to] = 0;
}
writeln(ans % mod);
}
return 0;
}
另附 O ( n 2 log n ) O(n^2\log n) O(n2logn)的做法:
a n s = ∑ u = 1 a ∑ v = 1 b ∑ w = 1 c ⌊ a u ⌋ ⌊ b v ⌋ ⌊ c w ⌋ [ ( u , v ) = 1 ] [ ( u , w ) = 1 ] [ ( v , w ) = 1 ] ans=\sum\limits_{u=1}^a\sum\limits_{v=1}^b\sum\limits_{w=1}^c\lfloor\frac{a}{u}\rfloor\lfloor\frac{b}{v}\rfloor\lfloor\frac{c}{w}\rfloor[(u,v)=1][(u,w)=1][(v,w)=1] ans=u=1∑av=1∑bw=1∑c⌊ua⌋⌊vb⌋⌊wc⌋[(u,v)=1][(u,w)=1][(v,w)=1]
= ∑ i = 1 a ∑ j = 1 b ∑ k = 1 c ⌊ a i ⌋ ⌊ b j ⌋ ⌊ c k ⌋ ∑ d ∣ ( i , j ) μ ( d ) [ ( i , k ) = 1 ] [ ( j , k ) = 1 ] =\sum\limits_{i=1}^a\sum\limits_{j=1}^b\sum\limits_{k=1}^c\lfloor\frac{a}{i}\rfloor\lfloor\frac{b}{j}\rfloor\lfloor\frac{c}{k}\rfloor\sum\limits_{d|(i,j)}\mu(d)[(i,k)=1][(j,k)=1] =i=1∑aj=1∑bk=1∑c⌊ia⌋⌊jb⌋⌊kc⌋d∣(i,j)∑μ(d)[(i,k)=1][(j,k)=1]
= ∑ d = 1 min ( a , b ) μ ( d ) ∑ d ∣ i a ∑ d ∣ j b ⌊ a i ⌋ ⌊ b j ⌋ ∑ k = 1 c ⌊ c k ⌋ [ ( i , k ) = 1 ] [ ( j , k ) = 1 ] =\sum\limits_{d=1}^{\min(a,b)}\mu(d)\sum\limits_{d|i}^a\sum\limits_{d|j}^b\lfloor\frac{a}{i}\rfloor\lfloor\frac{b}{j}\rfloor\sum\limits_{k=1}^c\lfloor\frac{c}{k}\rfloor[(i,k)=1][(j,k)=1] =d=1∑min(a,b)μ(d)d∣i∑ad∣j∑b⌊ia⌋⌊jb⌋k=1∑c⌊kc⌋[(i,k)=1][(j,k)=1]
= ∑ d = 1 min ( a , b ) μ ( d ) ∑ i = 1 ⌊ a d ⌋ ∑ j = 1 ⌊ b d ⌋ ⌊ a i d ⌋ ⌊ b j d ⌋ ∑ k = 1 c ⌊ c k ⌋ [ ( i d , k ) = 1 ] [ ( j d , k ) = 1 ] =\sum\limits_{d=1}^{\min(a,b)}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{d}\rfloor}\lfloor\frac{a}{id}\rfloor\lfloor\frac{b}{jd}\rfloor\sum\limits_{k=1}^c\lfloor\frac{c}{k}\rfloor[(id,k)=1][(jd,k)=1] =d=1∑min(a,b)μ(d)i=1∑⌊da⌋j=1∑⌊db⌋⌊ida⌋⌊jdb⌋k=1∑c⌊kc⌋[(id,k)=1][(jd,k)=1]
= ∑ k = 1 c ⌊ c k ⌋ ∑ d = 1 min ( a , b ) μ ( d ) [ ( d , k ) = 1 ] ∑ i = 1 ⌊ a d ⌋ ⌊ a i d ⌋ [ ( i , k ) = 1 ] ∑ j = 1 ⌊ b d ⌋ ⌊ b j d ⌋ [ ( j , k ) = 1 ] =\sum\limits_{k=1}^c\lfloor\frac{c}{k}\rfloor\sum\limits_{d=1}^{\min(a,b)}\mu(d)[(d,k)=1]\sum\limits_{i=1}^{\lfloor\frac{a}{d}\rfloor}\lfloor\frac{a}{id}\rfloor[(i,k)=1]\sum\limits_{j=1}^{\lfloor\frac{b}{d}\rfloor}\lfloor\frac{b}{jd}\rfloor[(j,k)=1] =k=1∑c⌊kc⌋d=1∑min(a,b)μ(d)[(d,k)=1]i=1∑⌊da⌋⌊ida⌋[(i,k)=1]j=1∑⌊db⌋⌊jdb⌋[(j,k)=1]
设函数 f ( x , y ) = ∑ i = 1 x ⌊ x i ⌋ [ ( i , y ) = 1 ] f(x,y)=\sum\limits_{i=1}^x\lfloor\frac{x}{i}\rfloor[(i,y)=1] f(x,y)=i=1∑x⌊ix⌋[(i,y)=1]那么
a n s = ∑ k = 1 c ⌊ c k ⌋ ∑ d = 1 min ( a , b ) μ ( d ) [ ( d , k ) = 1 ] f ( ⌊ a d ⌋ , k ) f ( ⌊ b d ⌋ , k ) ans=\sum\limits_{k=1}^c\lfloor\frac{c}{k}\rfloor\sum\limits_{d=1}^{\min(a,b)}\mu(d)[(d,k)=1]f(\lfloor\frac{a}{d}\rfloor,k)f(\lfloor\frac{b}{d}\rfloor,k) ans=k=1∑c⌊kc⌋d=1∑min(a,b)μ(d)[(d,k)=1]f(⌊da⌋,k)f(⌊db⌋,k)
枚举 k k k,枚举 d d d,这个 f f f直接 O ( x ) O(x) O(x)暴力就行了,因为 x = a / d x=a/d x=a/d所以枚举 d d d和计算 f f f的总复杂度是 O ( n log n ) O(n\log n) O(nlogn)(调和级数乘以 n n n),总复杂度 O ( n 2 log n ) O(n^2\log n) O(n2logn)