Computing discrete logarithms

The nmod_discrete_log_pohlig_hellman functions solve g^x \equiv y \pmod p for a prime p, where g is the least primitive root modulo p. The work that depends only on p (factoring p - 1, finding g, building lookup tables) is done once by nmod_discrete_log_pohlig_hellman_precompute_prime, after which each logarithm is computed by nmod_discrete_log_pohlig_hellman_run.

When to use#

The running time depends on the largest prime factors of p - 1, not on the size of p. If p - 1 has only small prime factors (if it is smooth), a logarithm costs a few exponentiations. If p - 1 has a large prime factor q > 10^4, the search for each base-q digit costs up to about q / 50 multiplications, since the lookup tables are limited to 100 entries. The return value of the precomputation is an estimate of the number of multiplications per logarithm, which can be used to decide whether a prime is suitable.

Usage#

1. Initialise an nmod_discrete_log_pohlig_hellman_t with nmod_discrete_log_pohlig_hellman_init.
2. Set it up for p with nmod_discrete_log_pohlig_hellman_precompute_prime. This can be called again with another prime; the old data is freed.
3. Call nmod_discrete_log_pohlig_hellman_run for any number of y \in [1, p). Each result x satisfies 0 \le x < p - 1.
4. Free the memory with nmod_discrete_log_pohlig_hellman_clear.

The base is chosen internally and returned by nmod_discrete_log_pohlig_hellman_primitive_root. To take logarithms to another base h, use that h^x = y if and only if x \log_g h \equiv \log_g y \pmod{p - 1}, a linear congruence that nmod_divides solves (see the second example).

Algorithm#

For each prime power q^e exactly dividing p - 1, the problem is mapped to the subgroup of order q^e by raising y to the power (p-1)/q^e. There the logarithm modulo q^e is found one base-q digit at a time, each digit by a baby-step giant-step search in a table of at most 100 powers of g^{(p-1)/q}. The results are combined by the Chinese remainder theorem.

For many logarithms to a fixed base, also modulo prime powers, the dlog module amortises a larger precomputation; see Discrete logarithm tables.

Examples#

Logarithms to the base of the least primitive root, which is 5 for p = 1000033:

#include <flint/flint.h>
#include <flint/nmod.h>

int main(void)
{
    nmod_discrete_log_pohlig_hellman_t L;
    ulong p = 1000033;                  /* p - 1 = 2^5 * 3 * 11 * 947 */
    ulong ys[3] = { 2, 1000, 1000032 };
    ulong g, x;
    int i;

    nmod_discrete_log_pohlig_hellman_init(L);
    nmod_discrete_log_pohlig_hellman_precompute_prime(L, p);
    g = nmod_discrete_log_pohlig_hellman_primitive_root(L);
    flint_printf("base g = %wu\n", g);

    for (i = 0; i < 3; i++)
    {
        x = nmod_discrete_log_pohlig_hellman_run(L, ys[i]);
        flint_printf("log(%wu) = %wu, g^x = %wu\n", ys[i], x,
                     nmod_pow_ui(g, x, L->mod));
    }

    nmod_discrete_log_pohlig_hellman_clear(L);
    return 0;
}

base g = 5
log(2) = 773396, g^x = 2
log(1000) = 320127, g^x = 1000
log(1000032) = 500016, g^x = 1000032

The logarithm of -1 = 1000032 is (p-1)/2 = 500016, as it must be.

Logarithms to a base h that is not a primitive root. A solution exists only if y lies in the subgroup generated by h:

#include <flint/flint.h>
#include <flint/nmod.h>

int main(void)
{
    nmod_discrete_log_pohlig_hellman_t L;
    nmod_t mod, ord;
    ulong p = 1000033, h = 9, lh, ly, x;
    ulong ys[2] = { 81, 3 };
    int i;

    nmod_init(&mod, p);
    nmod_init(&ord, p - 1);             /* exponents live modulo p - 1 */

    nmod_discrete_log_pohlig_hellman_init(L);
    nmod_discrete_log_pohlig_hellman_precompute_prime(L, p);
    lh = nmod_discrete_log_pohlig_hellman_run(L, h);

    for (i = 0; i < 2; i++)
    {
        ly = nmod_discrete_log_pohlig_hellman_run(L, ys[i]);

        /* h^x = y  <=>  x * log(h) = log(y) mod p - 1 */
        if (nmod_divides(&x, ly, lh, ord))
            flint_printf("%wu^%wu = %wu mod %wu\n", h, x, nmod_pow_ui(h, x, mod), p);
        else
            flint_printf("%wu is not a power of %wu mod %wu\n", ys[i], h, p);
    }

    nmod_discrete_log_pohlig_hellman_clear(L);
    return 0;
}

9^2 = 81 mod 1000033
3 is not a power of 9 mod 1000033