LOJ3069. 「2019 集训队互测 Day 1」整点计数(min_25筛)

时间:2023-11-11 13:09:32

题目链接

https://loj.ac/problem/3069

题解

复数真神奇。

一句话题意:令 \(f(x)\) 表示以原点 \((0, 0)\) 为圆心,半径为 \(x\) 的圆上的整点数量,求 \(\sum_\limits{i = 1}^N f(i)^k \bmod P\) 的值。

令 \(g(x) = \frac{f(x)}{4}\),那么我们需要求 \(\left(4^k\sum_\limits{i = 1}^Ng(i)^k\right) \bmod P\)。打表可得 \(g(x)\) 是一个积性函数,其证明如下:

首先显然有 \(g(1) = 1\)。

对于正整数 \(x \geq 2\),\(g(x)\) 有如下求法:

记 \(y = x^2\),假设 \(y\) 的唯一分解式中,形如 \(4n + 1(n \in \mathbb{N}^*)\) 的质数对应的指数分别为 \(\alpha_1, \alpha_2, \cdots, \alpha_m\),形如 \(4n + 3(n \in \mathbb{N})\) 的质数对应的指数分别为 \(\beta_1, \beta_2, \cdots, \beta_k\),那么有:

\[\begin{aligned}g(x) = \left(\prod_{i = 1}^m (\alpha_i + 1)\right) \left(\prod_{i = 1}^k [\beta_i 是偶数]\right)\end{aligned}
\]

关于上述 \(g(x)\) 的求法的具体原因及相关理论,可参考视频 隐藏在素数规律中的 π。包看包懂。

由于 \(y = x^2\),因此 \(y\) 的唯一分解式中任意质数的指数一定是偶数,即上述式子中后面部分 \(\left(\prod_\limits{i = 1}^k [\beta_i 是偶数]\right)\) 的值一定为 \(1\)。故有:

\[\begin{aligned}g(x) = \prod_{i = 1}^m (\alpha_i + 1)\end{aligned}
\]

对于互质的两个数 \(a, b\),显然 \(a^2, b^2\) 的唯一分解式中不会出现相同的质数,故此时有 \(g(ab) = g(a)g(b)\)。因此 \(g(x)\) 是一个积性函数。

\(g(x)\) 是积性函数那么 \(g(x)^k\) 显然也是积性函数。考虑使用 min_25 筛来求积性函数 \(g(x)^k\) 的前缀和。根据上面的理论,不难得出对于一个质数 \(p\),有

\[g(p^c)^k = \begin{cases}(2c + 1)^k, &p \bmod 4 = 1 \\ 1, &p \bmod 4 \neq 1 \end{cases}
\]

要注意的是,在算 \(g(p^c)\) 时,实际是对 \((p^c)^2\) 做唯一分解,因此当 \(p \bmod 4 = 1\) 时,\(g(p^c)\) 的值为 \(2c + 1\) 而不是 \(c + 1\)。

质数的任意次幂对应的函数值我们已经能够快速计算。接下来考虑如何求所有质数对应的函数值的前缀和。考虑先分别求出 \(1\sim n\) 内模 \(4\) 等于 \(1\) 的质数的数量以及模 \(4\) 等于 \(3\) 的质数的数量,不妨将两种质数对应的数量记为 \(S_1, S_3\),那么边界值(即将 \(2 \sim n\) 的所有数均视为质数时对应的值)显然是 \(S_1(n) = \left\lfloor\frac{n - 1}{4}\right\rfloor,\) \(S_3(n) = \left\lfloor\frac{n + 1}{4}\right\rfloor\)。考虑如何依次去掉以每一个 \(\geq 3\) 的质数作为最小质因子的所有数的贡献。注意到在模 \(4\) 意义下,有:当 \(a = 1\) 时,\(\begin{cases}a \times 1 = 1 \\ a \times 3 = 3\end{cases}\);当 \(a = 3\) 时,\(\begin{cases}a \times 1 = 3 \\ a \times 3 = 1\end{cases}\),因此对于每个 \(\geq 3\) 的质数 \(p\),\(S_1\) 与 \(S_3\) 的更新方式如下:

\[\begin{aligned}\text{if $p$ mod 4 = 1: }&\begin{cases} {S_1(x) \leftarrow S_1(x) - \left(S_1 (\left\lfloor\frac{x}{p}\right\rfloor) - S_1(p - 1)\right) \\S_3(x) \leftarrow S_3(x) - \left(S_3 (\left\lfloor\frac{x}{p}\right\rfloor) - S_3(p - 1)\right)} \end{cases} \\ \text{else: }&\begin{cases}{S_1(x) \leftarrow S_1(x) - \left(S_3 (\left\lfloor\frac{x}{p}\right\rfloor) - S_3(p - 1)\right) \\S_3(x) \leftarrow S_3(x) - \left(S_1 (\left\lfloor\frac{x}{p}\right\rfloor) - S_1(p - 1)\right)}\end{cases} \end{aligned}
\]

令 \(S(x)\) 为 \(1 \sim x\) 内的所有质数 \(p\) 对应的 \(g(p)^k\) 的值,那么有 \(S(x) = S_1(x) \times 3^k + S_3(x)\)。注意还应添加质数 \(2\) 的贡献。

时间复杂度为 \(O(\frac{n^{\frac{3}{4}}}{\log n})\)。

代码

#include<bits/stdc++.h>

using namespace std;

const int N = 345678;

long long n, k;
int mod, sq, ways[N];
vector<int> primes; class my_array {
public:
int a[N << 1]; int& operator [] (long long x) {
return x <= sq ? a[x] : a[sq + n / x];
}
}; my_array sump1, sump3, sump; void add(int& x, int y) {
x += y;
if (x >= mod) {
x -= mod;
}
} void sub(int& x, int y) {
x -= y;
if (x < 0) {
x += mod;
}
} int mul(int x, int y) {
return (long long) x * y % mod;
} int qpow(int v, long long p) {
int result = 1;
for (; p; p >>= 1, v = mul(v, v)) {
if (p & 1) {
result = mul(result, v);
}
}
return result;
} void sieve(int n) {
vector<bool> is_prime(n + 1, true);
for (int i = 2; i <= n; ++i) {
if (is_prime[i]) {
primes.push_back(i);
}
for (auto x : primes) {
if (i * x > n) {
break;
}
is_prime[i * x] = false;
if (i % x == 0) {
break;
}
}
}
} void min_25_sieve() {
vector<long long> values;
for (long long i = 1; i <= n; i = n / (n / i) + 1) {
values.push_back(n / i);
}
for (auto x : values) {
sump1[x] = (x - 1 >> 2) % mod;
sump3[x] = (x + 1 >> 2) % mod;
}
for (auto p : primes) {
if (p != 2) {
if (p % 4 == 1) {
for (auto x : values) {
if (x < (long long) p * p) {
break;
}
sub(sump1[x], (sump1[x / p] - sump1[p - 1] + mod) % mod);
sub(sump3[x], (sump3[x / p] - sump3[p - 1] + mod) % mod);
}
} else {
for (auto x : values) {
if (x < (long long) p * p) {
break;
}
sub(sump1[x], (sump3[x / p] - sump3[p - 1] + mod) % mod);
sub(sump3[x], (sump1[x / p] - sump1[p - 1] + mod) % mod);
}
}
}
}
for (auto x : values) {
sump[x] = (mul(sump1[x], ways[1]) + sump3[x] + (x >= 2)) % mod;
}
} int solve(long long n, int x) {
if (n <= 1 || (x <= primes.size() && primes[x - 1] > n)) {
return 0;
}
int result = (sump[n] - (x <= primes.size() ? sump[primes[x - 1] - 1] : sump[primes.back()]) + mod) % mod;
for (int i = x; i <= primes.size(); ++i) {
long long p = primes[i - 1], q = p;
if (p * q > n) {
break;
}
for (int j = 1; p * q <= n; ++j, q *= p) {
add(result, mul(ways[p % 4 == 1 ? j : 0], solve(n / q, i + 1)));
add(result, ways[p % 4 == 1 ? j + 1 : 0]);
}
}
return result;
} int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
cin >> n >> k >> mod;
for (int i = 0; i < 34; ++i) {
ways[i] = qpow(i << 1 | 1, k);
}
sq = sqrt(n);
sq = max(sq, 2);
sieve(sq);
min_25_sieve();
cout << mul(qpow(4, k), (1 + solve(n, 1)) % mod) << '\n';
return 0;
}