商の列挙
(math/enumerate-quotient.hpp)

• View this file on GitHub
• Last update: 2023-04-10 23:43:04+09:00
• Include: #include "math/enumerate-quotient.hpp"

Depends on

• math/isqrt.hpp

Required by

• 無平方数の数え上げ (multiplicative-function/count-square-free.hpp)
• 乗法的関数のprefix sum の列挙 (multiplicative-function/enumerate-sum-of-multiplicative-function.hpp)

Verified with

• verify/verify-unit-test/enumerate-quotient.test.cpp
• verify/verify-yosupo-math/yosupo-count-squarefrees.test.cpp
• verify/verify-yosupo-math/yosupo-counting-primes-4.test.cpp
• verify/verify-yosupo-math/yosupo-enumerate-quotient.test.cpp
• verify/verify-yosupo-math/yosupo-sum-of-totient-3.test.cpp
• verify/verify-yuki/yuki-2262.test.cpp
• verify/verify-yuki/yuki-2266.test.cpp

Code

#pragma once

#include "isqrt.hpp"

namespace EnumerateQuotientImpl {
long long fast_div(long long a, long long b) { return 1.0 * a / b; };
long long slow_div(long long a, long long b) { return a / b; };
}  // namespace EnumerateQuotientImpl

// { (q, l, r) : forall x in (l,r], floor(N/x) = q }
// を引数に取る関数f(q, l, r)を渡す。範囲が左に半開なのに注意
// 商は小さい方から走査する
template <typename T, typename F>
void enumerate_quotient(T N, const F& f) {
  T sq = isqrt(N);

#define FUNC(d)                       \
  T upper = N, quo = 0;               \
  while (upper > sq) {                \
    T thres = d(N, (++quo + 1));      \
    f(quo, thres, upper);             \
    upper = thres;                    \
  }                                   \
  while (upper > 0) {                 \
    f(d(N, upper), upper - 1, upper); \
    upper--;                          \
  }

  if (N <= 1e12) {
    FUNC(EnumerateQuotientImpl::fast_div);
  } else {
    FUNC(EnumerateQuotientImpl::slow_div);
  }
#undef FUNC
}

/**
 *  @brief 商の列挙
 */

#line 2 "math/enumerate-quotient.hpp"

#line 2 "math/isqrt.hpp"

#include <cmath>
using namespace std;

// floor(sqrt(n)) を返す (ただし n が負の場合は 0 を返す)
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;
}
#line 4 "math/enumerate-quotient.hpp"

namespace EnumerateQuotientImpl {
long long fast_div(long long a, long long b) { return 1.0 * a / b; };
long long slow_div(long long a, long long b) { return a / b; };
}  // namespace EnumerateQuotientImpl

// { (q, l, r) : forall x in (l,r], floor(N/x) = q }
// を引数に取る関数f(q, l, r)を渡す。範囲が左に半開なのに注意
// 商は小さい方から走査する
template <typename T, typename F>
void enumerate_quotient(T N, const F& f) {
  T sq = isqrt(N);

#define FUNC(d)                       \
  T upper = N, quo = 0;               \
  while (upper > sq) {                \
    T thres = d(N, (++quo + 1));      \
    f(quo, thres, upper);             \
    upper = thres;                    \
  }                                   \
  while (upper > 0) {                 \
    f(d(N, upper), upper - 1, upper); \
    upper--;                          \
  }

  if (N <= 1e12) {
    FUNC(EnumerateQuotientImpl::fast_div);
  } else {
    FUNC(EnumerateQuotientImpl::slow_div);
  }
#undef FUNC
}

/**
 *  @brief 商の列挙
 */

Back to top page