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? – bill s Jan 17 '17 at 5:51
• What exactly are you trying to do that has you computing such huge modular roots? – J. M. will be back soon Jan 19 '17 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)$. – caesar Jan 19 '17 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? – caesar Jan 19 '17 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. – J. M. will be back soon Jan 19 '17 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? – caesar Jan 19 '17 at 9:46
• "I need only a single root." - you really ought to have mentioned that in your question to begin with. :| – J. M. will be back soon Jan 19 '17 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? – caesar Jan 19 '17 at 22:54