题目链接
题解
复数真神奇。
一句话题意：令 (f(x)) 表示以原点 ((0, 0)) 为圆心，半径为 (x) 的圆上的整点数量，求 (sum_limits{i = 1}^N f(i)^k mod P) 的值。
令 (g(x) = frac{f(x)}{4})，那么我们需要求 (left(4^ksum_limits{i = 1}^Ng(i)^k
ight) mod 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})) 的质数对应的指数分别为 (eta_1, eta_2, cdots, eta_k)，那么有：
[egin{aligned}g(x) = left(prod_{i = 1}^m (alpha_i + 1) ight) left(prod_{i = 1}^k [eta_i 是偶数] ight)end{aligned} ]关于上述 (g(x)) 的求法的具体原因及相关理论，可参考视频 隐藏在素数规律中的 π。
包看包懂。由于 (y = x^2)，因此 (y) 的唯一分解式中任意质数的指数一定是偶数，即上述式子中后面部分 (left(prod_limits{i = 1}^k [eta_i 是偶数] ight)) 的值一定为 (1)。故有：
[egin{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 mod 4 = 1) 时，(g(p^c)) 的值为 (2c + 1) 而不是 (c + 1)。
质数的任意次幂对应的函数值我们已经能够快速计算。接下来考虑如何求所有质数对应的函数值的前缀和。考虑先分别求出 (1sim n) 内模 (4) 等于 (1) 的质数的数量以及模 (4) 等于 (3) 的质数的数量，不妨将两种质数对应的数量记为 (S_1, S_3)，那么边界值（即将 (2 sim n) 的所有数均视为质数时对应的值）显然是 (S_1(n) = leftlfloorfrac{n - 1}{4} ight floor,) (S_3(n) = leftlfloorfrac{n + 1}{4} ight floor)。考虑如何依次去掉以每一个 (geq 3) 的质数作为最小质因子的所有数的贡献。注意到在模 (4) 意义下，有：当 (a = 1) 时，(egin{cases}a imes 1 = 1 \ a imes 3 = 3end{cases})；当 (a = 3) 时，(egin{cases}a imes 1 = 3 \ a imes 3 = 1end{cases})，因此对于每个 (geq 3) 的质数 (p)，(S_1) 与 (S_3) 的更新方式如下：
令 (S(x)) 为 (1 sim x) 内的所有质数 (p) 对应的 (g(p)^k) 的值，那么有 (S(x) = S_1(x) imes 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) << '
';
return 0;
}