floor_monoid_product.hpp¶
Universal Euclidean algorithm for noncommutative monoid products along a floor path.
\[
\displaystyle h_i=\left\lfloor\frac{ai+b}{m}\right\rfloor,\qquad W=Y^{h_0}\prod_{i=0}^{N-1}\left(XY^{h_{i+1}-h_i}\right)
\]
Implementation¶
#ifndef NOYA_FLOOR_MONOID_PRODUCT_HPP
#define NOYA_FLOOR_MONOID_PRODUCT_HPP 1
/// @complexity Time: O(log(max(n,m))) monoid compositions plus requested powers.
/// Space: O(log(max(n,m))) recursion stack.
#include <cassert>
#include <utility>
namespace noya {
namespace floor_monoid_product_internal {
template <class Monoid, class UInt>
typename Monoid::value_type power(typename Monoid::value_type value,
UInt exponent) {
typename Monoid::value_type result = Monoid::unit();
while (exponent > 0) {
if (exponent & 1) {
result = Monoid::op(result, value);
}
value = Monoid::op(value, value);
exponent >>= 1;
}
return result;
}
} // namespace floor_monoid_product_internal
/// @brief Universal Euclidean algorithm: evaluate a noncommutative monoid
/// product along the lattice path below y=(a*x+b)/modulus.
///
/// The returned word is Y^h(0) product_{i=0}^{n-1}
/// (X Y^(h(i+1)-h(i))), where h(i)=floor((a*i+b)/modulus). All integer
/// parameters must be nonnegative, modulus must be positive, and intermediate
/// products such as a*n+b must fit UInt.
template <class Monoid, class UInt>
typename Monoid::value_type
floor_monoid_product(typename Monoid::value_type horizontal,
typename Monoid::value_type vertical, UInt n, UInt a,
UInt b, UInt modulus) {
assert(modulus > 0);
using value_type = typename Monoid::value_type;
using floor_monoid_product_internal::power;
UInt remaining_vertical = (a * n + b) / modulus;
value_type prefix = Monoid::unit();
value_type suffix = Monoid::unit();
while (true) {
UInt slope_quotient = a / modulus;
UInt intercept_quotient = b / modulus;
a %= modulus;
b %= modulus;
horizontal =
Monoid::op(horizontal, power<Monoid>(vertical, slope_quotient));
prefix = Monoid::op(prefix, power<Monoid>(vertical, intercept_quotient));
remaining_vertical -= slope_quotient * n + intercept_quotient;
if (remaining_vertical == 0) {
break;
}
assert(a > 0);
UInt split = (modulus * remaining_vertical - b - 1) / a + 1;
suffix = Monoid::op(
vertical, Monoid::op(power<Monoid>(horizontal, n - split), suffix));
b = modulus - b - 1 + a;
n = remaining_vertical - 1;
remaining_vertical = split;
std::swap(modulus, a);
std::swap(horizontal, vertical);
}
return Monoid::op(Monoid::op(prefix, power<Monoid>(horizontal, n)), suffix);
}
} // namespace noya
#endif // NOYA_FLOOR_MONOID_PRODUCT_HPP
#include <cassert>
#include <utility>
/// @complexity Time: O(log(max(n,m))) monoid compositions plus requested powers.
/// Space: O(log(max(n,m))) recursion stack.
namespace noya {
namespace floor_monoid_product_internal {
template <class Monoid, class UInt>
typename Monoid::value_type power(typename Monoid::value_type value,
UInt exponent) {
typename Monoid::value_type result = Monoid::unit();
while (exponent > 0) {
if (exponent & 1) {
result = Monoid::op(result, value);
}
value = Monoid::op(value, value);
exponent >>= 1;
}
return result;
}
} // namespace floor_monoid_product_internal
/// @brief Universal Euclidean algorithm: evaluate a noncommutative monoid
/// product along the lattice path below y=(a*x+b)/modulus.
///
/// The returned word is Y^h(0) product_{i=0}^{n-1}
/// (X Y^(h(i+1)-h(i))), where h(i)=floor((a*i+b)/modulus). All integer
/// parameters must be nonnegative, modulus must be positive, and intermediate
/// products such as a*n+b must fit UInt.
template <class Monoid, class UInt>
typename Monoid::value_type
floor_monoid_product(typename Monoid::value_type horizontal,
typename Monoid::value_type vertical, UInt n, UInt a,
UInt b, UInt modulus) {
assert(modulus > 0);
using value_type = typename Monoid::value_type;
using floor_monoid_product_internal::power;
UInt remaining_vertical = (a * n + b) / modulus;
value_type prefix = Monoid::unit();
value_type suffix = Monoid::unit();
while (true) {
UInt slope_quotient = a / modulus;
UInt intercept_quotient = b / modulus;
a %= modulus;
b %= modulus;
horizontal =
Monoid::op(horizontal, power<Monoid>(vertical, slope_quotient));
prefix = Monoid::op(prefix, power<Monoid>(vertical, intercept_quotient));
remaining_vertical -= slope_quotient * n + intercept_quotient;
if (remaining_vertical == 0) {
break;
}
assert(a > 0);
UInt split = (modulus * remaining_vertical - b - 1) / a + 1;
suffix = Monoid::op(
vertical, Monoid::op(power<Monoid>(horizontal, n - split), suffix));
b = modulus - b - 1 + a;
n = remaining_vertical - 1;
remaining_vertical = split;
std::swap(modulus, a);
std::swap(horizontal, vertical);
}
return Monoid::op(Monoid::op(prefix, power<Monoid>(horizontal, n)), suffix);
}
} // namespace noya