math/congruence_equations.cpp

• View this file on GitHub
• Last update: 2022-09-27 19:53:51+09:00

Depends on

• math/inv_gcd.cpp
• math/rem_euclid.cpp
• utility/rep.cpp

Code

#pragma once
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
#include "../utility/rep.cpp"
#include "inv_gcd.cpp"
#include "rem_euclid.cpp"

template <class T, class Result = T>
std::optional<std::pair<Result, Result>> congruence_equations(const std::vector<T>& m, const std::vector<T>& r) {
    using U = std::make_signed_t<Result>;
    assert(m.size() == r.size());
    const int n = m.size();
    U r0 = 0, m0 = 1;
    for (const int i : rep(0, n)) {
        assert(m[i] > 0);
        U r1 = rem_euclid(r[i], m[i]), m1 = m[i];
        if (m0 < m1) {
            std::swap(r0, r1);
            std::swap(m0, m1);
        }
        if (m0 % m1 == 0) {
            if (r0 % m1 != r1) return std::nullopt;
            continue;
        }
        const auto [g, inv_m] = inv_gcd(m0, m1);
        if ((r1 - r0) % g != 0) return std::nullopt;
        const U u1 = m1 / g;
        const U x = (r1 - r0) / g % u1 * inv_m % u1;
        r0 += x * m0;
        m0 *= u1;
        if (r0 < 0) r0 += m0;
    }
    return std::pair<Result, Result>(r0, m0);
}

#line 2 "math/congruence_equations.cpp"
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>
#line 2 "utility/rep.cpp"
#include <algorithm>

class Range {
    struct Iter {
        int itr;
        constexpr Iter(const int pos) noexcept : itr(pos) {}
        constexpr void operator++() noexcept { ++itr; }
        constexpr bool operator!=(const Iter& other) const noexcept { return itr != other.itr; }
        constexpr int operator*() const noexcept { return itr; }
    };
    const Iter first, last;

  public:
    explicit constexpr Range(const int first, const int last) noexcept : first(first), last(std::max(first, last)) {}
    constexpr Iter begin() const noexcept { return first; }
    constexpr Iter end() const noexcept { return last; }
};

constexpr Range rep(const int l, const int r) noexcept { return Range(l, r); }
constexpr Range rep(const int n) noexcept { return Range(0, n); }
#line 2 "math/rem_euclid.cpp"
#include <cassert>

template <class T> constexpr T rem_euclid(T value, const T& mod) {
    assert(mod > 0);
    return (value %= mod) >= 0 ? value : value + mod;
}
#line 5 "math/inv_gcd.cpp"

template <class T> constexpr std::pair<T, T> inv_gcd(const T& a, const T& b) {
    using U = std::make_signed_t<T>;
    U t = rem_euclid(a, b);
    if (t == 0) return {b, 0};
    U s = b, m0 = 0, m1 = 1;
    while (t != 0) {
        const U u = s / t;
        s -= t * u;
        m0 -= m1 * u;
        U tmp = s;
        s = t;
        t = tmp;
        tmp = m0;
        m0 = m1;
        m1 = tmp;
    }
    if (m0 < 0) m0 += b / s;
    return {(T)s, (T)m0};
}
#line 9 "math/congruence_equations.cpp"

template <class T, class Result = T>
std::optional<std::pair<Result, Result>> congruence_equations(const std::vector<T>& m, const std::vector<T>& r) {
    using U = std::make_signed_t<Result>;
    assert(m.size() == r.size());
    const int n = m.size();
    U r0 = 0, m0 = 1;
    for (const int i : rep(0, n)) {
        assert(m[i] > 0);
        U r1 = rem_euclid(r[i], m[i]), m1 = m[i];
        if (m0 < m1) {
            std::swap(r0, r1);
            std::swap(m0, m1);
        }
        if (m0 % m1 == 0) {
            if (r0 % m1 != r1) return std::nullopt;
            continue;
        }
        const auto [g, inv_m] = inv_gcd(m0, m1);
        if ((r1 - r0) % g != 0) return std::nullopt;
        const U u1 = m1 / g;
        const U x = (r1 - r0) / g % u1 * inv_m % u1;
        r0 += x * m0;
        m0 *= u1;
        if (r0 < 0) r0 += m0;
    }
    return std::pair<Result, Result>(r0, m0);
}

Back to top page