Since $p$ is prime, this equation only needs to be solved within finite field. Depending on $p$, if $p \equiv 1 \space mod \space 4$, then $\sqrt{-1}$ has good reduction on $\mathbb{F_p}$, gaussian integers are not needed. Otherwise $\sqrt{-1}$ exists in $\mathbb{F_{p^2}}$. A general method would be finding discrete logarithm in respective finite field, then solving linear congruence $log(z^n) \space\equiv\space n log(z)$in either $\mathbb{Z/(p-1)Z}$ or $\mathbb{Z/(p^2-1)Z}$. Index calculus algorithm can be used to crack discrete logarithm over finite field, with subexponential complexity wrt. $log(p)$.
Mathematica don't come with strong utilities dealing finite fields. Mod
and PowerMod
both work with gaussian integers. PowerMod
with rational arguments can be used to find roots modulo integers. However it won't work well with fields extensions. (Bug?)
In[1]:= PowerMod[2,1/8192,2^61-1]
Out[1]= 2305843009213562879
In[2]:= PowerMod[I, 1/8192,2^61-1]
Algebra`PolynomialRemainderModList::unipoly: -- Message text not found -- ({-I, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, <<7973>>, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1})
Interrupt>
Personally I'd just call PARI via some python bridge.
GaussianIntegerQ[n_] := AllTrue[ReIm[n], IntegerQ]
DiscreteLog[n_Integer, pr_Integer,
p_Integer /; PrimeQ[p] && Mod[p, 4] == 1] :=
ExternalEvaluate["Python",
TemplateApply[
"import cypari2;pari=cypari2.Pari();Mod=pari.Mod;int(pari.znlog(\
Mod(`n`,`p`),Mod(`pr`,`p`)))", <|"n" -> n, "pr" -> pr, "p" -> p|>]]
DiscreteLog[n_Integer, p_Integer /; PrimeQ[p] && Mod[p, 4] == 1] :=
DiscreteLog[n, PrimitiveRoot[p], p]
DiscreteLog[n_ /; GaussianIntegerQ[n], pr_ /; GaussianIntegerQ[pr],
p_Integer /; PrimeQ[p] && Mod[p, 4] == 3] :=
ExternalEvaluate["Python", TemplateApply["
import cypari2
pari=cypari2.Pari()
p=`p`
ii=pari.ffgen(pari.Mod(pari('x^2=1'),p))
n=`na`+`nb`*ii
pr=`pra`+`prb`*ii
int(pari.fflog(n,pr))",
With[{nn = ReIm@n, prr = ReIm@pr}, <|"na" -> nn[[1]],
"nb" -> nn[[2]], "pra" -> prr[[1]], "prb" -> prr[[2]],
"p" -> p|>]]]
GFp2PrimitiveRoot[p_Integer /; PrimeQ[p] && Mod[p, 4] == 3] :=
ToExpression[ExternalEvaluate["Python", TemplateApply["
import cypari2
pari=cypari2.Pari()
ii=pari.ffgen(pari.Mod(pari('x^2+1'),`p`))
str(pari.ffprimroot(ii))", <|"p" -> p|>]]] /. x -> I
DiscreteLog[n_ /; GaussianIntegerQ[n],
p_Integer /; PrimeQ[p] && Mod[p, 4] == 3] :=
ExternalEvaluate["Python", TemplateApply["
import cypari2
pari=cypari2.Pari()
ii=pari.ffgen(pari.Mod(pari('x^2+1'),`p`))
n=`na`+`nb`*ii
pr=pari.ffprimroot(ii)
int(pari.fflog(n,pr))",
With[{nn = ReIm@n}, <|"na" -> nn[[1]], "nb" -> nn[[2]],
"p" -> p|>]]]
Some tests:
p = 2^61 - 1;
pr = GFp2PrimitiveRoot[2^61 - 1];
logx = DiscreteLog[I, p]; (*8192 divides logx*)
y = PowerMod[pr, logx/8192, p]
This gives -782078743845564464 + 729159760183691423 I
in less than 1 second. To verity, PowerMod[y, 8192, p]
gives I
.