题意:给定 nnn,求 ∑i=1ngcd(⌊i3⌋,i) mod 998244353\displaystyle \sum_{i=1}^n \gcd\left( \left\lfloor\sqrt[3]{i}\right\rfloor,i\right) \bmod 998244353i=1∑ngcd(⌊3i⌋,i)mod998244353,其中 n≤1×1021n \leq 1\times 10^{21}n≤1×1021。
首先考虑分块,根据 i3\sqrt[3]{i}3i 的值可以将 [1,n][1,n][1,n] 的整数分成 n3\sqrt[3]{n}3n 块。记 m=⌊n3⌋m=\left \lfloor \sqrt[3]{n} \right \rfloorm=⌊3n⌋,原式可化为 (∑i=1m∑j=i3(i+1)3−1gcd(i,j))+∑i=m3ngcd(m,i)\displaystyle \left (\sum_{i=1}^m \sum_{j=i^3}^{(i+1)^3-1}\gcd(i,j) \right )+\sum_{i=m^3}^n \gcd(m,i)⎝⎛i=1∑mj=i3∑(i+1)3−1gcd(i,j)⎠⎞+i=m3∑ngcd(m,i)。首先来考虑固定 iii,整块的做法。
首先容易注意到,gcd(i,j)=gcd(j,i mod j)\gcd(i,j)=\gcd(j,i \bmod j)gcd(i,j)=gcd(j,imodj),因而对于一个固定的 ddd,∑i=d3d3+3d2+3dgcd(d,i)=∑i=d3d3+3d2+3dgcd(d,i mod d)\displaystyle \sum_{i=d^3}^{d^3+3d^2+3d} \gcd(d,i)=\sum_{i=d^3}^{d^3+3d^2+3d}\gcd(d,i \bmod d)i=d3∑d3+3d2+3dgcd(d,i)=i=d3∑d3+3d2+3dgcd(d,imodd)。而集合 [d3,d3+3d2+3d][d^3,d^3+3d^2+3d][d3,d3+3d2+3d] 可以拆分为 3d2+13d^2+13d2+1 个完整的 ∑i=1dgcd(i,d)\displaystyle \sum_{i=1}^d \gcd(i,d)i=1∑dgcd(i,d) 和最后的一项 gcd(d,d3+3d2+3d)=d\gcd(d,d^3+3d^2+3d)=dgcd(d,d3+3d2+3d)=d,因而原式前半部分可化为:
∑d=1m∑i=d3d3+3d2+3dgcd(i,d)=∑d=1m((3d+1)∑i=1dgcd(d,i))+d
\sum_{d=1}^{m}\sum_{i=d^3}^{d^3+3d^2+3d} \gcd(i,d)=\sum_{d=1}^{m} \left((3d+1)\sum_{i=1}^d \gcd(d,i)\right)+d
d=1∑mi=d3∑d3+3d2+3dgcd(i,d)=d=1∑m((3d+1)i=1∑dgcd(d,i))+d
对于最后一部分,式子也是形如 ∑i=m3ngcd(m,i)\displaystyle \sum_{i=m^3}^n \gcd(m,i)i=m3∑ngcd(m,i),考虑 ∑i=1ngcd(i,a)\displaystyle \sum_{i=1}^n \gcd(i,a)i=1∑ngcd(i,a) 的求法。
∑i=1ngcd(i,a)=∑d∣ad∑i=1n[gcd(i,a)=d]=∑d∣ad∑i=1⌊nd⌋∑ed∣(di),e∣adμ(e)=t=ed∑d∣ad∑i=1⌊nd⌋∑t∣di,t∣aμ(td)=∑t∣a⌊nt⌋∑d∣tdμ(td)=∑t∣a⌊nt⌋φ(t)
\begin{aligned}
\sum_{i=1}^n \gcd(i,a)=&\sum_{d|a} d\sum_{i=1}^n[\gcd(i,a)=d]\\
=&\sum_{d|a}d\sum_{i=1}^{\lfloor{\frac{n}{d}} \rfloor} \sum_{ed|(di),e|\frac{a}{d}} \mu(e)\\
\overset{\underset{t=ed}{}}=&\sum_{d|a}d\sum_{i=1}^{\lfloor{\frac{n}{d}} \rfloor} \sum_{t|di,t|a} \mu \left(\dfrac{t}{d}\right)\\
=&\sum_{t|a} \left\lfloor\frac{n}{t}\right \rfloor\sum_{d|t} d\mu\left(\dfrac{t}{d}\right)\\
=&\sum_{t|a} \left \lfloor\dfrac{n}{t} \right \rfloor\varphi(t)
\end{aligned}
i=1∑ngcd(i,a)===t=ed==d∣a∑di=1∑n[gcd(i,a)=d]d∣a∑di=1∑⌊dn⌋ed∣(di),e∣da∑μ(e)d∣a∑di=1∑⌊dn⌋t∣di,t∣a∑μ(dt)t∣a∑⌊tn⌋d∣t∑dμ(dt)t∣a∑⌊tn⌋φ(t)
若 a≠na \neq na=n,对于单个的 aaa 可以考虑在 O(a)\mathcal{O}(\sqrt{a})O(a) 的时间内枚举因子求得。对于大部分的情况——a=na=na=n,原式可化简为 ∑t∣nntφ(t)=φ∗id\displaystyle \sum_{t|n} \dfrac{n}{t} \varphi(t)=\varphi*idt∣n∑tnφ(t)=φ∗id,因而是一个积性函数,不妨设为 f(n)f(n)f(n),考虑欧拉筛线性的求得。
欧拉筛的难点在于 f(pk)→f(pk+1)f(p^k) \to f(p^{k+1})f(pk)→f(pk+1) 的转化。
f(pk)=∑t∣pkpktφ(t)=∑i=0kpk−iφ(pi)=pk+∑i=1kpk−ipi−1(p−1)=p+(p−1)∑i=1kpk−1=p+(p−1)kpk−1=p(f(pk−1)+φ(pk−1))
\begin{aligned}
f(p^k)=\sum_{t|p^k} \dfrac{p^k}{t} \varphi(t)=&\sum_{i=0}^k p^{k-i} \varphi(p^i)\\
=&p^k+\sum_{i=1}^k p^{k-i} p^{i-1}(p-1)\\
=&p+(p-1)\sum_{i=1}^k p^{k-1}\\
=&p+(p-1)kp^{k-1}\\
=&p(f(p^{k-1})+\varphi(p^{k-1}))
\end{aligned}
f(pk)=t∣pk∑tpkφ(t)=====i=0∑kpk−iφ(pi)pk+i=1∑kpk−ipi−1(p−1)p+(p−1)i=1∑kpk−1p+(p−1)kpk−1p(f(pk−1)+φ(pk−1))
因而维护一个 φ\varphiφ 数组,最小质因子 M\mathcal MM,在该数中最小质因子 T\mathcal TT 的次方数(例如,T12=22=4{\mathcal T}_{12}=2^2=4T12=22=4)即可。对于 p∤np \nmid np∤n,则 f(np)=f(n)f(p)=f(n)(2p−1)f(np)=f(n)f(p)=f(n)(2p-1)f(np)=f(n)f(p)=f(n)(2p−1);若 p∣np|np∣n,考虑 Tn=pi{\mathcal T}_n=p^iTn=pi。
f(np)=f(npi×pi+1)=f(npi)f(pi+1)=f(npi)×p(f(pi)+φ(pi)))=p[f(n)+f(npi)φ(pi)]
\begin{aligned}
f(np)=&f\left(\dfrac{n}{p^i} \times p^{i+1}\right)\\
=&f\left(\dfrac{n}{p^i}\right)f(p^{i+1})\\
=&f\left(\dfrac{n}{p^i}\right)\times p(f(p^i)+\varphi(p^i)))\\
=&p\left[f(n)+f\left(\dfrac{n}{p^i}\right) \varphi(p^i)\right]
\end{aligned}
f(np)====f(pin×pi+1)f(pin)f(pi+1)f(pin)×p(f(pi)+φ(pi)))p[f(n)+f(pin)φ(pi)]
因而同步的转移即可。
#include <bits/stdc++.h>
using namespace std;
const long long mod = 998244353;
const int N = 10000000;
char a[30];
int sum[N + 5];
int f[N + 5], prime[N + 5], phi[N + 5], minfac[N + 5], mintimes[N + 5];
int tot;
long long power(long long a, long long x)
{
long long ans = 1;
while(x)
{
if (x & 1)
ans = ans * a % mod;
a = a % a % mod;
x >>= 1;
}
return ans;
}
long long inv(long long a)
{
return power(a, mod - 2);
}
inline void sieve(int n)
{
phi[1] = f[1] = 1;
for (int i = 2; i <= n; i++)
{
if(!minfac[i])
{
minfac[i] = prime[++tot] = mintimes[i] = i;
phi[i] = i - 1;
f[i] = i + i - 1;
}
for (int j = 1; j <= tot && (long long)prime[j] * i <= n; j++)
{
long long num = prime[j] * i;
if (i % prime[j])
{
minfac[num] = mintimes[num] = prime[j];
phi[num] = (long long)phi[i] * (prime[j] - 1);
f[num] = (long long)f[i] * (2 * prime[j] - 1) % mod;
}
else
{
minfac[num] = prime[j];
mintimes[num] = (long long)prime[j] * mintimes[i];
phi[num] = (long long)phi[i] * prime[j];
f[num] = (f[i] + (long long)f[i / mintimes[i]] * phi[mintimes[i]] % mod) % mod * (long long)prime[j] % mod;
}
}
}
for (int i = 1; i <= n; i++)
sum[i] = ((3ll * i + 3) * f[i] % mod + i + sum[i - 1]) % mod;
}
long long g(long long n, long long m)
{
long long ans = 0;
for (int i = 1; i * i <= m;i++)
if (m % i == 0)
{
int j = m / i;
ans = (ans + phi[i] * (n / i) % mod) % mod;
if (i != j)
ans = (ans + phi[j] * (n / j) % mod) % mod;
}
return ans;
}
int main()
{
sieve(N);
int t;
scanf("%d", &t);
while(t--)
{
__int128_t n = 0;
scanf("%s", a);
int len = strlen(a);
for (int i = 0; i < len;i++)
n = n * 10 + a[i] - 48;
__int128_t left = 0, right = N, lim = 0;
while (left <= right)
{
__int128_t mid = (left + right) >> 1;
if (mid * mid * mid <= n)
{
lim = mid;
left = mid + 1;
}
else
right = mid - 1;
}
long long ans = (((n / lim - lim * lim) % mod * f[lim] % mod + lim) % mod + g(n % lim, lim)) % mod;
ans = (ans + sum[lim - 1]) % mod;
printf("%lld\n", ans);
}
return 0;
}