LibreOJ 6686.Stupid GCD

真 TM 毒瘤……
典型的式子好推不可实现题……

一个显然的想法是先枚举 r=i3r = \lfloor\sqrt[3]i\rfloor,然后再考虑 ii 的贡献。
i=1i3gcd(i3,i)=r=1i3i=1n[i3=r]gcd(r,i)\sum\limits_{i=1}^{\lfloor\sqrt[3]i\rfloor}\gcd(\lfloor\sqrt[3]i\rfloor,i) = \sum\limits_{r = 1}^{\lfloor\sqrt[3]i\rfloor}\sum\limits_{i=1}^n[\lfloor\sqrt[3]i\rfloor=r]\gcd(r,i)

又因为 i3=r\lfloor\sqrt[3]i\rfloor = r 等价于 r3i<(r+1)3r^3 \le i < (r+1)^3,相当于把这些数分成了 i3\lfloor\sqrt[3]i\rfloor 块。
i=1i3gcd(i3,i)=i=n33ngcd(i3,i)+r=1i31i=r3(r+1)31gcd(r,i)\sum\limits_{i=1}^{\lfloor\sqrt[3]i\rfloor}\gcd(\lfloor\sqrt[3]i\rfloor,i) = \sum\limits_{i={\lfloor\sqrt[3]n\rfloor}^3}^n\gcd(\lfloor\sqrt[3]i\rfloor,i) + \sum\limits_{r=1}^{\lfloor\sqrt[3]i\rfloor-1}\sum\limits_{i=r^3}^{(r+1)^3-1}\gcd(r,i)

根据基本反演套路易得 i=1ngcd(r,i)=drndφ(d)\sum\limits_{i=1}^n\gcd(r,i) = \sum\limits_{d|r}\lfloor\frac nd\rfloor\varphi(d)
于是两个和式分别算即可。

第一个和式需要 n3\lfloor\sqrt[3]n\rfloor 的约数的欧拉函数值,可以考虑类似杜教筛的思想,预处理 n2/9n^{2/9} 之内的值,然后别的值分解质因数硬算(只需要 n1/6n^{1/6} 以下的质数)。
复杂度比较难算,姑且认为是 O(n2/9)O(n^{2/9})(实际上跑不满)。
(实际上可以通过 DFS 做到 O(n1/6)O(n^{1/6}),因为都是 n3\lfloor\sqrt[3]n\rfloor 的约数,但是我懒得写)

第二个和式,根据上式得

注意到 (dr+1)31d=d3r3+3d2r2+3drd=d2r3+3dr2+3r\left\lfloor\frac{(dr+1)^3-1}{d}\right\rfloor = \left\lfloor\frac{d^3r^3 + 3d^2r^2 + 3dr}{d}\right\rfloor = d^2r^3 + 3dr^2 + 3r,且 (dr)31d=d3r31d=d2r31\left\lfloor\frac{(dr)^3-1}{d}\right\rfloor = \left\lfloor\frac{d^3r^3-1}{d}\right\rfloor = d^2r^3 - 1

然后稍微杜教筛一下就行了。

代码:

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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#include <cstdio>
#include <cstring>
using namespace std;
template<class T>
inline void read(T &x)
{
char ch = 0,w = 0;
for(;ch < '0' || ch > '9';w |= ch == '-',ch = getchar());
for(x = 0;ch >= '0' && ch <= '9';x = (x << 3) + (x << 1) + (ch ^ '0'),ch = getchar());
w && (x = -x);
}
const int MX = 5e6;
const int mod = 998244353;
const int inv2 = 499122177;
const int inv6 = 166374059;
__int128 n;
long long rt;
int vis[MX + 5],cnt,prime[MX + 5],phi[MX + 5],sum[MX + 5];
inline long long cubert(__int128 n)
{
__int128 l = 1,r = n,mid,ret;
while(l <= r)
mid = l + r >> 1,mid <= n / mid / mid ? (l = mid + 1,ret = mid) : (r = mid - 1);
return ret;
}
int w[3][MX + 5];
inline int &mem_calc(long long x)
{
return w[0][rt / x];
}
inline int &mem_sphi(long long x)
{
return w[1][rt / x];
}
inline int &mem_sidphi(long long x)
{
return w[2][rt / x];
}
int calc(long long n)
{
if(n <= MX)
return phi[n] - phi[n - 1];
if(~mem_calc(n))
return mem_calc(n);
long long ret = n,t = n;
for(register int i = 1;(long long)prime[i] * prime[i] <= n;++i)
{
!(t % prime[i]) && (ret -= ret / prime[i]);
for(;!(t % prime[i]);t /= prime[i]);
}
t > 1 && (ret -= ret / t);
return mem_calc(n) = ret % mod;
}
int sphi(long long n)
{
if(n <= MX)
return phi[n];
if(~mem_sphi(n))
return mem_sphi(n);
long long ret = (n % mod) * ((n + 1) % mod) % mod * inv2 % mod;
for(register long long l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret = (ret - (r - l + 1) % mod * sphi(n / l) % mod + mod) % mod;
}
return mem_sphi(n) = ret;
}
int sidphi(long long n)
{
if(n <= MX)
return sum[n];
if(~mem_sidphi(n))
return mem_sidphi(n);
long long ret = (n % mod) * ((n + 1) % mod) % mod * ((2 * n + 1) % mod) % mod * inv6 % mod;
for(register long long l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret = (ret - ((l + r) % mod) * ((r - l + 1) % mod) % mod * inv2 % mod * sidphi(n / l) % mod + mod) % mod;
}
return mem_sidphi(n) = ret;
}
int ans;
int main()
{
memset(w,-1,sizeof w),sum[1] = phi[1] = 1;
for(register int i = 2;i <= MX;++i)
{
if(!vis[i])
phi[prime[++cnt] = i] = i - 1;
for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
{
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
sum[i] = (sum[i - 1] + (long long)i * phi[i] % mod) % mod,phi[i] = (phi[i] + phi[i - 1]) % mod;
}
read(n),rt = cubert(n);
for(register int i = 1;(long long)i * i <= rt;++i)
if(!(rt % i))
{
ans = (ans + (long long)((int)(n / i % mod) - (int)(((__int128)rt * rt * rt - 1) / i % mod) + mod) * calc(i) % mod) % mod;
if(i ^ rt / i)
ans = (ans + (long long)((int)(n / (rt / i) % mod) - (int)(((__int128)rt * rt * rt - 1) / (rt / i) % mod) + mod) * calc(rt / i) % mod) % mod;
}
--rt;
for(register long long l = 1,r;l <= rt;l = r + 1)
{
r = rt / (rt / l);
ans = (ans + 3LL * (sidphi(r) - sidphi(l - 1) + mod) % mod * (rt / l % mod) % mod * ((rt / l + 1) % mod) % mod * ((rt / l * 2 + 1) % mod) % mod * inv6 % mod) % mod;
ans = (ans + (long long)(sphi(r) - sphi(l - 1) + mod) * (3LL * (rt / l % mod) % mod * ((rt / l + 1) % mod) % mod * inv2 % mod + (rt / l) % mod) % mod) % mod;
}
printf("%d\n",ans);
}