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

题目链接

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^k(x)\) 显然也是积性函数。考虑使用 min_25 筛来求积性函数 \(g^k(x)\) 的前缀和。根据上面的理论，不难得出对于一个质数 \(p\)，有

\[g^k(p^c) = \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^k(p)\) 的值，那么有 \(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;
}