Discrete logarithms modulo a multiprecision prime

These functions solve $\alpha^x \equiv y \pmod p$ for a prime $p$ of any size and a fixed primitive root $\alpha$, by the method of Pohlig and Hellman [PohHel1978]: the logarithm is computed modulo each prime power $q^e$ dividing $p - 1$, one base-$q$ digit at a time with Shanks' baby-step giant-step search in the subgroup of order $q$, and the results are combined by the Chinese remainder theorem. This is fast when every prime factor $q$ of $p - 1$ is small; the cost grows like $\sqrt{q}$ for the largest $q$. The method and its analysis are described in Discrete logarithms modulo a prime (Pohlig–Hellman); the word-size version is nmod_discrete_log_pohlig_hellman_run.

Usage#

1. Initialise an fmpz_mod_discrete_log_pohlig_hellman_t with fmpz_mod_discrete_log_pohlig_hellman_init.
2. Call fmpz_mod_discrete_log_pohlig_hellman_precompute_prime with the prime $p$. It factors $p - 1$, finds the least primitive root $\alpha$, which fmpz_mod_discrete_log_pohlig_hellman_primitive_root returns, and precomputes the tables. The return value estimates the number of multiplications modulo $p$ per logarithm, which a caller can use to decide whether $p$ is suitable.
3. Call fmpz_mod_discrete_log_pohlig_hellman_run for each nonzero residue $y$; the result is the logarithm in $[0, p - 1)$. The precomputed object is not modified, so several threads may share it.
4. Free the memory with fmpz_mod_discrete_log_pohlig_hellman_clear.

The function does not check that $p$ is prime or that $y$ is nonzero and reduced. The precomputation throws an exception if $p - 1$ has a prime factor that does not fit in a word, and the search throws one if it finds no logarithm, which can only happen for invalid input. The factorisation of $p - 1$ dominates the precomputation unless $p - 1$ is known to be smooth; fmpz_next_smooth_prime provides primes for which it is, with $p - 1$ a product of powers of 2, 5, 7, 11, 13, 17, 19 and 23.

The tables for a prime $q$ have at most 100 entries: a giant step is about $\sqrt{q}$ for $q$ below about $10^4$ and is doubled until the table is small enough for larger $q$, so the number of baby steps grows from $\sqrt{q}$ to about $q/50$ for large $q$.

Example#

The program takes the next smooth prime after $10^{30}$ (a 128-bit prime; the table of smooth primes has no entries between $2^{64}$ and $2^{127}$), computes a power of the primitive root and recovers the exponent.

#include <flint/fmpz.h>
#include <flint/fmpz_mod.h>

int main(void)
{
    fmpz_t p, x, y, e;
    fmpz_mod_ctx_t ctx;
    fmpz_mod_discrete_log_pohlig_hellman_t L;
    double cost;

    fmpz_init(p);
    fmpz_init(x);
    fmpz_init(y);
    fmpz_init(e);

    /* a prime p > 10^30 with p - 1 smooth */
    fmpz_set_str(x, "1000000000000000000000000000000", 10);
    if (!fmpz_next_smooth_prime(p, x))
        return 1;
    printf("p = "); fmpz_print(p); printf("\n");

    fmpz_mod_discrete_log_pohlig_hellman_init(L);
    cost = fmpz_mod_discrete_log_pohlig_hellman_precompute_prime(L, p);
    printf("base = ");
    fmpz_print(fmpz_mod_discrete_log_pohlig_hellman_primitive_root(L));
    printf(", cost estimate %s 2000\n", cost < 2000 ? "<" : ">=");

    /* y = base^e mod p, then recover e */
    fmpz_mod_ctx_init(ctx, p);
    fmpz_set_str(e, "123456789012345678901234567890", 10);
    fmpz_mod_pow_fmpz(y, fmpz_mod_discrete_log_pohlig_hellman_primitive_root(L), e, ctx);
    fmpz_mod_discrete_log_pohlig_hellman_run(x, L, y);
    printf("log = "); fmpz_print(x); printf("\n");

    fmpz_mod_ctx_clear(ctx);
    fmpz_mod_discrete_log_pohlig_hellman_clear(L);
    fmpz_clear(p);
    fmpz_clear(x);
    fmpz_clear(y);
    fmpz_clear(e);
    return 0;
}

p = 201315753943344863644036821430284967937
base = 3, cost estimate < 2000
log = 123456789012345678901234567890