dirichlet_prefix.hpp

SECTIONMath INCLUDEnoya/dirichlet_prefix.hpp

枚举 \(\left\lfloor n/i\right\rfloor\) 的所有不同值及对应整除块；为杜教筛、Min_25 类商值数组提供下标。

\[ \displaystyle q_i=\lfloor\frac{n}{i}\rfloor \]

Complexity: Time: O(n^(2/3)) arithmetic operations per convolution or division. Space: O(sqrt(n)).

AC 记录：dirichlet_convolution_and_prefix_sums, dirichlet_inverse_and_prefix_sums。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O(n^(2/3)) arithmetic operations per convolution or
/// division. Space: O(sqrt(n)).

#include "noya/integer_kth_root.hpp"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <vector>

namespace noya {

namespace dirichlet_prefix_internal {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

inline u64 integer_square_root(u64 val) {
  u64 rt = u64(std::sqrt(static_cast<long double>(val)));
  while (u128(rt + 1) * (rt + 1) <= val) {
    rt++;
  }
  while (u128(rt) * rt > val) {
    rt--;
  }
  return rt;
}

template <class T> class quotient_table {
public:
  explicit quotient_table(u64 lim) : lm_(lim) {
    assert(lm_ >= 1);
    sq_ = integer_square_root(lm_);
    sl_ = u128(sq_) * sq_ + sq_ <= lm_ ? sq_ : sq_ - 1;
    ts_ = int(sl_ + sq_ + 1);
    vs_.resize(ts_);
  }

  quotient_table(u64 lim, const std::vector<T> &pv) : quotient_table(lim) {
    assert(int(pv.size()) + 1 == ts_);
    std::copy(pv.begin(), pv.end(), vs_.begin() + 1);
  }

  int index_of(u64 val) const {
    assert(1 <= val && val <= lm_);
    int idx = val <= sl_ ? int(val) : ts_ - int(lm_ / val);
    assert(1 <= idx && idx < ts_);
    return idx;
  }

  u64 argument(int idx) const {
    assert(1 <= idx && idx < ts_);
    return u64(idx) <= sl_ ? u64(idx) : lm_ / u64(ts_ - idx);
  }

  std::vector<T> export_values() const {
    return std::vector<T>(vs_.begin() + 1, vs_.end());
  }

  u64 limit() const { return lm_; }
  u64 square_root() const { return sq_; }
  u64 small_limit() const { return sl_; }
  int table_size() const { return ts_; }

  T &operator[](int idx) { return vs_[idx]; }
  const T &operator[](int idx) const { return vs_[idx]; }

private:
  u64 lm_ = 0;
  u64 sq_ = 0;
  u64 sl_ = 0;
  int ts_ = 0;
  std::vector<T> vs_;
};

template <class T>
quotient_table<T> convolution(const quotient_table<T> &a,
                              const quotient_table<T> &b) {
  assert(a.limit() == b.limit());
  u64 lim = a.limit();
  quotient_table<T> res(lim);
  int ts = a.table_size();
  if (lim == 1) {
    res[1] = a[1] * b[1];
    return res;
  }

  std::vector<T> va(ts), vb(ts);
  for (int idx = 1; idx < ts; idx++) {
    va[idx] = a[idx] - a[idx - 1];
    vb[idx] = b[idx] - b[idx - 1];
  }

  u64 cr = integer_kth_root(lim, 3);
  u64 spl = cr * cr;
  u64 sl = a.small_limit();
  u64 sqr = a.square_root();

  for (u64 l = 1; l <= cr; l++) {
    res[a.index_of(l * l)] += va[l] * vb[l];

    if (l * (l + 1) <= sl) {
      u64 hi = sl / l;
      for (u64 r = l + 1; r <= hi; r++) {
        res[int(l * r)] += va[l] * vb[r] + va[r] * vb[l];
      }
    }

    u64 hi = std::min(spl / l, sl);
    for (u64 r = std::max(l, sl / l) + 1; r <= hi; r++) {
      res[a.index_of(l * r)] += va[l] * vb[r] + va[r] * vb[l];
    }

    if (lim / sqr <= spl / l) {
      u64 lo = lim / (spl / l + 1) + 1;
      for (u64 blk = lo; blk <= sqr; blk++) {
        int idl = ts - int(blk);
        res[a.index_of(l * (lim / blk))] += va[l] * vb[idl] + va[idl] * vb[l];
      }
    }
  }

  for (int idx = 1; idx < ts; idx++) {
    res[idx] += res[idx - 1];
  }
  for (u64 blk = 1; lim / blk > spl; blk++) {
    u64 mx = lim / blk;
    u64 rt = integer_square_root(mx);
    int ri = ts - int(blk);
    res[ri] = 0;
    for (u64 l = 1; l <= rt; l++) {
      int rin = a.index_of(mx / l);
      res[ri] += va[l] * b[rin] + vb[l] * a[rin];
    }
    res[ri] -= a[rt] * b[rt];
  }
  return res;
}

template <class T>
quotient_table<T> divide(const quotient_table<T> &num, quotient_table<T> den) {
  assert(num.limit() == den.limit());
  assert(den[1] != T(0));
  u64 lim = num.limit();
  quotient_table<T> quo(lim);
  int ts = den.table_size();
  if (lim == 1) {
    quo[1] = num[1] / den[1];
    return quo;
  }

  T ic = T(1) / den[1];
  for (int idx = 0; idx < ts; idx++) {
    den[idx] *= ic;
  }

  std::vector<T> dv(ts), qv(ts), rv(ts);
  for (int idx = 1; idx < ts; idx++) {
    dv[idx] = den[idx] - den[idx - 1];
    rv[idx] = num[idx] - num[idx - 1];
  }

  u64 cr = integer_kth_root(lim, 3);
  u64 spl = std::max(den.square_root(), cr * cr);
  qv[1] = num[1];
  for (int idx = 2; idx < ts; idx++) {
    u64 arg = den.argument(idx);
    if (arg > spl) {
      break;
    }
    qv[idx] = rv[idx] - qv[1] * dv[idx];
    if (arg * arg <= spl) {
      rv[den.index_of(arg * arg)] -= dv[idx] * qv[idx];
    }
    u64 hi = std::min<u64>(idx - 1, spl / arg);
    for (u64 oth = 2; oth <= hi; oth++) {
      rv[den.index_of(arg * oth)] -= dv[idx] * qv[oth] + dv[oth] * qv[idx];
    }
  }

  for (int idx = 1; idx < ts; idx++) {
    quo[idx] = quo[idx - 1] + qv[idx];
  }
  for (u64 blk = lim / (spl + 1); blk > 0; blk--) {
    int ri = ts - int(blk);
    u64 mx = lim / blk;
    u64 rt = integer_square_root(mx);
    quo[ri] = num[ri] - qv[1] * den[ri] + den[rt] * quo[rt];
    for (u64 l = 2; l <= rt; l++) {
      int rin = den.index_of(mx / l);
      quo[ri] -= dv[l] * quo[rin] + qv[l] * den[rin];
    }
  }
  for (int idx = 0; idx < ts; idx++) {
    quo[idx] *= ic;
  }
  return quo;
}

} // namespace dirichlet_prefix_internal

/// @brief Return Q_n={floor(n/i)} in increasing order. Consecutive equal
/// quotients are represented once.
inline std::vector<std::uint64_t> dirichlet_quotients(std::uint64_t lim) {
  dirichlet_prefix_internal::quotient_table<int> lay(lim);
  std::vector<std::uint64_t> res;
  res.reserve(lay.table_size() - 1);
  for (int idx = 1; idx < lay.table_size(); idx++) {
    res.push_back(lay.argument(idx));
  }
  return res;
}

/// @brief Given prefix sums of f and g at all increasing values in Q_n,
/// return the corresponding prefix sums of their Dirichlet convolution. Small
/// products are enumerated once around n^(1/3); for a large quotient x, the
/// hyperbola sum groups all equal floor(x/i) values through the Q_n index.
template <class T>
std::vector<T> dirichlet_convolution_prefix_sums(std::uint64_t lim,
                                                 const std::vector<T> &a,
                                                 const std::vector<T> &b) {
  using dirichlet_prefix_internal::convolution;
  using dirichlet_prefix_internal::quotient_table;
  return convolution(quotient_table<T>(lim, a), quotient_table<T>(lim, b))
      .export_values();
}

/// @brief Given prefix sums of f at all increasing values in Q_n, return the
/// prefix sums of its Dirichlet inverse. The recurrence f*g=delta is solved for
/// all small arguments first; each remaining large quotient is then recovered
/// by one grouped divisor-hyperbola equation using those known values.
template <class T>
std::vector<T> dirichlet_inverse_prefix_sums(std::uint64_t lim,
                                             const std::vector<T> &fun) {
  using dirichlet_prefix_internal::divide;
  using dirichlet_prefix_internal::quotient_table;
  quotient_table<T> ide(lim);
  for (int idx = 1; idx < ide.table_size(); idx++) {
    ide[idx] = T(1);
  }
  return divide(ide, quotient_table<T>(lim, fun)).export_values();
}

} // namespace noya

#ifndef NOYA_DIRICHLET_PREFIX_HPP
#define NOYA_DIRICHLET_PREFIX_HPP 1

/// @complexity Time: O(n^(2/3)) arithmetic operations per convolution or
/// division. Space: O(sqrt(n)).

#include "noya/integer_kth_root.hpp"

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <vector>

namespace noya {

namespace dirichlet_prefix_internal {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

inline u64 integer_square_root(u64 val) {
  u64 rt = u64(std::sqrt(static_cast<long double>(val)));
  while (u128(rt + 1) * (rt + 1) <= val) {
    rt++;
  }
  while (u128(rt) * rt > val) {
    rt--;
  }
  return rt;
}

template <class T> class quotient_table {
public:
  explicit quotient_table(u64 lim) : lm_(lim) {
    assert(lm_ >= 1);
    sq_ = integer_square_root(lm_);
    sl_ = u128(sq_) * sq_ + sq_ <= lm_ ? sq_ : sq_ - 1;
    ts_ = int(sl_ + sq_ + 1);
    vs_.resize(ts_);
  }

  quotient_table(u64 lim, const std::vector<T> &pv) : quotient_table(lim) {
    assert(int(pv.size()) + 1 == ts_);
    std::copy(pv.begin(), pv.end(), vs_.begin() + 1);
  }

  int index_of(u64 val) const {
    assert(1 <= val && val <= lm_);
    int idx = val <= sl_ ? int(val) : ts_ - int(lm_ / val);
    assert(1 <= idx && idx < ts_);
    return idx;
  }

  u64 argument(int idx) const {
    assert(1 <= idx && idx < ts_);
    return u64(idx) <= sl_ ? u64(idx) : lm_ / u64(ts_ - idx);
  }

  std::vector<T> export_values() const {
    return std::vector<T>(vs_.begin() + 1, vs_.end());
  }

  u64 limit() const { return lm_; }
  u64 square_root() const { return sq_; }
  u64 small_limit() const { return sl_; }
  int table_size() const { return ts_; }

  T &operator[](int idx) { return vs_[idx]; }
  const T &operator[](int idx) const { return vs_[idx]; }

private:
  u64 lm_ = 0;
  u64 sq_ = 0;
  u64 sl_ = 0;
  int ts_ = 0;
  std::vector<T> vs_;
};

template <class T>
quotient_table<T> convolution(const quotient_table<T> &a,
                              const quotient_table<T> &b) {
  assert(a.limit() == b.limit());
  u64 lim = a.limit();
  quotient_table<T> res(lim);
  int ts = a.table_size();
  if (lim == 1) {
    res[1] = a[1] * b[1];
    return res;
  }

  std::vector<T> va(ts), vb(ts);
  for (int idx = 1; idx < ts; idx++) {
    va[idx] = a[idx] - a[idx - 1];
    vb[idx] = b[idx] - b[idx - 1];
  }

  u64 cr = integer_kth_root(lim, 3);
  u64 spl = cr * cr;
  u64 sl = a.small_limit();
  u64 sqr = a.square_root();

  for (u64 l = 1; l <= cr; l++) {
    res[a.index_of(l * l)] += va[l] * vb[l];

    if (l * (l + 1) <= sl) {
      u64 hi = sl / l;
      for (u64 r = l + 1; r <= hi; r++) {
        res[int(l * r)] += va[l] * vb[r] + va[r] * vb[l];
      }
    }

    u64 hi = std::min(spl / l, sl);
    for (u64 r = std::max(l, sl / l) + 1; r <= hi; r++) {
      res[a.index_of(l * r)] += va[l] * vb[r] + va[r] * vb[l];
    }

    if (lim / sqr <= spl / l) {
      u64 lo = lim / (spl / l + 1) + 1;
      for (u64 blk = lo; blk <= sqr; blk++) {
        int idl = ts - int(blk);
        res[a.index_of(l * (lim / blk))] += va[l] * vb[idl] + va[idl] * vb[l];
      }
    }
  }

  for (int idx = 1; idx < ts; idx++) {
    res[idx] += res[idx - 1];
  }
  for (u64 blk = 1; lim / blk > spl; blk++) {
    u64 mx = lim / blk;
    u64 rt = integer_square_root(mx);
    int ri = ts - int(blk);
    res[ri] = 0;
    for (u64 l = 1; l <= rt; l++) {
      int rin = a.index_of(mx / l);
      res[ri] += va[l] * b[rin] + vb[l] * a[rin];
    }
    res[ri] -= a[rt] * b[rt];
  }
  return res;
}

template <class T>
quotient_table<T> divide(const quotient_table<T> &num, quotient_table<T> den) {
  assert(num.limit() == den.limit());
  assert(den[1] != T(0));
  u64 lim = num.limit();
  quotient_table<T> quo(lim);
  int ts = den.table_size();
  if (lim == 1) {
    quo[1] = num[1] / den[1];
    return quo;
  }

  T ic = T(1) / den[1];
  for (int idx = 0; idx < ts; idx++) {
    den[idx] *= ic;
  }

  std::vector<T> dv(ts), qv(ts), rv(ts);
  for (int idx = 1; idx < ts; idx++) {
    dv[idx] = den[idx] - den[idx - 1];
    rv[idx] = num[idx] - num[idx - 1];
  }

  u64 cr = integer_kth_root(lim, 3);
  u64 spl = std::max(den.square_root(), cr * cr);
  qv[1] = num[1];
  for (int idx = 2; idx < ts; idx++) {
    u64 arg = den.argument(idx);
    if (arg > spl) {
      break;
    }
    qv[idx] = rv[idx] - qv[1] * dv[idx];
    if (arg * arg <= spl) {
      rv[den.index_of(arg * arg)] -= dv[idx] * qv[idx];
    }
    u64 hi = std::min<u64>(idx - 1, spl / arg);
    for (u64 oth = 2; oth <= hi; oth++) {
      rv[den.index_of(arg * oth)] -= dv[idx] * qv[oth] + dv[oth] * qv[idx];
    }
  }

  for (int idx = 1; idx < ts; idx++) {
    quo[idx] = quo[idx - 1] + qv[idx];
  }
  for (u64 blk = lim / (spl + 1); blk > 0; blk--) {
    int ri = ts - int(blk);
    u64 mx = lim / blk;
    u64 rt = integer_square_root(mx);
    quo[ri] = num[ri] - qv[1] * den[ri] + den[rt] * quo[rt];
    for (u64 l = 2; l <= rt; l++) {
      int rin = den.index_of(mx / l);
      quo[ri] -= dv[l] * quo[rin] + qv[l] * den[rin];
    }
  }
  for (int idx = 0; idx < ts; idx++) {
    quo[idx] *= ic;
  }
  return quo;
}

} // namespace dirichlet_prefix_internal

/// @brief Return Q_n={floor(n/i)} in increasing order. Consecutive equal
/// quotients are represented once.
inline std::vector<std::uint64_t> dirichlet_quotients(std::uint64_t lim) {
  dirichlet_prefix_internal::quotient_table<int> lay(lim);
  std::vector<std::uint64_t> res;
  res.reserve(lay.table_size() - 1);
  for (int idx = 1; idx < lay.table_size(); idx++) {
    res.push_back(lay.argument(idx));
  }
  return res;
}

/// @brief Given prefix sums of f and g at all increasing values in Q_n,
/// return the corresponding prefix sums of their Dirichlet convolution. Small
/// products are enumerated once around n^(1/3); for a large quotient x, the
/// hyperbola sum groups all equal floor(x/i) values through the Q_n index.
template <class T>
std::vector<T> dirichlet_convolution_prefix_sums(std::uint64_t lim,
                                                 const std::vector<T> &a,
                                                 const std::vector<T> &b) {
  using dirichlet_prefix_internal::convolution;
  using dirichlet_prefix_internal::quotient_table;
  return convolution(quotient_table<T>(lim, a), quotient_table<T>(lim, b))
      .export_values();
}

/// @brief Given prefix sums of f at all increasing values in Q_n, return the
/// prefix sums of its Dirichlet inverse. The recurrence f*g=delta is solved for
/// all small arguments first; each remaining large quotient is then recovered
/// by one grouped divisor-hyperbola equation using those known values.
template <class T>
std::vector<T> dirichlet_inverse_prefix_sums(std::uint64_t lim,
                                             const std::vector<T> &fun) {
  using dirichlet_prefix_internal::divide;
  using dirichlet_prefix_internal::quotient_table;
  quotient_table<T> ide(lim);
  for (int idx = 1; idx < ide.table_size(); idx++) {
    ide[idx] = T(1);
  }
  return divide(ide, quotient_table<T>(lim, fun)).export_values();
}

} // namespace noya

#endif // NOYA_DIRICHLET_PREFIX_HPP

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <vector>

/// @complexity Time: O(n^(2/3)) arithmetic operations per convolution or
/// division. Space: O(sqrt(n)).

/// @complexity Time: O(k ceil(64/k)) = O(64) bounded multiplications for a
/// 64-bit input.  Space: O(1).

namespace noya {

/// @brief Return floor(val^(1/exp)) for an unsigned 64-bit integer.
/// Binary search uses an exponent-dependent upper bound.  The comparison
/// checks prd <= val / bas before multiplying by bas, so it is
/// exact and never relies on floating-point rounding or overflowing products.
inline std::uint64_t integer_kth_root(std::uint64_t val, int exp) {
  assert(1 <= exp && exp <= 64);
  if (exp == 1 || val <= 1) {
    return val;
  }
  auto chk = [&](std::uint64_t bas) {
    std::uint64_t prd = 1;
    for (int cnt = 0; cnt < exp; cnt++) {
      if (bas != 0 && prd > val / bas) {
        return false;
      }
      prd *= bas;
    }
    return true;
  };
  int ub = (64 + exp - 1) / exp;
  std::uint64_t low = 0;
  std::uint64_t hig = std::uint64_t(1) << ub;
  while (hig - low > 1) {
    std::uint64_t mid = low + (hig - low) / 2;
    (chk(mid) ? low : hig) = mid;
  }
  return low;
}

} // namespace noya

namespace noya {

namespace dirichlet_prefix_internal {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

inline u64 integer_square_root(u64 val) {
  u64 rt = u64(std::sqrt(static_cast<long double>(val)));
  while (u128(rt + 1) * (rt + 1) <= val) {
    rt++;
  }
  while (u128(rt) * rt > val) {
    rt--;
  }
  return rt;
}

template <class T> class quotient_table {
public:
  explicit quotient_table(u64 lim) : lm_(lim) {
    assert(lm_ >= 1);
    sq_ = integer_square_root(lm_);
    sl_ = u128(sq_) * sq_ + sq_ <= lm_ ? sq_ : sq_ - 1;
    ts_ = int(sl_ + sq_ + 1);
    vs_.resize(ts_);
  }

  quotient_table(u64 lim, const std::vector<T> &pv) : quotient_table(lim) {
    assert(int(pv.size()) + 1 == ts_);
    std::copy(pv.begin(), pv.end(), vs_.begin() + 1);
  }

  int index_of(u64 val) const {
    assert(1 <= val && val <= lm_);
    int idx = val <= sl_ ? int(val) : ts_ - int(lm_ / val);
    assert(1 <= idx && idx < ts_);
    return idx;
  }

  u64 argument(int idx) const {
    assert(1 <= idx && idx < ts_);
    return u64(idx) <= sl_ ? u64(idx) : lm_ / u64(ts_ - idx);
  }

  std::vector<T> export_values() const {
    return std::vector<T>(vs_.begin() + 1, vs_.end());
  }

  u64 limit() const { return lm_; }
  u64 square_root() const { return sq_; }
  u64 small_limit() const { return sl_; }
  int table_size() const { return ts_; }

  T &operator[](int idx) { return vs_[idx]; }
  const T &operator[](int idx) const { return vs_[idx]; }

private:
  u64 lm_ = 0;
  u64 sq_ = 0;
  u64 sl_ = 0;
  int ts_ = 0;
  std::vector<T> vs_;
};

template <class T>
quotient_table<T> convolution(const quotient_table<T> &a,
                              const quotient_table<T> &b) {
  assert(a.limit() == b.limit());
  u64 lim = a.limit();
  quotient_table<T> res(lim);
  int ts = a.table_size();
  if (lim == 1) {
    res[1] = a[1] * b[1];
    return res;
  }

  std::vector<T> va(ts), vb(ts);
  for (int idx = 1; idx < ts; idx++) {
    va[idx] = a[idx] - a[idx - 1];
    vb[idx] = b[idx] - b[idx - 1];
  }

  u64 cr = integer_kth_root(lim, 3);
  u64 spl = cr * cr;
  u64 sl = a.small_limit();
  u64 sqr = a.square_root();

  for (u64 l = 1; l <= cr; l++) {
    res[a.index_of(l * l)] += va[l] * vb[l];

    if (l * (l + 1) <= sl) {
      u64 hi = sl / l;
      for (u64 r = l + 1; r <= hi; r++) {
        res[int(l * r)] += va[l] * vb[r] + va[r] * vb[l];
      }
    }

    u64 hi = std::min(spl / l, sl);
    for (u64 r = std::max(l, sl / l) + 1; r <= hi; r++) {
      res[a.index_of(l * r)] += va[l] * vb[r] + va[r] * vb[l];
    }

    if (lim / sqr <= spl / l) {
      u64 lo = lim / (spl / l + 1) + 1;
      for (u64 blk = lo; blk <= sqr; blk++) {
        int idl = ts - int(blk);
        res[a.index_of(l * (lim / blk))] += va[l] * vb[idl] + va[idl] * vb[l];
      }
    }
  }

  for (int idx = 1; idx < ts; idx++) {
    res[idx] += res[idx - 1];
  }
  for (u64 blk = 1; lim / blk > spl; blk++) {
    u64 mx = lim / blk;
    u64 rt = integer_square_root(mx);
    int ri = ts - int(blk);
    res[ri] = 0;
    for (u64 l = 1; l <= rt; l++) {
      int rin = a.index_of(mx / l);
      res[ri] += va[l] * b[rin] + vb[l] * a[rin];
    }
    res[ri] -= a[rt] * b[rt];
  }
  return res;
}

template <class T>
quotient_table<T> divide(const quotient_table<T> &num, quotient_table<T> den) {
  assert(num.limit() == den.limit());
  assert(den[1] != T(0));
  u64 lim = num.limit();
  quotient_table<T> quo(lim);
  int ts = den.table_size();
  if (lim == 1) {
    quo[1] = num[1] / den[1];
    return quo;
  }

  T ic = T(1) / den[1];
  for (int idx = 0; idx < ts; idx++) {
    den[idx] *= ic;
  }

  std::vector<T> dv(ts), qv(ts), rv(ts);
  for (int idx = 1; idx < ts; idx++) {
    dv[idx] = den[idx] - den[idx - 1];
    rv[idx] = num[idx] - num[idx - 1];
  }

  u64 cr = integer_kth_root(lim, 3);
  u64 spl = std::max(den.square_root(), cr * cr);
  qv[1] = num[1];
  for (int idx = 2; idx < ts; idx++) {
    u64 arg = den.argument(idx);
    if (arg > spl) {
      break;
    }
    qv[idx] = rv[idx] - qv[1] * dv[idx];
    if (arg * arg <= spl) {
      rv[den.index_of(arg * arg)] -= dv[idx] * qv[idx];
    }
    u64 hi = std::min<u64>(idx - 1, spl / arg);
    for (u64 oth = 2; oth <= hi; oth++) {
      rv[den.index_of(arg * oth)] -= dv[idx] * qv[oth] + dv[oth] * qv[idx];
    }
  }

  for (int idx = 1; idx < ts; idx++) {
    quo[idx] = quo[idx - 1] + qv[idx];
  }
  for (u64 blk = lim / (spl + 1); blk > 0; blk--) {
    int ri = ts - int(blk);
    u64 mx = lim / blk;
    u64 rt = integer_square_root(mx);
    quo[ri] = num[ri] - qv[1] * den[ri] + den[rt] * quo[rt];
    for (u64 l = 2; l <= rt; l++) {
      int rin = den.index_of(mx / l);
      quo[ri] -= dv[l] * quo[rin] + qv[l] * den[rin];
    }
  }
  for (int idx = 0; idx < ts; idx++) {
    quo[idx] *= ic;
  }
  return quo;
}

} // namespace dirichlet_prefix_internal

/// @brief Return Q_n={floor(n/i)} in increasing order. Consecutive equal
/// quotients are represented once.
inline std::vector<std::uint64_t> dirichlet_quotients(std::uint64_t lim) {
  dirichlet_prefix_internal::quotient_table<int> lay(lim);
  std::vector<std::uint64_t> res;
  res.reserve(lay.table_size() - 1);
  for (int idx = 1; idx < lay.table_size(); idx++) {
    res.push_back(lay.argument(idx));
  }
  return res;
}

/// @brief Given prefix sums of f and g at all increasing values in Q_n,
/// return the corresponding prefix sums of their Dirichlet convolution. Small
/// products are enumerated once around n^(1/3); for a large quotient x, the
/// hyperbola sum groups all equal floor(x/i) values through the Q_n index.
template <class T>
std::vector<T> dirichlet_convolution_prefix_sums(std::uint64_t lim,
                                                 const std::vector<T> &a,
                                                 const std::vector<T> &b) {
  using dirichlet_prefix_internal::convolution;
  using dirichlet_prefix_internal::quotient_table;
  return convolution(quotient_table<T>(lim, a), quotient_table<T>(lim, b))
      .export_values();
}

/// @brief Given prefix sums of f at all increasing values in Q_n, return the
/// prefix sums of its Dirichlet inverse. The recurrence f*g=delta is solved for
/// all small arguments first; each remaining large quotient is then recovered
/// by one grouped divisor-hyperbola equation using those known values.
template <class T>
std::vector<T> dirichlet_inverse_prefix_sums(std::uint64_t lim,
                                             const std::vector<T> &fun) {
  using dirichlet_prefix_internal::divide;
  using dirichlet_prefix_internal::quotient_table;
  quotient_table<T> ide(lim);
  for (int idx = 1; idx < ide.table_size(); idx++) {
    ide[idx] = T(1);
  }
  return divide(ide, quotient_table<T>(lim, fun)).export_values();
}

} // namespace noya