「SDOI 2018」旧试题

Problem

Description

给定 A,B,C ,求

\sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} \tau\left(ijk\right)

答案对 10^{9}+7 取模。

每个测试点共 T 组数据。

Constraints

1\le T\le 10, 1\le A, B, C\le 10^{5}, 1\le \sum\max\left(A, B, C\right)\le 2\times 10^{5}

Solution

Analysis

先来推一波柿子 qwq

\begin{aligned} \DeclareMathOperator*{\lcm}{lcm} & \sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} \tau(ijk) \\ =& \sum_{i = 1}^{A} \sum_{j = 1}^{B} \sum_{k = 1}^{C} \sum_{u | i} \sum_{v | j} \sum_{w | k} [u \perp v] [v \perp w] [w \perp u] \\ =& \sum_{x = 1}^{A} \sum_{y = 1}^{B} \sum_{z = 1}^{C} [x \perp y] [y \perp z] [z \perp x] (A / x) (B / y) (C / z) \\ =& \sum_{x = 1}^{A} \sum_{y = 1}^{B} \sum_{z = 1}^{C} (A / x) (B / y) (C / z) \sum_{u | x, u | y} \mu(u) \sum_{v | y, v | z} \mu(v) \sum_{w | z, w | x} \mu(w) \\ =& \sum_{u = 1}^{A} \sum_{v = 1}^{B} \sum_{w = 1}^{C} \mu(u) \mu(v) \mu(w) \sum_{w | x, u | x} (A / x) \sum_{u | y,v | y} (B / y) \sum_{v | z,w | z} (C / z) \\ =& \sum_{u = 1}^{A} \sum_{v = 1}^{B} \sum_{w = 1}^{C} \mu(u) \mu(v) \mu(w) \sum_{\lcm(w, u) | x} (A / x) \sum_{\lcm(u, v) | y} (B / y) \sum_{\lcm(v, w) | z} (C / z) \end{aligned}

可知

\sum_{x | d} (n / d) = \sum_{ix \le n} (n / ix )= \sum_{ix \le n} ((n / x) / i) = \sum_{ix \le n} \tau(i)

\tau(n) 的前缀和为 f(n)

f(n)= \sum_{i = 1}^{n}\tau(i)

故答案即为

\sum_{u = 1}^{A} \sum_{v = 1}^{B} \sum_{w = 1}^{C} \mu(u) \mu(v) \mu(w) f(A / \lcm(w, u)) f(B / \lcm(u, v)) f(C / \lcm(v, w))

可以发现所有有贡献的三元组 (u, v, w) 必定满足:

\mu(u) \mu(v) \mu(w)\ne 0,\lcm(w, u)\le A,\lcm(u, v)\le B,\lcm(v, w)\le C

故对于每个质数 p ,枚举其被 u, v, w 整除的八种情况即可。

发现对于 u, v, w 中有至少两个能被 p 整除的四种情况,其贡献的绝对值是相同的。
因为当 u, v, w 中至少两个能被整除时, \lcm(w, u),\lcm(u, v),\lcm(v, w) 必定都能被 p 整除。
则四种情况贡献的绝对值均为:

f\left(A / \left(\lcm(w, u) / p\right)\right) f\left(B / \left(\lcm(u, v) / p\right)\right) f\left(C / \left(\lcm(v, w) / p\right)\right)

这四种情况贡献之和的系数即为 3 \mu^{2}(p)+ \mu^{3}(p) = 2
故只需枚举五种情况爆搜即可。

可以发现当运用此优化后,从大到小搜时状态会显著减少,故采用从大到小搜索的方式。
又发现若从大到小搜索,则当 \max(A, B, C)<p p 的取值并不会影响当前状态的答案。

则可以预处理一部分的 \tau(ijk) ,当 \max(A,B,C)<p \max(A,B,C) 处于预处理范围中时直接返回预处理的答案。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include<cstdio>
#include<algorithm>

using i64=long long;
constexpr int maxn(100000);
constexpr int U(100000);
constexpr int M(46);
constexpr int p(1000000007);

template<typename _Tp>
inline void inc(_Tp&x,const _Tp&y)
{x+=y;(p<=x)&&(x-=p);}
template<typename _Tp>
inline void dec(_Tp&x,const _Tp&y)
{x-=y;(x<0)&&(x+=p);}
template<typename _Tp>
inline _Tp mtp2(const _Tp&v)
{return(v>=p-v)?(v-p+v):(v+v);}

template<typename _Tp>
inline _Tp Max(const _Tp&x,const _Tp&y)
{return x>y?x:y;}
template<typename _Tp>
inline void chk(_Tp&x)
{(x<0)?x+=p:(p<=x)&&(x-=p);}

bool notp[U+1];
int pri[U/7];
int cnt;

int tau[maxn+1];
int stau[maxn+1];
int ans[M+1][M+1][M+1];

void Sieve(const int&n){
std::fill(tau+1,tau+n+1,1);

for(int i=2;i<=n;++i)if(!notp[i]){
pri[++cnt]=i;
for(int j=i,k=1;j<=n;j+=i,++k)
notp[j]=true,tau[j]+=tau[k];
}

std::reverse(pri+1,pri+cnt+1);

for(int i=1,v=0;i<=n;++i)
v+=tau[i],stau[i]=v;

for(int i=1;i<=M;++i)
for(int j=1;j<=M;++j)
for(int k=1;k<=M;++k)
ans[i][j][k]=(i64)(
ans[i-1][j][k] +ans[i][j-1][k] +ans[i][j][k-1]
-ans[i][j-1][k-1] -ans[i-1][j][k-1] -ans[i-1][j-1][k]
+ans[i-1][j-1][k-1]
+tau[i*j*k]
)%p;
}

int calc(const int&c,const int&A,const int&B,const int&C){
if(!pri[c])
return(i64)stau[A]*stau[B]%p*stau[C]%p;
int pn=pri[c];
int n=Max(Max(A,B),C);
if(pn>n&&n<=M)
return ans[A][B][C];
if(!(A&&B&&C))
return 0;

int v=calc(c+1,A,B,C);
if(pn<=A&&pn<=B)
dec(v,calc(c+1,A/pn,B/pn,C));
if(pn<=B&&pn<=C)
dec(v,calc(c+1,A,B/pn,C/pn));
if(pn<=C&&pn<=A)
dec(v,calc(c+1,A/pn,B,C/pn));
if(pn<=A&&pn<=B&&pn<=C){
int cv=calc(c+1,A/pn,B/pn,C/pn);
chk(cv<<=1);chk(v+=cv);
}

return v;
}

int main(){
Sieve(U);
int T;
for(scanf(" %d",&T);T;--T){
int A,B,C;
scanf(" %d %d %d",&A,&B,&C);
int n=Max(Max(A,B),C);
int c=1;
for(;n<pri[c];++c);
printf("%d\n",calc(c,A,B,C));
}return 0;
}