Nyaan's Library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: 素数カウント( $\mathrm{O}(\frac{N^{\frac{3}{4}}}{\log N})$・高速化版)
(multiplicative-function/prime-counting-faster.hpp)

素数の個数$\pi(N)$の高速計算

$\pi(N)$を$\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)$で計算するライブラリ。

概要

まず初めに、$f(x,n):=$($n$以下の自然数のうち$x$以下の素数を因数に持たない数の個数)とおく。

この時、

\[\pi(N) = f(\lfloor\sqrt{N}\rfloor,N) + \pi(\lfloor\sqrt{N}\rfloor) - 1\]

であるため、$f(\lfloor\sqrt{N}\rfloor,N)$を高速に求めれば$\pi(N)$が求まる。また、

\[f(x,n) = \begin{cases} f(x-1,n) & \mathrm{if}\ x\ \mathrm{is}\ \mathrm{not}\ \mathrm{prime} \newline f(x-1,n)-f(x-1,\lfloor\frac{n}{x}\rfloor) & \mathrm{otherwise} \end{cases}\]

となるが、ここで$\lfloor\frac{N}{ij}\rfloor=\lfloor\lfloor\frac{N}{i}\rfloor/j\rfloor$に着目すると、$f(\lfloor\sqrt{N}\rfloor,N)$は$f(x,n),1 \leq x \leq \sqrt{N},n = \lfloor\frac{N}{K}\rfloor$($K$は自然数)の線形和で表されるとわかる。

よって、$f(0,n)=n$を初期値として$x\leq \sqrt{N},n =\lfloor\frac{N}{K}\rfloor$を範囲とする、$x$が素数の時に更新するin-placeな動的計画法で$f(\lfloor\sqrt{N}\rfloor,N)$を求められる。計算量は

なので$\mathrm{O}\left(\frac{N}{\log N}\right)$となる。

さらに高速化する。$g(x,n) :=f(x,n) + \pi(\min(n,x))$とおくと、$g(x,n)$は$n$以下の数に$x$以下の素数でふるいを掛けた時に残ったもの($1$を含む)となる。ここから

\[g(\lfloor\sqrt{N}\rfloor,N) = \pi(N) - 1\]

であり、また$f(x,n)$の式に代入すると

\[g(x,n) = \begin{cases} g(x-1,n) & \mathrm{if}\ x\ \mathrm{is}\ \mathrm{not}\ \mathrm{prime}\ \cup\ n \lt x^2 \newline g(x-1,n)-g(x-1,\lfloor\frac{n}{x}\rfloor) + \pi(x) & \mathrm{otherwise} \end{cases}\]

を示せる。

$h(x,n) = g(x,n)-1$とおいて整理すると、

\[\pi(N) = h(\lfloor\sqrt{N}\rfloor,N)\] \[h(0,n) = n - 1\] \[h(x,n) = \begin{cases} h(x-1,n) & \mathrm{if}\ x\ \mathrm{is}\ \mathrm{not}\ \mathrm{prime}\ \cup\ n \lt x^2 \newline h(x-1,n)-h(x-1,\lfloor\frac{n}{x}\rfloor) + h(x-1,x-1) & \mathrm{otherwise} \end{cases}\]

となり、簡素な動的計画法が実行できる。計算量は$h(x,n)$の更新回数を数えればよく、

以上より全体の計算量は$\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)$となる。

補足

区間を3つに区切ってFenwick Treeを使用することで更なる計算量の改善が可能である。(ただし実装はかなり煩雑になる。)

$\mathrm{O}(N^\frac{2}{3})$での実装

なお、高速化をしたい場合は整数同士の除算をdouble型同士の除算で計算するのが最も高速化への寄与が大きい。なぜなら、64bit整数同士の除算は最悪ケースで80クロックと低速なのに対してdouble型同士は11クロックとはるかに高速であり、また$N=10^{12}$程度の大きさ同士の切り捨て除算ならばdouble型で計算しても誤差が発生しないことが規格で保証されているためである。実際にLibrary Checkerの提出で比較すると、除算を変えるだけで686ms→180msとおよそ4倍もの高速化となっている。

定数倍高速化したバージョン　提出(104ms)　他の人の提出を見ると異常高速化している人が何人も居て、魔境か？という気持ちに…

Verified with

Code

#pragma once



namespace PrimeCounting {
using i64 = long long;
static inline i64 my_div(i64 n, i64 p) { return double(n) / p; };

__attribute__((target("avx2"), optimize("O3", "unroll-loops"))) i64
prime_counting(i64 N) {
  i64 N2 = sqrt(N);
  i64 NdN2 = my_div(N, N2);

  vector<i64> hl(NdN2);
  for (int i = 1; i < NdN2; i++) hl[i] = my_div(N, i) - 1;

  vector<int> hs(N2 + 1);
  iota(begin(hs), end(hs), -1);

  for (int x = 2, pi = 0; x <= N2; ++x) {
    if (hs[x] == hs[x - 1]) continue;
    i64 x2 = i64(x) * x;
    i64 imax = min<i64>(NdN2, my_div(N, x2) + 1);
    i64 ix = x;
    for (i64 i = 1; i < imax; ++i) {
      hl[i] -= (ix < NdN2 ? hl[ix] : hs[my_div(N, ix)]) - pi;
      ix += x;
    }
    for (int n = N2; n >= x2; n--) {
      hs[n] -= hs[my_div(n, x)] - pi;
    }
    ++pi;
  }
  return hl[1];
}

}  // namespace PrimeCounting

/**
 * @brief 素数カウント( $\mathrm{O}(\frac{N^{\frac{3}{4}}}{\log N})$・高速化版)
 * @docs docs/multiplicative-function/prime-counting.md
 */
#line 2 "multiplicative-function/prime-counting-faster.hpp"



namespace PrimeCounting {
using i64 = long long;
static inline i64 my_div(i64 n, i64 p) { return double(n) / p; };

__attribute__((target("avx2"), optimize("O3", "unroll-loops"))) i64
prime_counting(i64 N) {
  i64 N2 = sqrt(N);
  i64 NdN2 = my_div(N, N2);

  vector<i64> hl(NdN2);
  for (int i = 1; i < NdN2; i++) hl[i] = my_div(N, i) - 1;

  vector<int> hs(N2 + 1);
  iota(begin(hs), end(hs), -1);

  for (int x = 2, pi = 0; x <= N2; ++x) {
    if (hs[x] == hs[x - 1]) continue;
    i64 x2 = i64(x) * x;
    i64 imax = min<i64>(NdN2, my_div(N, x2) + 1);
    i64 ix = x;
    for (i64 i = 1; i < imax; ++i) {
      hl[i] -= (ix < NdN2 ? hl[ix] : hs[my_div(N, ix)]) - pi;
      ix += x;
    }
    for (int n = N2; n >= x2; n--) {
      hs[n] -= hs[my_div(n, x)] - pi;
    }
    ++pi;
  }
  return hl[1];
}

}  // namespace PrimeCounting

/**
 * @brief 素数カウント( $\mathrm{O}(\frac{N^{\frac{3}{4}}}{\log N})$・高速化版)
 * @docs docs/multiplicative-function/prime-counting.md
 */
Back to top page