nt/factor_interval.hpp

View this file on GitHub
Last update: 2023-08-06 03:59:07+09:00
Include: #include "nt/factor_interval.hpp"
Link: View error logs on GitHub Actions

Depends on

nt/primetable.hpp

Verified with

test/5_atcoder/abc227g.test.cpp

Code

#include "nt/primetable.hpp"

// n が p を持つとき f(n, p) を呼ぶ
template <typename F>
void factor_interval(ll L, ll R, F f) {
  int n = R - L;
  auto primes = primetable(sqrt(R));
  vi A(n);
  iota(all(A), L);
  for (auto&& p: primes) {
    ll pp = 1;
    while (1) {
      if (pp > R / p) break;
      pp *= p;
      ll s = ceil(L, pp) * pp;
      while (s < R) {
        f(s, p);
        A[s - L] /= p;
        s += pp;
      }
    }
  }
  FOR(i, n) if (A[i] > 1) f(L + i, A[i]);
}

#line 2 "nt/primetable.hpp"

template <typename T = int>
vc<T> primetable(int LIM) {
  ++LIM;
  const int S = 32768;
  static int done = 2;
  static vc<T> primes = {2}, sieve(S + 1);

  if (done < LIM) {
    done = LIM;

    primes = {2}, sieve.assign(S + 1, 0);
    const int R = LIM / 2;
    primes.reserve(int(LIM / log(LIM) * 1.1));
    vc<pair<int, int>> cp;
    for (int i = 3; i <= S; i += 2) {
      if (!sieve[i]) {
        cp.eb(i, i * i / 2);
        for (int j = i * i; j <= S; j += 2 * i) sieve[j] = 1;
      }
    }
    for (int L = 1; L <= R; L += S) {
      array<bool, S> block{};
      for (auto& [p, idx]: cp)
        for (int i = idx; i < S + L; idx = (i += p)) block[i - L] = 1;
      FOR(i, min(S, R - L)) if (!block[i]) primes.eb((L + i) * 2 + 1);
    }
  }
  int k = LB(primes, LIM + 1);
  return {primes.begin(), primes.begin() + k};
}
#line 2 "nt/factor_interval.hpp"

// n が p を持つとき f(n, p) を呼ぶ
template <typename F>
void factor_interval(ll L, ll R, F f) {
  int n = R - L;
  auto primes = primetable(sqrt(R));
  vi A(n);
  iota(all(A), L);
  for (auto&& p: primes) {
    ll pp = 1;
    while (1) {
      if (pp > R / p) break;
      pp *= p;
      ll s = ceil(L, pp) * pp;
      while (s < R) {
        f(s, p);
        A[s - L] /= p;
        s += pp;
      }
    }
  }
  FOR(i, n) if (A[i] > 1) f(L + i, A[i]);
}