【系列】Library Checker – Math（Part 1）

Link

Counting Prime

Description

给定一个数 $N$，求 $[1,N]$ 中有多少个质数。
$1 \le N \le 10^{11}$

Solution

法一：暴力判质数。
对于每个在 $[1,N]$ 的 $i$，主要分为若干个判断流程：

• 如果 $i\le 1$，则这个数不是质数
• 枚举在 $[2,\sqrt i]$ 中的 $k$，如果 $i$ 能被 $k$ 整除，那么这个数不是质数
• 如果枚举了所有的 $k$ 仍然没有 $i$ 能被 $k$ 整除的情况，那么这个数就是质数

代码一：

#include <bits/stdc++.h>

using namespace std;

bool prime (long long x) {
    if (x <= 1) return false;
    for (long long i = 2; i * i <= x; i++)
        if (x % i == 0)
            return false;
    return true;
}

int main () {
    long long n;
    scanf("%lld", &n);
    long long cnt = 0;
    for (long long i = 1; i <= n; i++)
        if (prime(i))
            cnt++;
    printf("%lld", cnt);
    return 0;
}

时间复杂度为 $O(n\sqrt n)$，评测结果。
预算得分：$24\ pts$

法二：素数筛
线性筛质数，将判断质数的数组 $a_i$ 除了 $a_1$ 全初始化为 $\rm true$，然后如果 $a_k$ 是 $\rm true$，那么就把所有下标为 $k$ 的倍数的位置标为 $\rm false$，最后 $a_i$ 反映的就是 $i$ 是否为质数。

代码二：

#include <bits/stdc++.h>

using namespace std;

bool prime[1000000005];

void init (long long n) {
    memset(prime, true, sizeof(prime));
    prime[1] = false;
    for (long long i = 2; i <= n; i++) 
        if (prime[i]) 
            for (long long j = 2; i * j <= n; j++)
                prime[i * j] = false;
}

int main () {
    long long n;
    scanf("%lld", &n);
    init(n);
    long long cnt = 0;
    for (long long i = 1; i <= n; i++)
        if (prime[i])
            cnt++;
    printf("%lld", cnt);
    return 0;
}

时间复杂度为 $O(n \log \log n)$，我没开 $10^{11}$ 是因为这个垃圾评测器开 $10^{11}$ 就 CE。。。
评测结果
预测得分：$41\ pts$
没掉的那些分都是因为不能把数组开大导致的 RE。。。

法三：min_25 筛
这么恶心的判断质数只能用 min_25 筛了 ……
其实这个做法是我翻 AC 一览知道的（

#include <cstdio>
#include <cassert>
#include <cmath>
#include <vector>

using namespace std;

using i64 = long long;

int isqrt(i64 n) {
  return sqrtl(n);
}

__attribute__((target("avx"), optimize("O3", "unroll-loops")))
i64 prime_pi(const i64 N) {
  if (N <= 1) return 0;
  const int v = isqrt(N);
  vector<int> smalls(v + 1); for (int i = 2; i <= v; ++i) smalls[i] = (i + 1) / 2;
  int s = (v + 1) / 2;
  vector<int> roughs(s); for (int i = 0; i < s; ++i) roughs[i] = 2 * i + 1;
  vector<i64> larges(s); for (int i = 0; i < s; ++i) larges[i] = (N / (2 * i + 1) + 1) / 2;
  vector<bool> skip(v + 1);
  auto div = [] (i64 N, i64 d) -> int { return double(N) / d; };

  int pc = 0;
  for (int p = 3; p <= v; ++p) if (smalls[p] > smalls[p - 1]) {
    int q = p * p;
    pc++;
    if (i64(q) * q > N) break;
    skip[p] = true;
    for (int i = q; i <= v; i += 2 * p) skip[i] = true;
    int ns = 0;
    for (int k = 0; k < s; ++k) {
      int i = roughs[k];
      if (skip[i]) continue;
      i64 d = i64(i) * p;
      larges[ns] = larges[k] - (d <= v ? larges[smalls[d] - pc] : smalls[div(N, d)]) + pc;
      roughs[ns++] = i;
    }
    s = ns;
    for (int j = v / p; j >= p; --j) {
      int c = smalls[j] - pc;
      for (int i = j * p, e = min(i + p, v + 1); i < e; ++i) smalls[i] -= c;
    }
  }

  for (int k = 1; k < s; ++k) {
    const i64 M = N / roughs[k];
    i64 s = larges[k] - (pc + k - 1);
    for (int l = 1; l < k; ++l) {
      int p = roughs[l];
      if (i64(p) * p > M) break;
      s -= smalls[div(M, p)] - (pc + l - 1);
    }
    larges[0] -= s;
  }
  return larges[0];
}

int main(){
  i64 N;
  while (~scanf("%lld", &N)) {
    printf("%lld\n", prime_pi(N));
  }
}

应该就可以 A 了
时间复杂度：$O(n^{\frac{2}{3}})$

Enumerate Primes

Description

假设第 $i+1$ 个质数为 $p_i$，给定 $N,A,B$，求

• $\pi(N)$
• 满足 $p_{Ai+B} \le N$ 的 $p_k$ 的个数
• 满足 $p_{Ai+b} \le N$ 的所有 $p_k$

Solution

法一：暴力和线性筛（
不说了不说了，这 OJ 搞暴力必 TLE。

法二：min_25 筛

#include <algorithm>
#include <array>
#include <cassert>
#include <climits>
#include <complex>
#include <cstdio>
#include <cstring>
#include <deque>
#include <functional>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <string>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <vector>

// fast_cin.cpp

class fast_istream {
  bool is_space(char c) { return c < 0x21 || c > 0x7E; }
  template <typename T> void get_integer(T &var) {
    var = 0;
    T sign = 1;
    int c = getchar_unlocked();
    while (c < '0' || '9' < c) {
      if (c == '-') sign = -1;
      c = getchar_unlocked();
    }
    while ('0' <= c && c <= '9') {
      var = var * 10 + c - '0';
      c = getchar_unlocked();
    }
    var *= sign;
  }
  template <typename T> void get_real(T &var) {
    var = 0.0;
    T sign = 1.0;
    int c = getchar_unlocked();
    while ((c < '0' || '9' < c) && c != '.') {
      if (c == '-') sign = -1.0;
      c = getchar_unlocked();
    }
    while ('0' <= c && c <= '9') {
      var = var * 10.0 + c - '0';
      c = getchar_unlocked();
    }
    if (c == '.') {
      c = getchar_unlocked();
      T base = 1.0;
      while ('0' <= c && c <= '9') {
        base /= 10.0;
        var += base * (c - '0');
        c = getchar_unlocked();
      }
    }
    var *= sign;
  }

public:
  inline fast_istream &operator>>(int &var) {
    get_integer(var);
    return *this;
  }
  inline fast_istream &operator>>(long long &var) {
    get_integer(var);
    return *this;
  }
  inline fast_istream &operator>>(double &var) {
    get_real(var);
    return *this;
  }
  inline fast_istream &operator>>(long double &var) {
    get_real(var);
    return *this;
  }
  inline fast_istream &operator>>(std::string &var) {
    char c = getchar_unlocked();
    while (is_space(c)) {
      c = getchar_unlocked();
    }
    var = "";
    while (is_space(c)) {
      var.push_back(c);
      c = getchar_unlocked();
    }
    return *this;
  }
};

fast_istream fcin;

// fast_cout.cpp

class fast_ostream {
  char ch[32];
  template <typename T> inline void put_integer(const T &var) {
    T n = var;
    if (n == 0) {
      putchar_unlocked('0');
      return;
    }
    else if (n < 0) {
      putchar_unlocked('-');
      n = -n;
    }
    int count = 0;
    while (n != 0) {
      ch[count++] = n % 10 + '0';
      n /= 10;
    }
    while (count--) {
      putchar_unlocked(ch[count]);
    }
  }

public:
  inline fast_ostream &operator<<(const int &var) {
    put_integer(var);
    return *this;
  }
  inline fast_ostream &operator<<(const long long &var) {
    put_integer(var);
    return *this;
  }
  inline fast_ostream &operator<<(const std::string &var) {
    for (char c : var) putchar_unlocked(c);
    return *this;
  }
};

fast_ostream fcout;

const std::string fendl = "\n";

// typedef.cpp

using ll = long long;
using ld = long double;

// primes.cpp

std::vector<bool> sieve(int n) {
  std::vector<bool> is_prime(n, false);
  static const int s[16] =
    { 1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59 };
  for (int mx = 1; mx <= 60; ++mx) {
    for (int my = 1; my <= 60; my += 2) {
      int m = (4 * mx * mx + my * my) % 60;
      if (m == 1 || m == 13 || m == 17 || m == 29 ||
          m == 37 || m == 41 || m == 49 || m == 53)
        for (int x = mx; 4 * x * x < n; x += 60)
          for (int y = my, v; (v = 4 * x * x + y * y) < n; y += 60)
            is_prime[v] = !is_prime[v];
    }
  }
  for (int mx = 1; mx <= 60; mx += 2) {
    for (int my = 2; my <= 60; my += 2) {
      int m = (3 * mx * mx + my * my) % 60;
      if (m == 7 || m == 19 || m == 31 || m == 43)
        for (int x = mx; 3 * x * x < n; x += 60)
          for (int y = my, v; (v = 3 * x * x + y * y) < n; y += 60)
            is_prime[v] = !is_prime[v];
    }
  }
  for (int x = 2; 2 * x * x < n; x++) {
    for (int y = x - 1, v; y >= 1 && (v = 3 * x * x - y * y) < n; y -= 2) {
      int m = v % 60;
      if (m == 11 || m == 23 || m == 47 || m == 59)
        is_prime[v] = !is_prime[v];
    }
  }
  for (int w = 1, v; v = w / 16 * 60 + s[w % 16], v * v < n; ++w) {
    if (is_prime[v]) {
      for (ll x = 0, c; (c = v * v * ((x / 16) * 60 + s[x % 16])) < n; ++x) {
        is_prime[c] = false;
      }
    }
  }
  is_prime[2] = true;
  is_prime[3] = true;
  is_prime[5] = true;
  return is_prime;
}

// yosupo-enumerate_primes.cpp

const int t[16] =
  { 2, 3, 5 };
const int s[16] =
  { 1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59 };

int main() {
  int n, a, b;
  fcin >> n >> a >> b;
  auto is_prime = sieve(n + 1);
  int cnt = 0;
  std::vector<int> res;
  res.reserve(1000000);
  for (int i = 0; i < 3; ++i) {
    if (n >= t[i]) {
      ++cnt;
      if (--b < 0) { res.push_back(t[i]); b += a; }
    }
  }
  for (int v = 0; v < n / 60 + 1; ++v) {
    for (int w = 0; w < 16; ++w) {
      int x = v * 60 + s[w];
      if (x > n) break;
      if (is_prime[x]) {
        ++cnt;
        if (--b < 0) { res.push_back(x); b += a; }
      }
    }
  }
  int size = res.size();
  fcout << cnt << " " << size << fendl;
  for (int i = 0; i < size; ++i) {
    fcout << res[i];
    if (i == size - 1) fcout << fendl;
    else fcout << " ";
  }
  return 0;
}

也是从 AC 一览看的

至此，我们的 Prime 部分就完工啦！