min_of_mod_linear.hpp¶
Return an index attaining the minimum of (multiplier * x + addition) mod modulus for 0 <= x < count, together with that minimum. Continued-fraction neighbors describe every index at which the prefix minimum decreases; those indices form arithmetic runs, so only one run per Euclidean-algorithm step is inspected.
Verified by min_of_mod_of_linear.
\[
\displaystyle \min_{0\le x<n}(ax+b)\bmod m
\]
Implementation¶
#ifndef NOYA_MIN_OF_MOD_LINEAR_HPP
#define NOYA_MIN_OF_MOD_LINEAR_HPP 1
/// @complexity Time: O(log modulus). Space: O(log modulus).
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <utility>
#include <vector>
namespace noya {
namespace min_of_mod_linear_detail {
inline std::pair<std::vector<std::int64_t>, std::vector<std::int64_t>>
record_segments(std::int64_t multiplier, std::int64_t addition,
std::int64_t modulus) {
std::vector<std::int64_t> boundaries{0};
std::vector<std::int64_t> steps;
std::int64_t divisor = std::gcd(multiplier, modulus);
multiplier /= divisor;
addition /= divisor;
modulus /= divisor;
std::int64_t left_numerator = 0;
std::int64_t left_denominator = 1;
std::int64_t right_numerator = 1;
std::int64_t right_denominator = 1;
std::int64_t left_determinant = modulus - multiplier;
std::int64_t right_determinant = multiplier;
std::int64_t index = 0;
std::int64_t value = addition;
while (value != 0) {
std::int64_t quotient = right_determinant / left_determinant;
right_determinant %= left_determinant;
if (right_determinant == 0) {
--quotient;
right_determinant = left_determinant;
}
right_numerator += quotient * left_numerator;
right_denominator += quotient * left_denominator;
while (true) {
quotient = std::max<std::int64_t>(
0, (left_determinant - value + right_determinant - 1) /
right_determinant);
if (left_determinant - quotient * right_determinant <= 0) {
break;
}
left_determinant -= quotient * right_determinant;
left_numerator += quotient * right_numerator;
left_denominator += quotient * right_denominator;
quotient = value / left_determinant;
value -= quotient * left_determinant;
index += left_denominator * quotient;
boundaries.push_back(index);
steps.push_back(left_denominator);
}
quotient = left_determinant / right_determinant;
left_determinant -= quotient * right_determinant;
left_numerator += quotient * right_numerator;
left_denominator += quotient * right_denominator;
}
return {boundaries, steps};
}
} // namespace min_of_mod_linear_detail
/// @brief Return an index attaining the minimum of
/// `(multiplier * x + addition) mod modulus` for `0 <= x < count`, together
/// with that minimum. Continued-fraction neighbors describe every index at
/// which the prefix minimum decreases; those indices form arithmetic runs, so
/// only one run per Euclidean-algorithm step is inspected.
inline std::pair<std::uint64_t, std::uint64_t>
min_of_mod_linear_argument(std::uint64_t count, std::uint64_t modulus,
std::uint64_t multiplier,
std::uint64_t addition) {
assert(count >= 1 && modulus >= 1);
assert(multiplier < modulus && addition < modulus);
assert(count <= 1'000'000'000ULL && modulus <= 1'000'000'000ULL);
auto [boundaries, steps] = min_of_mod_linear_detail::record_segments(
std::int64_t(multiplier), std::int64_t(addition),
std::int64_t(modulus));
std::int64_t index = 0;
for (std::size_t run = 0; run + 1 < boundaries.size(); run++) {
std::int64_t left = boundaries[run];
std::int64_t right = boundaries[run + 1];
if (std::uint64_t(right) < count) {
index = right;
continue;
}
index = left + std::int64_t((count - 1 - std::uint64_t(left)) /
std::uint64_t(steps[run])) *
steps[run];
break;
}
std::uint64_t value = std::uint64_t(
(static_cast<unsigned __int128>(multiplier) * std::uint64_t(index) +
addition) %
modulus);
return {std::uint64_t(index), value};
}
/// @brief Return the minimum residue of a linear function on
/// `0 <= x < count` using the continued-fraction prefix-minimum decomposition.
inline std::uint64_t min_of_mod_linear(std::uint64_t count,
std::uint64_t modulus,
std::uint64_t multiplier,
std::uint64_t addition) {
return min_of_mod_linear_argument(count, modulus, multiplier, addition)
.second;
}
} // namespace noya
#endif // NOYA_MIN_OF_MOD_LINEAR_HPP
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <utility>
#include <vector>
/// @complexity Time: O(log modulus). Space: O(log modulus).
namespace noya {
namespace min_of_mod_linear_detail {
inline std::pair<std::vector<std::int64_t>, std::vector<std::int64_t>>
record_segments(std::int64_t multiplier, std::int64_t addition,
std::int64_t modulus) {
std::vector<std::int64_t> boundaries{0};
std::vector<std::int64_t> steps;
std::int64_t divisor = std::gcd(multiplier, modulus);
multiplier /= divisor;
addition /= divisor;
modulus /= divisor;
std::int64_t left_numerator = 0;
std::int64_t left_denominator = 1;
std::int64_t right_numerator = 1;
std::int64_t right_denominator = 1;
std::int64_t left_determinant = modulus - multiplier;
std::int64_t right_determinant = multiplier;
std::int64_t index = 0;
std::int64_t value = addition;
while (value != 0) {
std::int64_t quotient = right_determinant / left_determinant;
right_determinant %= left_determinant;
if (right_determinant == 0) {
--quotient;
right_determinant = left_determinant;
}
right_numerator += quotient * left_numerator;
right_denominator += quotient * left_denominator;
while (true) {
quotient = std::max<std::int64_t>(
0, (left_determinant - value + right_determinant - 1) /
right_determinant);
if (left_determinant - quotient * right_determinant <= 0) {
break;
}
left_determinant -= quotient * right_determinant;
left_numerator += quotient * right_numerator;
left_denominator += quotient * right_denominator;
quotient = value / left_determinant;
value -= quotient * left_determinant;
index += left_denominator * quotient;
boundaries.push_back(index);
steps.push_back(left_denominator);
}
quotient = left_determinant / right_determinant;
left_determinant -= quotient * right_determinant;
left_numerator += quotient * right_numerator;
left_denominator += quotient * right_denominator;
}
return {boundaries, steps};
}
} // namespace min_of_mod_linear_detail
/// @brief Return an index attaining the minimum of
/// `(multiplier * x + addition) mod modulus` for `0 <= x < count`, together
/// with that minimum. Continued-fraction neighbors describe every index at
/// which the prefix minimum decreases; those indices form arithmetic runs, so
/// only one run per Euclidean-algorithm step is inspected.
inline std::pair<std::uint64_t, std::uint64_t>
min_of_mod_linear_argument(std::uint64_t count, std::uint64_t modulus,
std::uint64_t multiplier,
std::uint64_t addition) {
assert(count >= 1 && modulus >= 1);
assert(multiplier < modulus && addition < modulus);
assert(count <= 1'000'000'000ULL && modulus <= 1'000'000'000ULL);
auto [boundaries, steps] = min_of_mod_linear_detail::record_segments(
std::int64_t(multiplier), std::int64_t(addition),
std::int64_t(modulus));
std::int64_t index = 0;
for (std::size_t run = 0; run + 1 < boundaries.size(); run++) {
std::int64_t left = boundaries[run];
std::int64_t right = boundaries[run + 1];
if (std::uint64_t(right) < count) {
index = right;
continue;
}
index = left + std::int64_t((count - 1 - std::uint64_t(left)) /
std::uint64_t(steps[run])) *
steps[run];
break;
}
std::uint64_t value = std::uint64_t(
(static_cast<unsigned __int128>(multiplier) * std::uint64_t(index) +
addition) %
modulus);
return {std::uint64_t(index), value};
}
/// @brief Return the minimum residue of a linear function on
/// `0 <= x < count` using the continued-fraction prefix-minimum decomposition.
inline std::uint64_t min_of_mod_linear(std::uint64_t count,
std::uint64_t modulus,
std::uint64_t multiplier,
std::uint64_t addition) {
return min_of_mod_linear_argument(count, modulus, multiplier, addition)
.second;
}
} // namespace noya