问题引入
求解 x a ≡ b ( m o d m ) x^a \equiv b(mod~~m) xa≡b(mod m)。
离散对数
离散对数是基于同余运算和原根的对数运算,当模 m m m有原根时,设 g g g为其中一个原根,则当 x ≡ g k ( m o d m ) x \equiv g^k(mod~~m) x≡gk(mod m)时,可以对两边求对数,这样问题变成了 l o g g x ≡ k ( m o d φ ( m ) ) log_gx \equiv k(mod~~\varphi(m)) loggx≡k(mod φ(m))。
当 m m m为质数时
设 m m m为质数 p p p,由于 p p p是质数,那么 p p p一定存在一个原根 g g g,根据质数原根的性质,因此对于模 p p p意义下的任何数 x ( 0 ≤ x < p ) x(0\leq x <p) x(0≤x<p)有且仅有一个数 i ( 0 ≤ i < p − 1 ) i(0 \leq i <p-1) i(0≤i<p−1)满足 x ≡ g i x \equiv g^i x≡gi。
方法一
令 x = g c x=g^c x=gc( g g g是 p p p的一个原根且 g , c g,c g,c一定存在),问题转化为求解 ( g c ) a ≡ b ( m o d p ) (g^c)^a \equiv b(mod~~p) (gc)a≡b(mod p),然后可以得到:
( g a ) c ≡ b (g^a)^c \equiv b (ga)c≡b
这样问题就变成了一个朴素 B S G S BSGS BSGS问题,可以在 O ( p ) O(\sqrt{p}) O(p)求出 c c c,这样得到了原方程的一个特解 x 0 ≡ g c ( m o d p ) x_0\equiv g^c(mod~~p) x0≡gc(mod p)
方法二
令 x = g c x=g^c x=gc,并且设 b = g t b=g^t b=gt,那么得到 g a c ≡ g t ( m o d p ) g^{ac} \equiv g^t(mod~~p) gac≡gt(mod p),两边同时取离散对数得到:
a c ≡ t ( m o d φ ( p ) ) ac \equiv t (mod~~\varphi(p)) ac≡t(mod φ(p))
这时可以通过 B S G S BSGS BSGS求解 g t ≡ b ( m o d p ) g^t \equiv b(mod~~p) gt≡b(mod p)得到 t t t,那么就转化成了线性同余方程组的问题,这样最后得到的也是一个特解 x 0 ≡ g c ( m o d p ) x_0 \equiv g^c(mod~~p) x0≡gc(mod p)
求出所有解
已知 g φ ( p ) ≡ 1 ( m o d p ) g^{\varphi(p)}\equiv 1(mod~~p) gφ(p)≡1(mod p),那么可以得到:
$x^a \equiv g^{ac+\varphi§t}\equiv b(mod~~p), \forall t \in \Z $
于是所有的解为:
x ≡ g c + φ ( p ) t a ( m o d p ) , ∀ t ∈ Z , a ∣ φ ( p ) t x \equiv g^{c+\frac{\varphi(p)t}{a}}(mod~~p), \forall t \in \Z,a|\varphi(p)t x≡gc+aφ(p)t(mod p),∀t∈Z,a∣φ(p)t
显然有 a g c d ( a , φ ( p ) ) ∣ t \frac{a}{gcd(a,\varphi(p))}|t gcd(a,φ(p))a∣t,设 t = a g c d ( a , φ ( p ) ) i t=\frac{a}{gcd(a,\varphi(p))}i t=gcd(a,φ(p))ai,那么得到原问题的所有解为:
x ≡ g c + φ ( p ) i g c d ( a , φ ( p ) ) ( m o d p ) , ∀ i ∈ Z x \equiv g^{c+\frac{\varphi(p)i}{gcd(a,\varphi(p))}}(mod~~p), \forall i \in \Z x≡gc+gcd(a,φ(p))φ(p)i(mod p),∀i∈Z
代码
以方法一为例,极限数据范围,求出的是 [ 0 , p ) [0,p) [0,p)范围内的所有解,即一共 p p p个上述 i i i的取值范围是 [ 0 , g c d ( a , φ ( p ) ) [0,gcd(a,\varphi(p)) [0,gcd(a,φ(p))。
bitset<maxn> vis;
unordered_map<ll, ll> mp;
ll prime[maxn], fac[maxn], ans[maxn];
int cnt, tot;
void euler() {
vis.reset(), cnt = 0;
for (int i = 2; i < maxn; i++) {
if (!vis[i]) prime[++cnt] = i;
for (int j = 1; j <= cnt && 1LL * i * prime[j] < maxn; j++) {
vis[i * prime[j]] = 1;
if (i % prime[j] == 0) break;
}
}
}
void getFac(ll n) {
tot = 0;
for (int i = 1; prime[i] * prime[i] <= n; i++) {
if (n % prime[i] == 0) {
fac[++tot] = prime[i];
while (n % prime[i] == 0) n /= prime[i];
}
}
if (n > 1) fac[++tot] = n;
}
ll slow_mul(ll a, ll b, ll p) {
ll ans = 0;
a %= p;
while (b) {
if (b & 1) ans = (ans + a) % p;
a = (a + a) % p;
b >>= 1;
}
return ans;
}
ll qkp(ll x, ll n, ll p) {
ll ans = 1;
x %= p;
while (n) {
if (n & 1) ans = slow_mul(ans, x, p);
x = slow_mul(x, x, p);
n >>= 1;
}
return ans;
}
inline bool check(ll g, ll p) {
for (int i = 1; i <= tot; i++) {
ll t = (p - 1) / fac[i];
if (qkp(g, t, p) == 1) return 0;
}
return 1;
}
ll getRoot(ll p) {
getFac(p - 1);
for (ll g = 2; g < p; g++) {
if (check(g, p)) return g;
}
}
inline ll gcd(ll a, ll b) {
return b == 0 ? a : gcd(b, a % b);
}
inline ll exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) {
x = 1, y = 0;
return a;
}
ll g = exgcd(b, a % b, y, x);
y -= a / b * x;
return g;
}
ll getInv(ll a, ll p) {
return qkp(a, p - 2, p);
}
ll bsgs(ll a, ll b, ll p) {
mp.clear();
a %= p, b %= p;
if (b == 1 || p == 1) return 0;
int m = sqrt(p) + 1;
ll cur = b;
for (int i = 0; i < m; i++) {
if (!mp.count(cur)) mp[cur] = i;
cur = slow_mul(cur, a, p);
}
ll x = qkp(a, m, p), ans = -1, s = 1;
for (int i = 1; i <= m; i++) {
s = slow_mul(s, x, p);
if (mp.count(s)) {
ans = 1LL * i * m - mp[s];
break;
}
}
return ans;
}
bool solve(ll a, ll b, ll p) {
ll g = getRoot(p);
ll t = bsgs(g, b, p), phi = p - 1;
ll d = gcd(a, phi);
if (t < 0 || t % d) return 0;
ll x, y;
exgcd(a, phi, x, y);
t /= d, phi /= d;
x = (x * t % phi + phi) % phi;
ans[0] = qkp(g, x, p);
for (int i = 1; i < d; i++) {
x += phi;
ans[i] = qkp(g, x, p);
}
sort(ans, ans + d);
for (int i = 0; i < d; i++) printf("%lld\n", ans[i]);
return 1;
}
当 m m m为合数时
求出 [ 0 , m ) [0,m) [0,m)范围内的解,模板中 m o d mod mod即为 m m m
ll A, B, mod;
ll pow(ll x, ll y, ll mod = 0, ll ans = 1) {
if (mod) {
for (; y; y >>= 1, x = x * x % mod)
if (y & 1) ans = ans * x % mod;
} else {
for (; y; y >>= 1, x = x * x)
if (y & 1) ans = ans * x;
}
return ans;
}
struct factor {
ll prime[20], expo[20], pk[20], tot;
void factor_lleger(ll n) {
tot = 0;
for (ll i = 2; i * i <= n; ++i)
if (n % i == 0) {
prime[tot] = i, expo[tot] = 0, pk[tot] = 1;
do ++expo[tot], pk[tot] *= i; while ((n /= i) % i == 0);
++tot;
}
if (n > 1) prime[tot] = n, expo[tot] = 1, pk[tot++] = n;
}
ll phi(ll id) const {
return pk[id] / prime[id] * (prime[id] - 1);
}
} mods, _p;
ll p_inverse(ll x, ll id) {
assert(x % mods.prime[id] != 0);
return pow(x, mods.phi(id) - 1, mods.pk[id]);
}
void exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) x = 1, y = 0;
else exgcd(b, a % b, y, x), y -= a / b * x;
}
ll inverse(ll x, ll mod) {
assert(__gcd(x, mod) == 1);
ll ret, tmp;
exgcd(x, mod, ret, tmp), ret %= mod;
return ret + (ret >> 31 & mod);
}
vector<ll> sol[20];
void solve_2(ll id, ll k) {
ll mod = 1 << k;
if (k == 0) {
sol[id].emplace_back(0);
return;
} else {
solve_2(id, k - 1);
vector<ll> t;
for (ll s : sol[id]) {
if (!((pow(s, A) ^ B) & mod - 1))
t.emplace_back(s);
if (!((pow(s | 1 << k - 1, A) ^ B) & mod - 1))
t.emplace_back(s | 1 << k - 1);
}
swap(sol[id], t);
}
}
ll BSGS(ll B, ll g, ll mod) {
unordered_map<ll, ll> map;
ll L = ceil(sqrt(mod)), t = 1;
for (ll i = 1; i <= L; ++i) {
t = t * g % mod;
map[B * t % mod] = i;
}
ll now = 1;
for (ll i = 1; i <= L; ++i) {
now = now * t % mod;
if (map.count(now)) return i * L - map[now];
}
assert(0);
}
ll find_primitive_root(ll id) {
ll phi = mods.phi(id);
_p.factor_lleger(phi);
auto check = [&](ll g) {
for (ll i = 0; i < _p.tot; ++i)
if (pow(g, phi / _p.prime[i], mods.pk[id]) == 1)
return 0;
return 1;
};
for (ll g = 2; g < mods.pk[id]; ++g) if (check(g)) return g;
assert(0);
}
void division(ll id, ll a, ll b, ll mod) { // ax = b (mod M)
ll M = mod, g = __gcd(__gcd(a, b), mod);
a /= g, b /= g, mod /= g;
if (__gcd(a, mod) > 1) return;
ll t = b * inverse(a, mod) % mod;
for (; t < M; t += mod) sol[id].emplace_back(t);
}
void solve_p(ll id, ll B = ::B) {
ll p = mods.prime[id], e = mods.expo[id], mod = mods.pk[id];
if (B % mod == 0) {
ll q = pow(p, (e + A - 1) / A);
for (ll t = 0; t < mods.pk[id]; t += q)
sol[id].emplace_back(t);
} else if (B % p != 0) {
ll phi = mods.phi(id);
ll g = find_primitive_root(id), z = BSGS(B, g, mod);
division(id, A, z, phi);
for (ll &x : sol[id]) x = pow(g, x, mod);
} else {
ll q = 0;
while (B % p == 0) B /= p, ++q;
ll pq = pow(p, q);
if (q % A != 0) return;
mods.expo[id] -= q, mods.pk[id] /= pq;
solve_p(id, B);
mods.expo[id] += q, mods.pk[id] *= pq;
if (!sol[id].size()) return;
ll s = pow(p, q - q / A);
ll t = pow(p, q / A);
ll u = pow(p, e - q);
vector<ll> res;
for (ll y : sol[id]) {
for (ll i = 0; i < s; ++i)
res.emplace_back((i * u + y) * t);
}
swap(sol[id], res);
}
}
vector<ll> allans;
void dfs(ll dep, ll ans, ll mod) {
if (dep == mods.tot) {
allans.emplace_back(ans);
return;
}
ll p = mods.pk[dep], k = p_inverse(mod % p, dep);
for (ll a : sol[dep]) {
ll nxt = (a - ans % p + p) * k % p * mod + ans;
dfs(dep + 1, nxt, mod * p);
}
}
int solve() {
mods.factor_lleger(mod);
allans.clear();
for (ll i = mods.tot - 1; ~i; --i) {
sol[i].clear();
mods.prime[i] == 2 ? solve_2(i, mods.expo[i]) : solve_p(i);
if (!sol[i].size()) return 0;
}
dfs(0, 0, 1);
sort(allans.begin(), allans.end());
cout << allans.size() << '\n';
for (auto i : allans) cout << i << ' ';
cout << '\n';
return 1;
}