Skip to content

integer_kth_root.hpp

SECTIONMath INCLUDEnoya/integer_kth_root.hpp

Return floor(value^(1/exponent)) for an unsigned 64-bit integer. Binary search uses an exponent-dependent upper bound. The comparison multiplies only while the next factor is at most value/current, so it is exact and never relies on floating-point rounding or overflowing products.

Verified by kth_root_integer.

\[ \displaystyle r=\lfloor x^{1/k}\rfloor \]

Implementation

View on GitHub

#ifndef NOYA_INTEGER_KTH_ROOT_HPP
#define NOYA_INTEGER_KTH_ROOT_HPP 1

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

#include <cassert>
#include <cstdint>

namespace noya {

/// @brief Return floor(value^(1/exponent)) for an unsigned 64-bit integer.
/// Binary search uses an exponent-dependent upper bound.  The comparison
/// multiplies only while the next factor is at most value/current, so it is
/// exact and never relies on floating-point rounding or overflowing products.
inline std::uint64_t integer_kth_root(std::uint64_t value, int exponent) {
  assert(1 <= exponent && exponent <= 64);
  if (exponent == 1 || value <= 1) {
    return value;
  }
  auto power_at_most = [&](std::uint64_t base) {
    std::uint64_t product = 1;
    for (int count = 0; count < exponent; count++) {
      if (base != 0 && product > value / base) {
        return false;
      }
      product *= base;
    }
    return true;
  };
  int upper_bit = (64 + exponent - 1) / exponent;
  std::uint64_t low = 0;
  std::uint64_t high = std::uint64_t(1) << upper_bit;
  while (high - low > 1) {
    std::uint64_t middle = low + (high - low) / 2;
    (power_at_most(middle) ? low : high) = middle;
  }
  return low;
}

} // namespace noya

#endif // NOYA_INTEGER_KTH_ROOT_HPP
#include <cassert>
#include <cstdint>

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

namespace noya {

/// @brief Return floor(value^(1/exponent)) for an unsigned 64-bit integer.
/// Binary search uses an exponent-dependent upper bound.  The comparison
/// multiplies only while the next factor is at most value/current, so it is
/// exact and never relies on floating-point rounding or overflowing products.
inline std::uint64_t integer_kth_root(std::uint64_t value, int exponent) {
  assert(1 <= exponent && exponent <= 64);
  if (exponent == 1 || value <= 1) {
    return value;
  }
  auto power_at_most = [&](std::uint64_t base) {
    std::uint64_t product = 1;
    for (int count = 0; count < exponent; count++) {
      if (base != 0 && product > value / base) {
        return false;
      }
      product *= base;
    }
    return true;
  };
  int upper_bit = (64 + exponent - 1) / exponent;
  std::uint64_t low = 0;
  std::uint64_t high = std::uint64_t(1) << upper_bit;
  while (high - low > 1) {
    std::uint64_t middle = low + (high - low) / 2;
    (power_at_most(middle) ? low : high) = middle;
  }
  return low;
}

} // namespace noya