题目链接
题解
复数真神奇。
一句话题意：令 \(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)\) 时，实际是对 \((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\) 的更新方式如下：
\]
令 \(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;
}