rcpl

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub ruthen71/rcpl

:heavy_check_mark: 線形ディオファントス方程式 ($ ax + by = c $) / 線形合同式 ( $ ax \equiv b \pmod m $ )
(math/linear_diophantine.hpp)

線形合同式は $ \bmod m $ で $ a^{-1} $ を求めることでは解が求められない場合がある

例えば、$ 9x \equiv 3 \pmod {15} $ は $ 9^{-1} $ となる値を持たないが、$ x \equiv 2 \pmod 5 $ が解である

このような場合、拡張ユークリッドの互除法を用いることで解が求められる

Depends on

Required by

Verified with

Code

#pragma once

#include <tuple>
#include "math/extended_gcd.hpp"

// Reference: https://ja.wikipedia.org/wiki/ベズーの等式
// ax + by = c を解く (a, b >= 0)
// return {x, y, gcd(a, b), has_solution}
// 解が存在するとき
// (1) a = 0, b = 0, c = 0 のとき
//   x, y は任意
//   {0, 0, gcd(a, b) = 0, 1} を返す
// (2) a = 0, c が b の倍数のとき
//   x は任意, y = c / b
//   {0, c / b, gcd(a, b) = b, 1} を返す
// (3) b = 0, c が a の倍数のとき
//   y は任意, x = c / b
//   {c / a, 0, gcd(a, b) = a, 1} を返す
// (4) a > 0, b > 0, c % gcd(a, b) = 0 のとき
//   x = x' + k * (b / gcd(a, b)), y = y' - k * (a / gcd(a, b))
//   {x', y', gcd(a, b), 1} を返す
// 解が存在しないとき
// {-1, -1, -1, 0} を返す
template <class T> std::tuple<T, T, T, int> linear_diophantine(T a, T b, T c) {
    assert(a >= 0 and b >= 0);
    if (a == 0 and b == 0) {
        if (c == 0) return {0, 0, 0, 1};
        return {-1, -1, -1, 0};
    }
    if (a == 0) {
        // by = c
        if (c % b == 0) return {0, c / b, b, 1};
        return {-1, -1, -1, 0};
    }
    if (b == 0) {
        // ax = c
        if (c % a == 0) return {c / a, 0, a, 1};
        return {-1, -1, -1, 0};
    }
    // as + bt = gcd(a, b) から ax + by = c を求める
    // x = s * (c / gcd(a, b)), y = t * (c / gcd(a, b)) よりも x, y が小さくなる?
    // c = c' + a * dx + b * dy となるように c' を求める
    // (a, b は gcd(a, b) の倍数なので c' は gcd(a, b) の倍数)
    // x = dx + s * (c' / gcd(a, b)), y = dy + t * (c' / gcd(a, b)) が解となる
    // ax + by = a * dx + b * dy + (as + bt) * (c' / gcd(a, b)) = a * dx + b * dy + c' = c
    auto [s, t, g] = extended_gcd(a, b);
    if (c % g != 0) return {-1, -1, -1, 0};
    T dx = c / a;
    c -= dx * a;
    T dy = c / b;
    c -= dy * b;
    T x = dx + s * (c / g);
    T y = dy + t * (c / g);
    return {x, y, g, 1};
}

// 線形合同式 ax = b (mod m) を解く (m > 0)
// 解が存在する場合 x = x' (mod h) となるときの最小の x' と h を返す
// 解が存在しない場合 {-1, -1} を返す
template <class T> std::pair<T, T> linear_congruence(T a, T b, T m) {
    assert(m > 0);
    a = (a % m + m) % m;
    b = (b % m + m) % m;
    // ax = b (mod m) <=> ax + my = b となる (x, y) が存在
    auto [x, y, g, is_ok] = linear_diophantine(a, m, b);
    if (!is_ok) return {-1, -1};
    T h = m / g;
    x = (x % h + h) % h;
    return {x, h};
}
#line 2 "math/linear_diophantine.hpp"

#include <tuple>
#line 2 "math/extended_gcd.hpp"

#line 4 "math/extended_gcd.hpp"
// find (x, y) s.t. ax + by = gcd(a, b)
// a, b >= 0
// return {x, y, gcd(a, b)}
template <class T> std::tuple<T, T, T> extended_gcd(T a, T b) {
    if (b == 0) return {1, 0, a};
    auto [y, x, g] = extended_gcd(b, a % b);
    return {x, y - (a / b) * x, g};
}
#line 5 "math/linear_diophantine.hpp"

// Reference: https://ja.wikipedia.org/wiki/ベズーの等式
// ax + by = c を解く (a, b >= 0)
// return {x, y, gcd(a, b), has_solution}
// 解が存在するとき
// (1) a = 0, b = 0, c = 0 のとき
//   x, y は任意
//   {0, 0, gcd(a, b) = 0, 1} を返す
// (2) a = 0, c が b の倍数のとき
//   x は任意, y = c / b
//   {0, c / b, gcd(a, b) = b, 1} を返す
// (3) b = 0, c が a の倍数のとき
//   y は任意, x = c / b
//   {c / a, 0, gcd(a, b) = a, 1} を返す
// (4) a > 0, b > 0, c % gcd(a, b) = 0 のとき
//   x = x' + k * (b / gcd(a, b)), y = y' - k * (a / gcd(a, b))
//   {x', y', gcd(a, b), 1} を返す
// 解が存在しないとき
// {-1, -1, -1, 0} を返す
template <class T> std::tuple<T, T, T, int> linear_diophantine(T a, T b, T c) {
    assert(a >= 0 and b >= 0);
    if (a == 0 and b == 0) {
        if (c == 0) return {0, 0, 0, 1};
        return {-1, -1, -1, 0};
    }
    if (a == 0) {
        // by = c
        if (c % b == 0) return {0, c / b, b, 1};
        return {-1, -1, -1, 0};
    }
    if (b == 0) {
        // ax = c
        if (c % a == 0) return {c / a, 0, a, 1};
        return {-1, -1, -1, 0};
    }
    // as + bt = gcd(a, b) から ax + by = c を求める
    // x = s * (c / gcd(a, b)), y = t * (c / gcd(a, b)) よりも x, y が小さくなる?
    // c = c' + a * dx + b * dy となるように c' を求める
    // (a, b は gcd(a, b) の倍数なので c' は gcd(a, b) の倍数)
    // x = dx + s * (c' / gcd(a, b)), y = dy + t * (c' / gcd(a, b)) が解となる
    // ax + by = a * dx + b * dy + (as + bt) * (c' / gcd(a, b)) = a * dx + b * dy + c' = c
    auto [s, t, g] = extended_gcd(a, b);
    if (c % g != 0) return {-1, -1, -1, 0};
    T dx = c / a;
    c -= dx * a;
    T dy = c / b;
    c -= dy * b;
    T x = dx + s * (c / g);
    T y = dy + t * (c / g);
    return {x, y, g, 1};
}

// 線形合同式 ax = b (mod m) を解く (m > 0)
// 解が存在する場合 x = x' (mod h) となるときの最小の x' と h を返す
// 解が存在しない場合 {-1, -1} を返す
template <class T> std::pair<T, T> linear_congruence(T a, T b, T m) {
    assert(m > 0);
    a = (a % m + m) % m;
    b = (b % m + m) % m;
    // ax = b (mod m) <=> ax + my = b となる (x, y) が存在
    auto [x, y, g, is_ok] = linear_diophantine(a, m, b);
    if (!is_ok) return {-1, -1};
    T h = m / g;
    x = (x % h + h) % h;
    return {x, h};
}
Back to top page