Lucy DP
(number-theory/lucy-dp.hpp)

View this file on GitHub
Last update: 2026-02-28 01:08:20+09:00
Include: #include "number-theory/lucy-dp.hpp"

正整数 $N$ が与えられたとき，完全乗法的関数で prefix の和が簡単に計算できる $f$ について $\sum_{p:\text{prime},p\leq N}f(p)$ を $O(N^{2/3}/(\log N)^{1/3})$ 時間で求められる．

以下の情報を与える必要がある．

正整数 $k$ に対する $\sum_{i=1}^{\lfloor N/k\rfloor}f(i)$ の値（$O(\sqrt{N})$ 個）．
$\sqrt{N}$ 以下の素数に対する $f$ の値．
- 実際は上の総和から計算できる．

例題：

https://atcoder.jp/contests/jsc2024-final/tasks/jsc2024_final_b

$f$ の例：

$f(n)=n^k$：$k=0$ のとき $\sum_{p:\text{prime}}f(p)$ は $N$ 以下の素数の個数である．
Dirichlet 指標：$4$ で割ったあまりが $1$ の素数の和なども求められる．

アルゴリズム：$O(N^{3/4}/\log N)$

この節は [1] を大幅に参考にしている．

$S(n)=\sum_{p:\text{prime},p\leq n}f(p)$ とおく．

エラトステネスの篩の要領で計算する．

$S_x(n)$ を $\sum_{i=2}^{n}f(i)$ から $x$ 以下の素因数をもつ合成数についての $f$ の値を減じたものとする．$S_x$ は $x$ 以下の数について篩った時点での累積和であり，$x\geq\lfloor\sqrt{n}\rfloor$ ならば $S(n)=S_x(n)$ である．また $S_1(n)=\sum_{i=2}^{n}f(i)$ である．

以下最小素因数を $\operatorname{lpf}(n)$ で表すことにする．

$S_x$ を $S_{x-1}$ で表す． $x$ が素数でないならば篩うことはなく $S_x=S_{x-1}$ である． $x$ が素数であるとき，エラトステネスの篩では $\operatorname{lpf}(i)=x$ かつ $i\neq x$ なる $i$ が篩われる．特に $x^2$ 未満の数は篩われないため以下が成り立つ． $S_x(n)=\begin{cases} S_{x-1}(n)-\left(S_{x-1}(\lfloor n/x \rfloor)-S_{x-1}(x-1)\right)f(x)&(x:\text{prime}, x^2\leq n)\\ S_{x-1}(n)&(\text{otherwise}) \end{cases}$

$S(N)$ を求めるためには $S_x(\lfloor N/\rfloor)$ の形の部分だけ管理すればよく，$n$ について降順に計算することで in-place に更新できる． $\lfloor N/\rfloor$ のとりうる値の種類数は $O(\sqrt{N})$ であるため空間計算量は $O(\sqrt{N})$ になる．

時間計算量は $\lfloor N/*\rfloor$ の取り得る値すべてについての，その平方根以下の素数の個数の総和のオーダーである．素数定理 $\pi(n)\approx\operatorname{Li}(n),\quad \operatorname{Li}(x)=\int_{2}^{x}\frac{1}{\log t}dt\approx\frac{x}{\log x}$ を用いれば次のように見積もれる．

\[\begin{align*} \sum_{n=2}^{\lfloor\sqrt{N}\rfloor}\pi(\sqrt{n}) +\sum_{n=1}^{\lfloor\sqrt{N}\rfloor}\pi\left(\sqrt{\frac{N}{n}}\right) %&\approx\int_{2}^{\sqrt{N}}\operatorname{Li}(\sqrt{x})dx+\int_{1}^{\sqrt{N}}\operatorname{Li}\left(\sqrt{N/x}\right)dx\\ &\approx\int_{2}^{\sqrt{N}}\frac{\sqrt{x}}{\log\sqrt{x}}dx+\int_{1}^{\sqrt{N}}\frac{\sqrt{N/x}}{\log\sqrt{N/x}}dx\\ &\leq\int_{2}^{\sqrt{N}}\frac{N^{1/4}}{\frac{1}{2}\log x}dx+\int_{1}^{\sqrt{N}}\frac{\sqrt{N/x}}{\log N^{1/4}}dx\\ &\approx N^{1/4}\cdot 2\operatorname{Li}(\sqrt{N})+\frac{1}{\log N^{1/4}}\cdot\frac{N^{3/4}}{2}\\ &=O\left(\frac{N^{3/4}}{\log N}\right) \end{align*}\]

以上の内容をアルゴリズムの形でまとめておく． $Q={\lfloor N/i\rfloor \mid 1\leq i\leq N}$ とする．

$S[q]=\sum_{i=2}^{q}f(i)$ と初期化．

$p^2\leq N$ なる素数 $p$ に対して以下を行う．

$p^2\leq q$ なる $q\in Q$ について降順に $S[q]$ を $S[q]-(S[\lfloor q/p\rfloor]-S[p-1])f(p)$ で更新する．

return $S$．

高速化：$\tilde O(N^{2/3})$

$\tilde O(N^{2/3})$ 時間に高速化できる．

[2] では $O(N^{2/3}(\log N)^{1/3})$ 時間 / $\tilde O(N^{2/3})$ 空間のアルゴリズムが紹介されている．そこまで大きな改善ではなく空間計算量は悪化しているが実装はシンプルである．

$O(N^{2/3}/(\log N)^{1/3})$ 時間

[1] の内容．

計算量を評価しなおす．素数 $p$ の計算量への寄与は $\max(\sqrt{N}-p^2,0)+\min(N/p^2,\sqrt{N})$ であるので，時間計算量は $\int_{2}^{\sqrt{N}}\left(\max(\sqrt{N}-x^2,0)+\min(N/x^2,\sqrt{N})\right)\frac{1}{\log x} dx$ と近似できる．ここで $\frac{1}{\log x}$ は $\frac{d\operatorname{Li}(x)}{dx}=\frac{1}{\log x}$ に由来する．

$p\leq N^{1/6}$ なる $p$ に対応する部分は以下のように評価できる． $\begin{align*} &\int_{2}^{N^{1/6}}\left(\max(\sqrt{N}-x^2,0)+\min(N/x^2,\sqrt{N})\right)\frac{1}{\log x}dx\\ &=\int_{2}^{N^{1/6}}\left((\sqrt{N}-x^2)+\sqrt{N}\right)\frac{1}{\log x}dx\\ &\leq 2\sqrt{N}\operatorname{Li}(N^{1/6})\\ &=O\left(\frac{N^{2/3}}{\log N}\right) \end{align*}$

また $N^{1/3}\leq p\leq N^{1/2}$ なる $p$ に対応する部分は以下のように評価できる． $\begin{align*} &\int_{N^{1/3}}^{N^{1/2}}\left(\max(\sqrt{N}-x^2,0)+\min(N/x^2,\sqrt{N})\right)\frac{1}{\log x}dx\\ &=\int_{N^{1/3}}^{N^{1/2}}\frac{N}{x^2}\cdot\frac{1}{\log x}dx\\ &\leq\int_{N^{1/3}}^{N^{1/2}}\frac{N}{x^2}\cdot\frac{1}{\log N^{1/3}}dx\\ &=O\left(\frac{N^{2/3}}{\log N}\right)\\ \end{align*}$

従ってこれらの範囲の $p$ については普通に Lucy DP の計算をしても $O\left(\frac{N^{2/3}}{\log N}\right)$ 時間で計算できる．

$N^{1/6}\lt p\lt N^{1/3}$ なる $p$ について考える． $q\geq N^{2/3}$ なる $q\in Q$ は $O(N^{1/3})$ 個なのでその遷移は愚直に行うことにして，$N^{1/3}\lt p^2\leq q\lt N^{2/3}$ なる $q$ についての遷移を再考する．

(TODO: やる)

$S[q]=\sum_{i=2}^{q}f(i)$ と初期化．

$p^2\leq N$ なる素数 $p$ に対して以下を行う．

$p^2\leq q$ なる $q\in Q$ について降順に $S[q]$ を $S[q]-(S[\lfloor q/p\rfloor]-S[p-1]-1)f(p)$ で更新する．

return $S$．

参考資料

[1] : https://rsk0315.github.io/slides/prime-counting.pdf
[2] : https://gbroxey.github.io/blog/2023/04/09/lucy-fenwick.html

Depends on

Code

#pragma once

#include "math/util.hpp"
#include "math/prime-sieve.hpp"

template <class T>
pair<vector<long long>, vector<T>> LucyDP(long long n, function<T(ll)> point_value, function<T(ll)> prefix_sum) {
  using ll = long long;
  const ll sq = Math::isqrt(n);

  vector<ll> qs(sq * 2 - (sq * sq == n));
  for (int i = 0; i < sq; i++) qs[i] = i + 1;
  for (int i = 0; i < sq; i++) qs[qs.size() - 1 - i] = n / (i + 1);
  vector<T> s(qs.size());
  auto v1 = point_value(1);
  for (int i = 0; i < qs.size(); i++) s[i] = prefix_sum(qs[i]) - v1;

  auto ps = PrimeSieve(sq);
  for (ll p : ps) {
    auto v = point_value(p);
    for (int i = (int)qs.size() - 1; i >= 0; i--) {
      ll q = qs[i];
      if (p * p > q) break;
      ll x = q / p;
      int j = x <= sq ? x - 1 : (int)s.size() - n / x;
      s[i] -= (s[j] - s[p - 2]) * v;
    }
  }

  return {qs, s};
}

/**
 * @brief Lucy DP
 * @docs docs/number-theory/lucy-dp.md
 */

#line 2 "number-theory/lucy-dp.hpp"

#line 2 "math/util.hpp"

namespace Math {
template <class T>
T safe_mod(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  a %= b;
  return a >= 0 ? a : a + b;
}
template <class T>
T floor(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  return a >= 0 ? a / b : (a + 1) / b - 1;
}
template <class T>
T ceil(T a, T b) {
  assert(b != 0);
  if (b < 0) a = -a, b = -b;
  return a > 0 ? (a - 1) / b + 1 : a / b;
}
long long isqrt(long long n) {
  if (n <= 0) return 0;
  long long x = sqrt(n);
  while ((x + 1) * (x + 1) <= n) x++;
  while (x * x > n) x--;
  return x;
}
// return g=gcd(a,b)
// a*x+b*y=g
// - b!=0 -> 0<=x<|b|/g
// - b=0  -> ax=g
template <class T>
T ext_gcd(T a, T b, T& x, T& y) {
  T a0 = a, b0 = b;
  bool sgn_a = a < 0, sgn_b = b < 0;
  if (sgn_a) a = -a;
  if (sgn_b) b = -b;
  if (b == 0) {
    x = sgn_a ? -1 : 1;
    y = 0;
    return a;
  }
  T x00 = 1, x01 = 0, x10 = 0, x11 = 1;
  while (b != 0) {
    T q = a / b, r = a - b * q;
    x00 -= q * x01;
    x10 -= q * x11;
    swap(x00, x01);
    swap(x10, x11);
    a = b, b = r;
  }
  x = x00, y = x10;
  if (sgn_a) x = -x;
  if (sgn_b) y = -y;
  if (b0 != 0) {
    a0 /= a, b0 /= a;
    if (b0 < 0) a0 = -a0, b0 = -b0;
    T q = x >= 0 ? x / b0 : (x + 1) / b0 - 1;
    x -= b0 * q;
    y += a0 * q;
  }
  return a;
}
constexpr long long inv_mod(long long x, long long m) {
  x %= m;
  if (x < 0) x += m;
  long long a = m, b = x;
  long long y0 = 0, y1 = 1;
  while (b > 0) {
    long long q = a / b;
    swap(a -= q * b, b);
    swap(y0 -= q * y1, y1);
  }
  if (y0 < 0) y0 += m / a;
  return y0;
}
long long pow_mod(long long x, long long n, long long m) {
  x = (x % m + m) % m;
  long long y = 1;
  while (n) {
    if (n & 1) y = y * x % m;
    x = x * x % m;
    n >>= 1;
  }
  return y;
}
constexpr long long pow_mod_constexpr(long long x, long long n, int m) {
  if (m == 1) return 0;
  unsigned int _m = (unsigned int)(m);
  unsigned long long r = 1;
  unsigned long long y = x % m;
  if (y >= m) y += m;
  while (n) {
    if (n & 1) r = (r * y) % _m;
    y = (y * y) % _m;
    n >>= 1;
  }
  return r;
}
constexpr bool is_prime_constexpr(int n) {
  if (n <= 1) return false;
  if (n == 2 || n == 7 || n == 61) return true;
  if (n % 2 == 0) return false;
  long long d = n - 1;
  while (d % 2 == 0) d /= 2;
  constexpr long long bases[3] = {2, 7, 61};
  for (long long a : bases) {
    long long t = d;
    long long y = pow_mod_constexpr(a, t, n);
    while (t != n - 1 && y != 1 && y != n - 1) {
      y = y * y % n;
      t <<= 1;
    }
    if (y != n - 1 && t % 2 == 0) {
      return false;
    }
  }
  return true;
}
template <int n>
constexpr bool is_prime = is_prime_constexpr(n);
};  // namespace Math
#line 2 "math/prime-sieve.hpp"

vector<int> PrimeSieve(int n) {
  vector<bool> f(n + 1, false);
  for (int p = 2; p * p <= n; p += (p & 1) + 1) {
    if (f[p]) continue;
    for (int i = p * p; i <= n; i += p) f[i] = true;
  }
  vector<int> ps;
  for (int p = 2; p <= n; p += (p & 1) + 1)
    if (!f[p]) ps.push_back(p);
  return ps;
}
/**
 * @brief 素数篩
 */
#line 5 "number-theory/lucy-dp.hpp"

template <class T>
pair<vector<long long>, vector<T>> LucyDP(long long n, function<T(ll)> point_value, function<T(ll)> prefix_sum) {
  using ll = long long;
  const ll sq = Math::isqrt(n);

  vector<ll> qs(sq * 2 - (sq * sq == n));
  for (int i = 0; i < sq; i++) qs[i] = i + 1;
  for (int i = 0; i < sq; i++) qs[qs.size() - 1 - i] = n / (i + 1);
  vector<T> s(qs.size());
  auto v1 = point_value(1);
  for (int i = 0; i < qs.size(); i++) s[i] = prefix_sum(qs[i]) - v1;

  auto ps = PrimeSieve(sq);
  for (ll p : ps) {
    auto v = point_value(p);
    for (int i = (int)qs.size() - 1; i >= 0; i--) {
      ll q = qs[i];
      if (p * p > q) break;
      ll x = q / p;
      int j = x <= sq ? x - 1 : (int)s.size() - n / x;
      s[i] -= (s[j] - s[p - 2]) * v;
    }
  }

  return {qs, s};
}

/**
 * @brief Lucy DP
 * @docs docs/number-theory/lucy-dp.md
 */

Lucy DP (number-theory/lucy-dp.hpp)