tetration_mod.hpp

SECTIONMath INCLUDEnoya/tetration_mod.hpp

Compute a power tower of height copies of base modulo modulus, with 0^0 = 1. Recursing through Euler-totient moduli keeps the exponent small; an extra "large" bit distinguishes a reduced exponent from a genuinely small one when the base is not coprime to the modulus.

Verified by tetration_mod.

\[ \displaystyle a\uparrow\uparrow h \bmod m \]

Implementation

View on GitHub

SourceExpanded

#ifndef NOYA_TETRATION_MOD_HPP
#define NOYA_TETRATION_MOD_HPP 1

/// @complexity Time: O(sqrt(modulus) log modulus).  Space: O(log modulus).

#include <cassert>
#include <cstdint>

namespace noya {

namespace tetration_mod_detail {

inline std::uint64_t totient(std::uint64_t value) {
  std::uint64_t result = value;
  for (std::uint64_t factor = 2; factor * factor <= value; factor++) {
    if (value % factor == 0) {
      result = result / factor * (factor - 1);
      do {
        value /= factor;
      } while (value % factor == 0);
    }
  }
  if (value != 1) {
    result = result / value * (value - 1);
  }
  return result;
}

inline std::uint64_t power(std::uint64_t base, std::uint64_t exponent,
                           std::uint64_t modulus, bool &large) {
  if (base == 1 || exponent == 0) {
    return 1;
  }
  std::uint64_t result = 1;
  large = large || base >= modulus;
  base %= modulus;
  while (true) {
    if (exponent & 1) {
      unsigned __int128 product =
          static_cast<unsigned __int128>(result) * base;
      if (product >= modulus) {
        large = true;
      }
      result = std::uint64_t(product % modulus);
    }
    exponent >>= 1;
    if (exponent == 0) {
      break;
    }
    unsigned __int128 square =
        static_cast<unsigned __int128>(base) * base;
    if (square >= modulus) {
      large = true;
    }
    base = std::uint64_t(square % modulus);
  }
  return result;
}

inline std::uint64_t tower(std::uint64_t base, std::uint64_t height,
                           std::uint64_t modulus, bool &large) {
  if (base == 0) {
    return (~height) & 1ULL;
  }
  if (modulus == 1) {
    large = true;
    return 1;
  }
  if (base == 1 || height == 0) {
    return 1;
  }
  if (height == 1) {
    large = large || base >= modulus;
    return base % modulus + std::uint64_t(large) * modulus;
  }
  std::uint64_t exponent =
      tower(base, height - 1, totient(modulus), large);
  std::uint64_t result = power(base, exponent, modulus, large);
  return result + std::uint64_t(large) * modulus;
}

} // namespace tetration_mod_detail

/// @brief Compute a power tower of `height` copies of `base` modulo
/// `modulus`, with `0^0 = 1`.  Recursing through Euler-totient moduli keeps
/// the exponent small; an extra "large" bit distinguishes a reduced exponent
/// from a genuinely small one when the base is not coprime to the modulus.
inline std::uint64_t tetration_mod(std::uint64_t base, std::uint64_t height,
                                   std::uint64_t modulus) {
  assert(modulus >= 1);
  bool large = false;
  return tetration_mod_detail::tower(base, height, modulus, large) % modulus;
}

} // namespace noya

#endif // NOYA_TETRATION_MOD_HPP

#include <cassert>
#include <cstdint>

/// @complexity Time: O(sqrt(modulus) log modulus).  Space: O(log modulus).

namespace noya {

namespace tetration_mod_detail {

inline std::uint64_t totient(std::uint64_t value) {
  std::uint64_t result = value;
  for (std::uint64_t factor = 2; factor * factor <= value; factor++) {
    if (value % factor == 0) {
      result = result / factor * (factor - 1);
      do {
        value /= factor;
      } while (value % factor == 0);
    }
  }
  if (value != 1) {
    result = result / value * (value - 1);
  }
  return result;
}

inline std::uint64_t power(std::uint64_t base, std::uint64_t exponent,
                           std::uint64_t modulus, bool &large) {
  if (base == 1 || exponent == 0) {
    return 1;
  }
  std::uint64_t result = 1;
  large = large || base >= modulus;
  base %= modulus;
  while (true) {
    if (exponent & 1) {
      unsigned __int128 product =
          static_cast<unsigned __int128>(result) * base;
      if (product >= modulus) {
        large = true;
      }
      result = std::uint64_t(product % modulus);
    }
    exponent >>= 1;
    if (exponent == 0) {
      break;
    }
    unsigned __int128 square =
        static_cast<unsigned __int128>(base) * base;
    if (square >= modulus) {
      large = true;
    }
    base = std::uint64_t(square % modulus);
  }
  return result;
}

inline std::uint64_t tower(std::uint64_t base, std::uint64_t height,
                           std::uint64_t modulus, bool &large) {
  if (base == 0) {
    return (~height) & 1ULL;
  }
  if (modulus == 1) {
    large = true;
    return 1;
  }
  if (base == 1 || height == 0) {
    return 1;
  }
  if (height == 1) {
    large = large || base >= modulus;
    return base % modulus + std::uint64_t(large) * modulus;
  }
  std::uint64_t exponent =
      tower(base, height - 1, totient(modulus), large);
  std::uint64_t result = power(base, exponent, modulus, large);
  return result + std::uint64_t(large) * modulus;
}

} // namespace tetration_mod_detail

/// @brief Compute a power tower of `height` copies of `base` modulo
/// `modulus`, with `0^0 = 1`.  Recursing through Euler-totient moduli keeps
/// the exponent small; an extra "large" bit distinguishes a reduced exponent
/// from a genuinely small one when the base is not coprime to the modulus.
inline std::uint64_t tetration_mod(std::uint64_t base, std::uint64_t height,
                                   std::uint64_t modulus) {
  assert(modulus >= 1);
  bool large = false;
  return tetration_mod_detail::tower(base, height, modulus, large) % modulus;
}

} // namespace noya