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

题目链接

https://loj.ac/problem/3069

题解

复数真神奇。

一句话题意：令 (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)^k = egin{cases}(2c + 1)^k, &p mod 4 = 1 \ 1, &p mod 4 eq 1 end{cases} ]

要注意的是，在算 (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) 的更新方式如下：

[egin{aligned} ext{if $p$ mod 4 = 1: }&egin{cases} {S_1(x) leftarrow S_1(x) - left(S_1 (leftlfloorfrac{x}{p} ight floor) - S_1(p - 1) ight) \S_3(x) leftarrow S_3(x) - left(S_3 (leftlfloorfrac{x}{p} ight floor) - S_3(p - 1) ight)} end{cases} \ ext{else: }&egin{cases}{S_1(x) leftarrow S_1(x) - left(S_3 (leftlfloorfrac{x}{p} ight floor) - S_3(p - 1) ight) \S_3(x) leftarrow S_3(x) - left(S_1 (leftlfloorfrac{x}{p} ight floor) - S_1(p - 1) ight)}end{cases} end{aligned} ]

令 (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;
}