Extended Garner's algorithm (tools/extended_garner.hpp)

• View this file on GitHub
• View document part on GitHub
• Last update: 2025-02-10 01:30:13+09:00
• Include: #include "tools/extended_garner.hpp"
• Link: https://qiita.com/drken/items/ae02240cd1f8edfc86fd

template <typename Iterator, typename ModType>
std::optional<std::pair<long long, long long>> extended_garner(Iterator begin, Iterator end, ModType M);

template <typename M, typename Iterator>
std::optional<std::pair<M, M>> extended_garner(Iterator begin, Iterator end);

Given an array of the pairs $(r_i, m_i)$ with length $n$, it solves the modular equation system

\[\begin{align*} x \equiv r_i \pmod{m_i}, \forall i \in \{0, 1, \ldots, n - 1\} \end{align*}\]

If there is no solution, it returns std::nullopt. Otherwise, all the solutions can be written as the form $x \equiv y \pmod{z}$, using integers $y, z \, (0 \leq y < z = \mathrm{lcm}(m_i))$. It returns this $(y \pmod{m}, z \pmod{m})$ as a pair. If $n = 0$, it returns $(0, 1)$.

Constraints

• $m_i \geq 1$
• $M \geq 1$

Time Complexity

• $O(n^2 \log(\min(m_i)))$

References

• 中国剰余定理 (CRT) の解説と、それを用いる問題のまとめ - Qiita

License

• unknown

Author

• drken

Depends on

• std::abs(x) extended for my library (tools/abs.hpp)
• Extended Euclidean algorithm (tools/extgcd.hpp)
• Garner's algorithm (tools/garner.hpp)
• $x^{-1} \pmod{M}$ (tools/inv_mod.hpp)
• std::is_integral extended for tools::int128_t and tools::uint128_t (tools/is_integral.hpp)
• Minimum non-negative reminder (tools/mod.hpp)

Verified with

• tests/extended_garner.test.cpp

Code

#ifndef TOOLS_EXTENDED_GARNER_HPP
#define TOOLS_EXTENDED_GARNER_HPP

#include <optional>
#include <utility>
#include <cstdint>
#include <vector>
#include <cstddef>
#include "tools/mod.hpp"
#include "tools/garner.hpp"

// Source: https://qiita.com/drken/items/ae02240cd1f8edfc86fd
// License: unknown
// Author: drken

namespace tools {

  template <typename Iterator, typename ModType>
  ::std::optional<::std::pair<long long, long long>> extended_garner(const Iterator& begin, const Iterator& end, const ModType& mod) {
    ::std::vector<::std::pair<long long, long long>> v;
    for (auto it = begin; it != end; ++it) {
      v.emplace_back(::tools::mod(it->first, it->second), it->second);
    }

    for (::std::size_t i = 0; i < v.size(); ++i) {
      for (::std::size_t j = 0; j < i; ++j) {
        long long g = ::std::gcd(v[i].second, v[j].second);

        if ((v[i].first - v[j].first) % g != 0) return ::std::nullopt;

        v[i].second /= g;
        v[j].second /= g;

        long long gi = ::std::gcd(v[i].second, g);
        long long gj = g / gi;

        do {
          g = ::std::gcd(gi, gj);
          gi *= g;
          gj /= g;
        } while (g != 1);

        v[i].second *= gi;
        v[j].second *= gj;

        v[i].first %= v[i].second;
        v[j].first %= v[j].second;
      }
    }

    return ::std::optional<::std::pair<long long, long long>>(::tools::garner(v.begin(), v.end(), mod));
  }

  template <typename M, typename Iterator>
  ::std::optional<::std::pair<M, M>> extended_garner(const Iterator& begin, const Iterator& end) {
    const auto result = ::tools::extended_garner(begin, end, M::mod());
    if (!result) return ::std::nullopt;
    return ::std::make_optional<::std::pair<M, M>>(M::raw(result->first), M::raw(result->second));
  }
}

#endif

#line 1 "tools/extended_garner.hpp"



#include <optional>
#include <utility>
#include <cstdint>
#include <vector>
#include <cstddef>
#line 1 "tools/mod.hpp"



#include <cassert>
#include <type_traits>
#line 1 "tools/is_integral.hpp"



#line 5 "tools/is_integral.hpp"

namespace tools {
  template <typename T>
  struct is_integral : ::std::is_integral<T> {};

  template <typename T>
  inline constexpr bool is_integral_v = ::tools::is_integral<T>::value;
}


#line 7 "tools/mod.hpp"

namespace tools {

  template <typename M, typename N> requires (
    ::tools::is_integral_v<M> && !::std::is_same_v<::std::remove_cv_t<M>, bool> &&
    ::tools::is_integral_v<N> && !::std::is_same_v<::std::remove_cv_t<N>, bool>)
  constexpr ::std::common_type_t<M, N> mod(const M a, const N b) noexcept {
    assert(b != 0);

    using UM = ::std::make_unsigned_t<M>;
    using UN = ::std::make_unsigned_t<N>;
    const UM ua = a >= 0 ? a : static_cast<UM>(-(a + 1)) + 1;
    const UN ub = b >= 0 ? b : static_cast<UN>(-(b + 1)) + 1;
    auto r = ua % ub;
    if (a < 0 && r > 0) {
      r = ub - r;
    }
    return r;
  }
}


#line 1 "tools/garner.hpp"



#line 1 "tools/inv_mod.hpp"



#line 1 "tools/extgcd.hpp"



#include <tuple>
#line 7 "tools/extgcd.hpp"
#include <algorithm>
#line 1 "tools/abs.hpp"



namespace tools {
  constexpr float abs(const float x) {
    return x < 0 ? -x : x;
  }
  constexpr double abs(const double x) {
    return x < 0 ? -x : x;
  }
  constexpr long double abs(const long double x) {
    return x < 0 ? -x : x;
  }
  constexpr int abs(const int x) {
    return x < 0 ? -x : x;
  }
  constexpr long abs(const long x) {
    return x < 0 ? -x : x;
  }
  constexpr long long abs(const long long x) {
    return x < 0 ? -x : x;
  }
  constexpr unsigned int abs(const unsigned int x) {
    return x;
  }
  constexpr unsigned long abs(const unsigned long x) {
    return x;
  }
  constexpr unsigned long long abs(const unsigned long long x) {
    return x;
  }
}


#line 9 "tools/extgcd.hpp"

namespace tools {

  template <typename T>
  ::std::tuple<T, T, T> extgcd(T prev_r, T r) {
    const bool prev_r_is_neg = prev_r < T(0);
    const bool r_is_neg = r < T(0);

    if (prev_r_is_neg) prev_r = -prev_r;
    if (r_is_neg) r = -r;

    #ifndef NDEBUG
    const T a = prev_r;
    const T b = r;
    #endif

    T prev_s(1);
    T prev_t(0);
    T s(0);
    T t(1);
    while (r != T(0)) {
      const T q = prev_r / r;
      ::std::tie(prev_r, r) = ::std::make_pair(r, prev_r - q * r);
      ::std::tie(prev_s, s) = ::std::make_pair(s, prev_s - q * s);
      ::std::tie(prev_t, t) = ::std::make_pair(t, prev_t - q * t);
    }

    if (prev_r_is_neg) prev_s = -prev_s;
    if (r_is_neg) prev_t = -prev_t;

    assert(::tools::abs(prev_s) <= ::std::max(b / prev_r / T(2), T(1)));
    assert(::tools::abs(prev_t) <= ::std::max(a / prev_r / T(2), T(1)));
    return ::std::make_tuple(prev_s, prev_t, prev_r);
  }
}


#line 7 "tools/inv_mod.hpp"

namespace tools {

  template <typename T1, typename T2>
  constexpr T2 inv_mod(const T1 x, const T2 m) {
    const auto [x0, y0, gcd] = ::tools::extgcd(x, m);
    assert(gcd == 1);
    return ::tools::mod(x0, m);
  }
}


#line 9 "tools/garner.hpp"

// Source: https://qiita.com/drken/items/ae02240cd1f8edfc86fd
// License: unknown
// Author: drken

namespace tools {

  template <typename Iterator, typename ModType>
  ::std::pair<long long, long long> garner(const Iterator& begin, const Iterator& end, const ModType& mod) {
    ::std::vector<long long> b, m;
    for (auto it = begin; it != end; ++it) {
      b.push_back(::tools::mod(it->first, it->second));
      m.push_back(it->second);
    }

    auto lcm = 1LL;
    for (::std::size_t i = 0; i < b.size(); ++i) {
      (lcm *= m[i]) %= mod;
    }

    m.push_back(mod);
    ::std::vector<long long> coeffs(m.size(), 1);
    ::std::vector<long long> constants(m.size(), 0);
    for (::std::size_t k = 0; k < b.size(); ++k) {
      long long t = ::tools::mod((b[k] - constants[k]) * ::tools::inv_mod(coeffs[k], m[k]), m[k]);
      for (::std::size_t i = k + 1; i < m.size(); ++i) {
        (constants[i] += t * coeffs[i]) %= m[i];
        (coeffs[i] *= m[k]) %= m[i];
      }
    }

    return ::std::make_pair(constants.back(), lcm);
  }

  template <typename M, typename Iterator>
  ::std::pair<M, M> garner(const Iterator& begin, const Iterator& end) {
    const auto [y, z] = ::tools::garner(begin, end, M::mod());
    return ::std::make_pair(M::raw(y), M::raw(z));
  }
}


#line 11 "tools/extended_garner.hpp"

// Source: https://qiita.com/drken/items/ae02240cd1f8edfc86fd
// License: unknown
// Author: drken

namespace tools {

  template <typename Iterator, typename ModType>
  ::std::optional<::std::pair<long long, long long>> extended_garner(const Iterator& begin, const Iterator& end, const ModType& mod) {
    ::std::vector<::std::pair<long long, long long>> v;
    for (auto it = begin; it != end; ++it) {
      v.emplace_back(::tools::mod(it->first, it->second), it->second);
    }

    for (::std::size_t i = 0; i < v.size(); ++i) {
      for (::std::size_t j = 0; j < i; ++j) {
        long long g = ::std::gcd(v[i].second, v[j].second);

        if ((v[i].first - v[j].first) % g != 0) return ::std::nullopt;

        v[i].second /= g;
        v[j].second /= g;

        long long gi = ::std::gcd(v[i].second, g);
        long long gj = g / gi;

        do {
          g = ::std::gcd(gi, gj);
          gi *= g;
          gj /= g;
        } while (g != 1);

        v[i].second *= gi;
        v[j].second *= gj;

        v[i].first %= v[i].second;
        v[j].first %= v[j].second;
      }
    }

    return ::std::optional<::std::pair<long long, long long>>(::tools::garner(v.begin(), v.end(), mod));
  }

  template <typename M, typename Iterator>
  ::std::optional<::std::pair<M, M>> extended_garner(const Iterator& begin, const Iterator& end) {
    const auto result = ::tools::extended_garner(begin, end, M::mod());
    if (!result) return ::std::nullopt;
    return ::std::make_optional<::std::pair<M, M>>(M::raw(result->first), M::raw(result->second));
  }
}

Back to top page