decimal_power_mod.hpp

SECTIONMath INCLUDEnoya/math/decimal_power_mod.hpp

指数以超长十进制字符串给出时计算模幂，不要求底数与模数互质。

\[ \displaystyle a^e\equiv r\;(\bmod m) \]

Complexity: Time: \(O(d)\) modular multiplications for \(d\) decimal exponent digits. Space: \(O(1)\) beyond the exponent string.

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

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

/// @complexity Time: O(d) modular multiplications for d decimal exponent digits.
/// Space: O(1) beyond the exponent string.

#include <array>
#include <cassert>
#include <cstdint>
#include <string_view>

namespace noya {

/// @brief Compute bas^exp modulo an arbitrary positive 64-bit modulus
/// from a decimal exponent of any length, without coprimality assumptions.
inline std::uint64_t decimal_power_mod(std::uint64_t bas, std::string_view exp,
                                       std::uint64_t mod) {
  assert(mod >= 1);
  assert(!exp.empty());
  for ([[maybe_unused]] char dig : exp) {
    assert('0' <= dig && dig <= '9');
  }
  using u128 = unsigned __int128;
  auto mul = [&](std::uint64_t a, std::uint64_t b) {
    return std::uint64_t(u128(a) * b % mod);
  };
  bas %= mod;
  auto first = exp.find_first_not_of('0');
  // Preserve 0^0 = 1 (mod mod), including an all-zero decimal string.
  if (first == std::string_view::npos) {
    return 1 % mod;
  }
  exp.remove_prefix(first);
  if (bas <= 1) {
    return bas;
  }
  if (bas == mod - 1) {
    return ((exp.back() - '0') & 1) ? bas : 1;
  }
  // Every decimal string of at most 19 digits fits uint64_t.
  if (exp.size() <= 19) {
    std::uint64_t exponent = 0;
    for (char dig : exp) {
      exponent = exponent * 10 + dig - '0';
    }
    std::uint64_t res = 1;
    while (exponent > 0) {
      if (exponent & 1) {
        res = mul(res, bas);
      }
      exponent >>= 1;
      if (exponent) {
        bas = mul(bas, bas);
      }
    }
    return res;
  }
  std::array<std::uint64_t, 10> powers{1 % mod, bas};
  for (int i = 2; i < 10; i++) {
    powers[i] = mul(powers[i - 1], bas);
  }
  std::uint64_t res = powers[exp.front() - '0'];
  exp.remove_prefix(1);
  if (res == 0) {
    return 0;
  }
  for (char dig : exp) {
    // Append digit d: a^(10*e+d) = (a^e)^10 * a^d.
    auto x2 = mul(res, res);
    auto x4 = mul(x2, x2);
    auto x8 = mul(x4, x4);
    res = mul(mul(x8, x2), powers[dig - '0']);
    // Further decimal digits keep a zero residue zero.
    if (res == 0) {
      return 0;
    }
  }
  return res;
}

} // namespace noya

 // NOYA_DECIMAL_POWER_MOD_HPP

#ifndef NOYA_DECIMAL_POWER_MOD_HPP
#define NOYA_DECIMAL_POWER_MOD_HPP 1

/// @complexity Time: O(d) modular multiplications for d decimal exponent digits.
/// Space: O(1) beyond the exponent string.

#include <array>
#include <cassert>
#include <cstdint>
#include <string_view>

namespace noya {

/// @brief Compute bas^exp modulo an arbitrary positive 64-bit modulus
/// from a decimal exponent of any length, without coprimality assumptions.
inline std::uint64_t decimal_power_mod(std::uint64_t bas, std::string_view exp,
                                       std::uint64_t mod) {
  assert(mod >= 1);
  assert(!exp.empty());
  for ([[maybe_unused]] char dig : exp) {
    assert('0' <= dig && dig <= '9');
  }
  using u128 = unsigned __int128;
  auto mul = [&](std::uint64_t a, std::uint64_t b) {
    return std::uint64_t(u128(a) * b % mod);
  };
  bas %= mod;
  auto first = exp.find_first_not_of('0');
  // Preserve 0^0 = 1 (mod mod), including an all-zero decimal string.
  if (first == std::string_view::npos) {
    return 1 % mod;
  }
  exp.remove_prefix(first);
  if (bas <= 1) {
    return bas;
  }
  if (bas == mod - 1) {
    return ((exp.back() - '0') & 1) ? bas : 1;
  }
  // Every decimal string of at most 19 digits fits uint64_t.
  if (exp.size() <= 19) {
    std::uint64_t exponent = 0;
    for (char dig : exp) {
      exponent = exponent * 10 + dig - '0';
    }
    std::uint64_t res = 1;
    while (exponent > 0) {
      if (exponent & 1) {
        res = mul(res, bas);
      }
      exponent >>= 1;
      if (exponent) {
        bas = mul(bas, bas);
      }
    }
    return res;
  }
  std::array<std::uint64_t, 10> powers{1 % mod, bas};
  for (int i = 2; i < 10; i++) {
    powers[i] = mul(powers[i - 1], bas);
  }
  std::uint64_t res = powers[exp.front() - '0'];
  exp.remove_prefix(1);
  if (res == 0) {
    return 0;
  }
  for (char dig : exp) {
    // Append digit d: a^(10*e+d) = (a^e)^10 * a^d.
    auto x2 = mul(res, res);
    auto x4 = mul(x2, x2);
    auto x8 = mul(x4, x4);
    res = mul(mul(x8, x2), powers[dig - '0']);
    // Further decimal digits keep a zero residue zero.
    if (res == 0) {
      return 0;
    }
  }
  return res;
}

} // namespace noya

#endif // NOYA_DECIMAL_POWER_MOD_HPP

/// @complexity Time: O(d) modular multiplications for d decimal exponent digits.
/// Space: O(1) beyond the exponent string.

#include <array>
#include <cassert>
#include <cstdint>
#include <string_view>

namespace noya {

/// @brief Compute bas^exp modulo an arbitrary positive 64-bit modulus
/// from a decimal exponent of any length, without coprimality assumptions.
inline std::uint64_t decimal_power_mod(std::uint64_t bas, std::string_view exp,
                                       std::uint64_t mod) {
  assert(mod >= 1);
  assert(!exp.empty());
  for ([[maybe_unused]] char dig : exp) {
    assert('0' <= dig && dig <= '9');
  }
  using u128 = unsigned __int128;
  auto mul = [&](std::uint64_t a, std::uint64_t b) {
    return std::uint64_t(u128(a) * b % mod);
  };
  bas %= mod;
  auto first = exp.find_first_not_of('0');
  // Preserve 0^0 = 1 (mod mod), including an all-zero decimal string.
  if (first == std::string_view::npos) {
    return 1 % mod;
  }
  exp.remove_prefix(first);
  if (bas <= 1) {
    return bas;
  }
  if (bas == mod - 1) {
    return ((exp.back() - '0') & 1) ? bas : 1;
  }
  // Every decimal string of at most 19 digits fits uint64_t.
  if (exp.size() <= 19) {
    std::uint64_t exponent = 0;
    for (char dig : exp) {
      exponent = exponent * 10 + dig - '0';
    }
    std::uint64_t res = 1;
    while (exponent > 0) {
      if (exponent & 1) {
        res = mul(res, bas);
      }
      exponent >>= 1;
      if (exponent) {
        bas = mul(bas, bas);
      }
    }
    return res;
  }
  std::array<std::uint64_t, 10> powers{1 % mod, bas};
  for (int i = 2; i < 10; i++) {
    powers[i] = mul(powers[i - 1], bas);
  }
  std::uint64_t res = powers[exp.front() - '0'];
  exp.remove_prefix(1);
  if (res == 0) {
    return 0;
  }
  for (char dig : exp) {
    // Append digit d: a^(10*e+d) = (a^e)^10 * a^d.
    auto x2 = mul(res, res);
    auto x4 = mul(x2, x2);
    auto x8 = mul(x4, x4);
    res = mul(mul(x8, x2), powers[dig - '0']);
    // Further decimal digits keep a zero residue zero.
    if (res == 0) {
      return 0;
    }
  }
  return res;
}

} // namespace noya

 // NOYA_DECIMAL_POWER_MOD_HPP