Unverified Commit 7ba1f7ae authored by bjorn-martinsson's avatar bjorn-martinsson Committed by GitHub

More general ModLog (#175)

parent 3f5fd9a3
/**
* Author: Andrei Navumenka, chilli
* Date: 2019-11-08
* License: Unlicense
* Source: https://github.com/indy256/codelibrary/blob/master/cpp/numbertheory/discrete_log.cpp
* Description: Returns the smallest $x \ge 0$ s.t. $a^x = b \pmod m$. a and m must be coprime.
* Author: Bjorn Martinsson
* Date: 2020-06-03
* License: CC0
* Source: own work
* Description: Returns the smallest $x > 0$ s.t. $a^x = b \pmod m$, or
* $-1$ if no such $x$ exists. Note that Modlog(a,1,m) can be used to
* calculate the order of $a$.
* Time: $O(\sqrt m)$
* Status: tested for all 0 <= a,x,m < 200.
* Status: tested for all 0 <= a,x < 500 and 0 < m < 500.
*
* Details: This algorithm uses the baby-step giant-step method to
* find (i,j) such that a^(n * i) = b * a^j (mod m), where n > sqrt(m)
* and 0 < i, j <= n. If a and m are coprime then a^j has a modular
* inverse, which means that a^(i * n - j) = b (mod m$).
*
* However this particular implementation of baby-step giant-step works even
* without assuming a and m are coprime, using the following idea:
*
* Assume p^x is a prime divisor of m. Then we have 3 cases
* 1. b is divisible by p^x
* 2. b is divisible only by some p^y, 0<y<x
* 3. b is not divisible by p
* The important thing to note is that in case 2, modLog(a,b,m) (if
* it exists) cannot be > sqrt(m), (technically it cannot be >= log2(m)).
* So once all exponenents of a that are <= sqrt(m) has been checked, you
* cannot have case 2. Case 2 is the only tricky case.
*
* So the modification allowing for non-coprime input invloves checking all
* exponents of a that are <= n, and then handling the non-tricky cases by
* a simple gcd(a^n,m) == gcd(b,m) check.
*/
#pragma once
ll modLog(ll a, ll b, ll m) {
assert(__gcd(a, m) == 1);
ll n = (ll) sqrt(m) + 1, e = 1, x = 1, res = LLONG_MAX;
unordered_map<ll, ll> f;
rep(i,0,n) e = e * a % m;
rep(i,0,n) x = x * e % m, f.emplace(x, i + 1);
rep(i,0,n) if (f.count(b = b * a % m))
res = min(res, f[b] * n - i - 1);
return res;
ll n = (ll) sqrt(m) + 1, e = 1, f = 1, j = 1;
unordered_map<ll, ll> A;
while (j <= n && (e = f = e * a % m) != b % m)
A[e * b % m] = j++;
if (e == b % m) return j;
if (__gcd(m, e) == __gcd(m, b))
rep(i,2,n+2) if (A.count(e = e * f % m))
return n * i - A[e];
return -1;
}
......@@ -3,29 +3,24 @@
#include "../../content/number-theory/ModLog.h"
int main() {
const int lim = 100;
rep(m,1,lim) {
rep(a,0,lim) {
if (__gcd(a, m) != 1) continue;
vector<ll> ans(m, -1);
ll b = 1 % m;
rep(x,0,m) {
if (ans[b] == -1) ans[b] = x;
b = b * a % m;
}
rep(b,0,lim) {
ll res = modLog(a, b, m);
ll b2 = b % m;
if (ans[b2] == -1) assert(res == LLONG_MAX);
else {
if (ans[b2] != res) {
cerr << "FAIL" << endl;
cerr << "Expected log(" << a << ", " << b << ", " << m << ") = " << ans[b2] << ", found " << res << endl;
return 1;
}
}
}
}
}
cout<<"Tests passed!"<<endl;
const int lim = 100;
rep(m,1,lim) {
rep(a,0,lim) {
vector<ll> ans(m, -1);
ll b = a % m;
rep(x,1,max(m,2)) {
if (ans[b] == -1) ans[b] = x;
b = b * a % m;
}
rep(b,0,m) {
ll res = modLog(a, b, m);
if (ans[b] != res) {
cerr << "FAIL" << endl;
cerr << "Expected log(" << a << ", " << b << ", " << m << ") = " << ans[b] << ", found " << res << endl;
return 1;
}
}
}
}
cout<<"Tests passed!"<<endl;
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment