subcrip
/
cp-templates
1
0
Fork 0
Code

Update exgcd.cc

This commit is contained in:
Ariel 2024-03-11 20:12:32 +08:00 committed by GitHub
parent 96f5877abb
commit 568886b27d
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
1 changed files with 26 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
number/exgcd.cc
View File

@ -1,19 +1,23 @@
namespace Exgcd { namespace Exgcd {
template <typename T> T abs(T x) { return x < 0 ? -x : x; }
template <typename T>
struct exgcd_solution_t { struct exgcd_solution_t {
ll x, y, gcd; T x, y, gcd;
}; };
template <typename T>
struct diophantine_solution_t { struct diophantine_solution_t {
exgcd_solution_t x_min, y_min; exgcd_solution_t<T> x_min, y_min;
ll range; T range;
}; };
// solve `ax + by = gcd(a, b)` // solve `ax + by = gcd(a, b)`
optional<exgcd_solution_t> exgcd(ll a, ll b) { template <typename T>
optional<exgcd_solution_t<T>> exgcd(T a, T b) {
if (a < 0 || b < 0 || a == 0 && b == 0) return nullopt; if (a < 0 || b < 0 || a == 0 && b == 0) return nullopt;
ll x, y, g; T x, y, g;
function<void(ll, ll)> __exgcd = [&__exgcd, &x, &y, &g] (ll a, ll b) -> void { function<void(T, T)> __exgcd = [&__exgcd, &x, &y, &g] (T a, T b) -> void {
if (b == 0) { if (b == 0) {
g = a, x = 1, y = 0; g = a, x = 1, y = 0;
} else { } else {
@ -26,7 +30,8 @@ namespace Exgcd {
return {{ x, y, g }}; return {{ x, y, g }};
}; };
optional<ll> inverse(ll a, ll b) { template <typename T>
optional<T> inverse(T a, T b) {
auto raw = exgcd(a, b); auto raw = exgcd(a, b);
if (raw == nullopt || raw.value().gcd != 1) { if (raw == nullopt || raw.value().gcd != 1) {
return nullopt; return nullopt;
@ -36,14 +41,15 @@ namespace Exgcd {
} }
// solve { x = a_i (mod n_i) } if n_i's are coprime // solve { x = a_i (mod n_i) } if n_i's are coprime
optional<ll> crt(const vector<pll>& equations) { template <typename T>
ll prod = 1; optional<T> crt(const vector<pair<T, T>>& equations) {
T prod = 1;
for (auto&& [a, n] : equations) { for (auto&& [a, n] : equations) {
prod *= n; prod *= n;
} }
ll res = 0; T res = 0;
for (auto&& [a, n] : equations) { for (auto&& [a, n] : equations) {
ll m = prod / n; T m = prod / n;
auto m_rev = inverse(m, n); auto m_rev = inverse(m, n);
if (m_rev == nullopt) return nullopt; if (m_rev == nullopt) return nullopt;
res = mod(res + a * mod(m * m_rev.value(), prod), prod); res = mod(res + a * mod(m * m_rev.value(), prod), prod);
@ -51,22 +57,24 @@ namespace Exgcd {
return res; return res;
} }
// find minimal non-negative integral solutions of `ax + by = c` // find minimal non-negative integral solutions of `ax + by = c`. It's not guaranteed that the other variable is non-negative.
optional<diophantine_solution_t> diophantine(ll a, ll b, ll c, bool force_positive = false) { template <typename T>
optional<diophantine_solution_t<T>> diophantine(T a, T b, T c, bool force_positive = false) {
if (a < 0 || b < 0 || a == 0 && b == 0) return nullopt; if (a < 0 || b < 0 || a == 0 && b == 0) return nullopt;
auto raw = exgcd(a, b).value(); auto raw = exgcd(a, b).value();
if (c % raw.gcd) { if (c % raw.gcd) {
return nullopt; return nullopt;
} else { } else {
ll x = raw.x * c / raw.gcd, y = raw.y * c / raw.gcd; T x = raw.x * c / raw.gcd, y = raw.y * c / raw.gcd;
ll kx = force_positive ? (x <= 0 ? (-x) * raw.gcd / b + 1 : 1 - (x + b / raw.gcd - 1) * raw.gcd / b) : (x <= 0 ? ((-x) + b / raw.gcd - 1) * raw.gcd / b : (- x * raw.gcd / b)); T kx = force_positive ? (x <= 0 ? (-x) * raw.gcd / b + 1 : 1 - (x + b / raw.gcd - 1) * raw.gcd / b) : (x <= 0 ? ((-x) + b / raw.gcd - 1) * raw.gcd / b : (- x * raw.gcd / b));
ll ky = force_positive ? (y <= 0 ? (- 1 - (-y) * raw.gcd / a) : (y + a / raw.gcd - 1) * raw.gcd / a - 1) : (y <= 0 ? (- ((-y) + a / raw.gcd - 1) * raw.gcd / a) : y * raw.gcd / a); T ky = force_positive ? (y <= 0 ? (- 1 - (-y) * raw.gcd / a) : (y + a / raw.gcd - 1) * raw.gcd / a - 1) : (y <= 0 ? (- ((-y) + a / raw.gcd - 1) * raw.gcd / a) : y * raw.gcd / a);
return {{ { x + b * kx / raw.gcd , y - a * kx / raw.gcd , raw.gcd }, { x + b * ky / raw.gcd , y - a * ky / raw.gcd, raw.gcd }, abs(kx - ky) + 1 }}; return {{ { x + b * kx / raw.gcd , y - a * kx / raw.gcd , raw.gcd }, { x + b * ky / raw.gcd , y - a * ky / raw.gcd, raw.gcd }, abs(kx - ky) + 1 }};
} }
} }
// find the minimal non-negative integral solution of `ax = b (mod n)` // find the minimal non-negative integral solution of `ax = b (mod n)`
optional<ll> congruential(ll a, ll b, ll n) { template <typename T>
optional<T> congruential(T a, T b, T n) {
if (a == 0) { if (a == 0) {
if (b != 0) return nullopt; if (b != 0) return nullopt;
return 0; return 0;