Skip to content

Log mod

c++

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
long long mul_modp(long long a, long long b){
    // a*b mod p
    return ((a % MOD) * (b % MOD)) % MOD;
}

long long pow_modp(long long a, long long n){
    // a^n mod p
    long long res = 1;
    while (n > 0) {
        if (n & 1) res = mul_modp(res, a);
        a = mul_modp(a, a);
        n >>= 1;
    }
    return res;
}

long long modinv(long long a){
    // a^{-1} = 1/a mod p (拡張Euclidの互除法)
    long long b = MOD, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= MOD;
    if (u < 0) u += MOD;
    return u;
}

long long log_modp(long long a, long long b) {
    // a^x ≡ b (mod. m) となる最小の正の整数 x を求める
    a %= MOD, b %= MOD;

    // calc sqrt{M}
    long long lo = -1, hi = MOD;
    while (hi - lo > 1) {
        long long mid = (lo + hi) / 2;
        if (mid * mid >= MOD) hi = mid;
        else lo = mid;
    }
    long long sqrtM = hi;

    // {a^0, a^1, a^2, ..., a^sqrt(m)} 
    map<long long, long long> apow;
    long long amari = 1;
    for (long long r = 0; r < sqrtM; ++r) {
        if (!apow.count(amari)) apow[amari] = r;
        (amari *= a) %= MOD;
    }

    // check each A^p
    long long A = pow_modp(modinv(a), sqrtM);
    amari = b;
    for (long long q = 0; q < sqrtM; ++q) {
        if (apow.count(amari)) {
            long long res = q * sqrtM + apow[amari];
            if (res > 0) return res;
        }
        (amari *= A) %= MOD;
    }

    // no solutions
    return -1;
}