floor_monoid_product.hpp¶
沿 \(\left\lfloor\frac{ai+b}{m}\right\rfloor\) 描述的格路径计算非交换幺半群乘积；用于欧几里得递归可压缩的序列积。
\[
\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)
\]
Complexity: Time: \(O(\log (a + \mathrm{mod} + 1))\) Euclidean steps, each using \(O(\log (n + h + 2))\) monoid compositions, where \(h = \lfloor (a \cdot n + b) / \mathrm{mod}\rfloor\). Space: \(O(1)\) monoid values and integer state, excluding Monoid::op storage.
Source: maspy
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(log(a + mod + 1)) Euclidean steps, each using
/// O(log(n + h + 2)) monoid compositions, where h=floor((a*n+b)/mod).
/// Space: O(1) monoid values and integer state, excluding Monoid::op storage.
#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 val, UInt exp) {
typename Monoid::value_type res = Monoid::unit();
while (exp > 0) {
if (exp & 1) {
res = Monoid::op(res, val);
}
exp >>= 1;
if (exp > 0) val = Monoid::op(val, val);
}
return res;
}
} // namespace floor_monoid_product_internal
/// @brief Universal Euclidean algorithm: evaluate a noncommutative monoid
/// product along the lattice path below y=(a*x+b)/mod.
///
/// The returned word is ver^h(0) product_{i=0}^{n-1}
/// (hor ver^(h(i+1)-h(i))), where h(i)=floor((a*i+b)/mod). All integer
/// parameters must be nonnegative, mod must be positive, and a*n+b must fit UInt.
template <class Monoid, class UInt>
typename Monoid::value_type
floor_monoid_product(typename Monoid::value_type hor,
typename Monoid::value_type ver, UInt n, UInt a, UInt b,
UInt mod) {
assert(mod > 0);
using value_type = typename Monoid::value_type;
using floor_monoid_product_internal::power;
UInt rv = (a * n + b) / mod;
value_type pre = Monoid::unit();
value_type suf = Monoid::unit();
while (true) {
UInt iq = b / mod;
pre = Monoid::op(pre, power<Monoid>(ver, iq));
if (n == 0) break;
UInt sq = a / mod;
a %= mod;
b %= mod;
hor = Monoid::op(hor, power<Monoid>(ver, sq));
rv -= sq * n + iq;
if (rv == 0) {
break;
}
assert(a > 0);
UInt spl = (mod * rv - b - 1) / a + 1;
suf = Monoid::op(ver, Monoid::op(power<Monoid>(hor, n - spl), suf));
// Remove one vertical step before swapping axes, avoiding b+a overflow.
pre = Monoid::op(pre, hor);
b = mod - b - 1;
n = rv - 1;
rv = spl - 1;
std::swap(mod, a);
std::swap(hor, ver);
}
return Monoid::op(Monoid::op(pre, power<Monoid>(hor, n)), suf);
}
} // namespace noya
// NOYA_FLOOR_MONOID_PRODUCT_HPP
#ifndef NOYA_FLOOR_MONOID_PRODUCT_HPP
#define NOYA_FLOOR_MONOID_PRODUCT_HPP 1
/// @complexity Time: O(log(a + mod + 1)) Euclidean steps, each using
/// O(log(n + h + 2)) monoid compositions, where h=floor((a*n+b)/mod).
/// Space: O(1) monoid values and integer state, excluding Monoid::op storage.
#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 val, UInt exp) {
typename Monoid::value_type res = Monoid::unit();
while (exp > 0) {
if (exp & 1) {
res = Monoid::op(res, val);
}
exp >>= 1;
if (exp > 0) val = Monoid::op(val, val);
}
return res;
}
} // namespace floor_monoid_product_internal
/// @brief Universal Euclidean algorithm: evaluate a noncommutative monoid
/// product along the lattice path below y=(a*x+b)/mod.
///
/// The returned word is ver^h(0) product_{i=0}^{n-1}
/// (hor ver^(h(i+1)-h(i))), where h(i)=floor((a*i+b)/mod). All integer
/// parameters must be nonnegative, mod must be positive, and a*n+b must fit UInt.
template <class Monoid, class UInt>
typename Monoid::value_type
floor_monoid_product(typename Monoid::value_type hor,
typename Monoid::value_type ver, UInt n, UInt a, UInt b,
UInt mod) {
assert(mod > 0);
using value_type = typename Monoid::value_type;
using floor_monoid_product_internal::power;
UInt rv = (a * n + b) / mod;
value_type pre = Monoid::unit();
value_type suf = Monoid::unit();
while (true) {
UInt iq = b / mod;
pre = Monoid::op(pre, power<Monoid>(ver, iq));
if (n == 0) break;
UInt sq = a / mod;
a %= mod;
b %= mod;
hor = Monoid::op(hor, power<Monoid>(ver, sq));
rv -= sq * n + iq;
if (rv == 0) {
break;
}
assert(a > 0);
UInt spl = (mod * rv - b - 1) / a + 1;
suf = Monoid::op(ver, Monoid::op(power<Monoid>(hor, n - spl), suf));
// Remove one vertical step before swapping axes, avoiding b+a overflow.
pre = Monoid::op(pre, hor);
b = mod - b - 1;
n = rv - 1;
rv = spl - 1;
std::swap(mod, a);
std::swap(hor, ver);
}
return Monoid::op(Monoid::op(pre, power<Monoid>(hor, n)), suf);
}
} // namespace noya
#endif // NOYA_FLOOR_MONOID_PRODUCT_HPP
/// @complexity Time: O(log(a + mod + 1)) Euclidean steps, each using
/// O(log(n + h + 2)) monoid compositions, where h=floor((a*n+b)/mod).
/// Space: O(1) monoid values and integer state, excluding Monoid::op storage.
#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 val, UInt exp) {
typename Monoid::value_type res = Monoid::unit();
while (exp > 0) {
if (exp & 1) {
res = Monoid::op(res, val);
}
exp >>= 1;
if (exp > 0) val = Monoid::op(val, val);
}
return res;
}
} // namespace floor_monoid_product_internal
/// @brief Universal Euclidean algorithm: evaluate a noncommutative monoid
/// product along the lattice path below y=(a*x+b)/mod.
///
/// The returned word is ver^h(0) product_{i=0}^{n-1}
/// (hor ver^(h(i+1)-h(i))), where h(i)=floor((a*i+b)/mod). All integer
/// parameters must be nonnegative, mod must be positive, and a*n+b must fit UInt.
template <class Monoid, class UInt>
typename Monoid::value_type
floor_monoid_product(typename Monoid::value_type hor,
typename Monoid::value_type ver, UInt n, UInt a, UInt b,
UInt mod) {
assert(mod > 0);
using value_type = typename Monoid::value_type;
using floor_monoid_product_internal::power;
UInt rv = (a * n + b) / mod;
value_type pre = Monoid::unit();
value_type suf = Monoid::unit();
while (true) {
UInt iq = b / mod;
pre = Monoid::op(pre, power<Monoid>(ver, iq));
if (n == 0) break;
UInt sq = a / mod;
a %= mod;
b %= mod;
hor = Monoid::op(hor, power<Monoid>(ver, sq));
rv -= sq * n + iq;
if (rv == 0) {
break;
}
assert(a > 0);
UInt spl = (mod * rv - b - 1) / a + 1;
suf = Monoid::op(ver, Monoid::op(power<Monoid>(hor, n - spl), suf));
// Remove one vertical step before swapping axes, avoiding b+a overflow.
pre = Monoid::op(pre, hor);
b = mod - b - 1;
n = rv - 1;
rv = spl - 1;
std::swap(mod, a);
std::swap(hor, ver);
}
return Monoid::op(Monoid::op(pre, power<Monoid>(hor, n)), suf);
}
} // namespace noya
// NOYA_FLOOR_MONOID_PRODUCT_HPP