How to solve equations in Gaussian integers modulo p

How to solve a complex equation of the form:

$$z^n \equiv i \pmod p$$

where $z$ is a Gaussian integer, $i$ is the imaginary number, $n, p \in \mathbb{Z}^+$ and $p$ is prime.

I am dealing with quite big numbers, for example: $z^{8192} \equiv i~\pmod {2^{61}-1}$. This is just an example, i.e. $p$ need not be a Mersenne prime.

I tried the binomial expansion of $z^n=(a+bi)^n$ and solving the real and imaginary parts individually, but it takes a very long time (it's been 4 hours and Mathematica hasn't solved for the real part yet).

I tried this code:

Clear[z, n, p, a, b, r, t, eq1, eq2, system];
z = a + b*I;
n = 4096;
p = 18446744069414584321;
exp = Expand[z^n];
r = Expand[Re[exp]];
t = Expand[Im[exp]];
eq1 = Refine[r, a ∈ Integers && b ∈ Integers];
eq2 = Refine[t, a ∈ Integers && b ∈ Integers];
system = eq1 == 0 && eq2 == 1;
sol1 = Solve[system, {a, b}, Modulus -> p]


This code is very slow for big numbers. Any suggestions to improve it?

• Have you tried it for small numbers? What code have you tried? Commented Jan 17, 2017 at 5:51
• What exactly are you trying to do that has you computing such huge modular roots? Commented Jan 19, 2017 at 9:29
• I need to verify that the discrete Galois transform can work with the prime number indicated above. I tried the transform with small primes $(257, 8191)$ and small $(n=2)$. It worked fine. Now I need to use it in larger settings shown above for $(n=4k, 8k, 16k, 32k, 64k)$. Commented Jan 19, 2017 at 9:40

Here is a rather naïve implementation; I'm sure there are more clever/efficient approaches:

modRoot[z_, {p_Integer?PrimeQ, r_Integer?Positive}, m_Integer?PrimeQ] :=
Nest[Select[Union[Flatten[(\[FormalX] + I \[FormalY] /.
Solve[{Sum[(-1)^j Binomial[p, 2 j] \[FormalX]^(p - 2 j)
\[FormalY]^(2 j), {j, 0, Quotient[p, 2]}] == Re[#],
Sum[(-1)^j Binomial[p, 2 j + 1] \[FormalX]^(p - 2 j - 1)
\[FormalY]^(2 j + 1), {j, 0, Quotient[p - 1, 2]}] ==
Im[#]},
{\[FormalX], \[FormalY]}, Modulus -> m]) & /@ #]], NumberQ] &, z, r]

gaussianPowerModList[z_, r_Rational, m_Integer?PrimeQ] :=
Fold[modRoot[#, #2, m] &,
{PowerMod[z, Numerator[r], m]}, FactorInteger[Denominator[r]]]


It takes a while on your particular example:

AbsoluteTiming[sols = gaussianPowerModList[I, 1/2^13, 2^61 - 1];]
{521.019, Null}

Length[sols]
8192

Take[sols, 3]
{146709914666800 + 1089528222731080675 I,
403020492323839 + 87412432372606530 I,
1417152643549164 + 256379023454276512 I}

PowerMod[%, 2^13, 2^61 - 1]
{I, I, I}


A slightly less difficult example:

gaussianPowerModList[1 - I, 1/(2 3 5), 2^19 - 1]
{26299 + 165240 I, 68075 + 434841 I, 94374 + 75794 I, 429913 + 448493 I,
456212 + 89446 I, 497988 + 359047 I}

PowerMod[%, 2 3 5, 2^19 - 1]
{1 - I, 1 - I, 1 - I, 1 - I, 1 - I, 1 - I}

• Calling the function sols1 = gaussianPowerModList[I, 8192, 18446744069414584321];] runs out of memory after more than 10 hours of execution. Is there any way it can be computed more efficiently? Commented Jan 19, 2017 at 9:41
• You'll notice that you get exactly 8192 roots since you have a prime modulus. Do you really need all those 8192 roots? If so, you could modify the method here so that they are computed a few at a time instead of in a single blow. Commented Jan 19, 2017 at 9:44
• I need only a single root. But I have no idea how to instruct solve to exit after finding a single solution? Commented Jan 19, 2017 at 9:46
• "I need only a single root." - you really ought to have mentioned that in your question to begin with. :| Commented Jan 19, 2017 at 9:48
• I am sorry I forgot to mention that @J.M.. One strange thing is: why it computes the roots modulo $(2^{61}-1)$ fast while working with that prime is very slow and requires large memory? Commented Jan 19, 2017 at 22:54

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]

AlgebraPolynomialRemainderModList::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`.