This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub ruthen71/rcpl
#include "math/linear_diophantine.hpp"
線形合同式は $ \bmod m $ で $ a^{-1} $ を求めることでは解が求められない場合がある
例えば、$ 9x \equiv 3 \pmod {15} $ は $ 9^{-1} $ となる値を持たないが、$ x \equiv 2 \pmod 5 $ が解である
このような場合、拡張ユークリッドの互除法を用いることで解が求められる
#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}; }