# Fast calculation of discrete logarithms

## Does Mathematica have any built-in fast algorithms for calculating discrete logarithms over $(\mathbb{Z}_p)^\times$ (the group of integers modulo $p$)?

Essentially, for a fixed large prime p, a generator g, and an integer y, I want to quickly compute the x such that PowerMod[g, x, p] == y.

Currently the only built-in method I've found is MultiplicativeOrder[g, p, y]; however, for large values of p (let's say p = 10000223 for now), MultiplicativeOrder can take over a second to calculate a single value (although for some values it returns almost instantly).

Module[{p = 10000223, g, n = 20},
g = PrimitiveRoot[p];
Grid @ MapAt[NumberForm[#, {∞, 3}] &,
Table[AbsoluteTiming[
PowerMod[g, MultiplicativeOrder[g, p, i], p]], {i, n}], {;; , 1}]
]

(*
0.000   1
1.198   2
0.801   3
1.098   4
0.000   5
0.691   6
0.668   7
0.968   8
0.287   9
1.217   10
1.217   11
0.569   12
0.776   13
0.557   14
0.796   15
0.855   16
1.059   17
0.166   18
0.110   19
1.062   20
*)


This is unacceptably slow. Although my p is large by some measures, cryptographically speaking it is minuscule at only 24 bits long.

## Custom function for with precomputation for fixed $p$

In fact, generating a table of all discrete logs with respect to one generator takes just a few seconds:

Module[{p = 10000223, g},
g = PrimitiveRoot[p];
AbsoluteTiming[
Ordering @ NestList[Mod[# g, p] &, g, p - 2];
];
]

(* {2.26656, Null} *)


We can actually use that one table to calculate discrete logs to any base:

\begin{align*} x &= \log_{g'} y\\ y &= (g')^x \bmod p \\ y &= \left(g^{\log_g g'}\right)^x \bmod p \\ y &= g^{x\log_g g'\bmod p-1} \bmod p \\ x\log_g g' &= \log_g y \bmod p-1 \\ x &= (\log_g y)(\log_g g')^{-1} \bmod p-1 \end{align*}

Both the logarithms are just a table lookup, and computing the modular inverse is extremely fast. With that in mind, here's the code I'm actually using to calculate discrete logarithms:

p = 10000223;
g = PrimitiveRoot[p];
lookup = Ordering @ NestList[Mod[# g, p] &, g, p - 2];

discreteLog[b_, y_] := Mod[lookup[[y]] PowerMod[lookup[[b]], -1, p - 1], p - 1]


Testing against MultiplicativeOrder:

Table[
With[{y = RandomInteger[{1, p - 1}], g = randomGenerator[]},
MultiplicativeOrder[g, p, y] == discreteLog[g, y]
],
{10}
]

(* {True, True..., True} *)


This implementation of discreteLog is extremely fast (about 300k per second), using just two table lookups, one modular inverse (computed quickly with e.g. the extended Euclidian algorithm) and one modular multiplication.

However, if I move to a slightly larger $p$ (say $p\approx 10^9$ or larger) then the lookup table becomes unwieldy. This is where I'd like a relatively fast builtin function.

• +1, good question. There's no faster "built-in" I'm aware of, and problem is O(ugly) I'm sure you're aware... – ciao Sep 27 '15 at 2:54
• @ciao Basically I was hoping for something like an implementation of baby-step/giant-step or Pollard rho. I could do them myself, but probably not quite as fast as a native implementation. I doubt my lookup table will actually outgrow my memory, but I thought it worth asking. – 2012rcampion Sep 27 '15 at 3:09
• Forgot to add - might take a gander at the FiniteFields package. Perhaps the cost of setting a field list to use with FieldInd might be worthwhile (not holding my breath, though...) – ciao Sep 27 '15 at 5:12
• @ciao Neat, I forgot about that package. Unfortunately it looks like before I can use FieldInd I need to set PowerListQ[GF[p]] = True, which seems to perform a similar precomputation step, and it actually does it much slower (probably because it's designed to handle fields with more complex structure). – 2012rcampion Sep 27 '15 at 5:22
• Yeah, I recall the setup being a bit of a time/memory pig... – ciao Sep 27 '15 at 5:29