1
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • $\begingroup$ Have you tried it for small numbers? What code have you tried? $\endgroup$
    – bill s
    Jan 17, 2017 at 5:51
  • $\begingroup$ What exactly are you trying to do that has you computing such huge modular roots? $\endgroup$
    – J. M.'s eventual burnout ♦
    Jan 19, 2017 at 9:29
  • $\begingroup$ 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)$. $\endgroup$
    – caesar
    Jan 19, 2017 at 9:40

2 Answers 2

Reset to default
3
$\begingroup$

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}
Improve this answer
$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$
    – caesar
    Jan 19, 2017 at 9:41
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s eventual burnout ♦
    Jan 19, 2017 at 9:44
  • $\begingroup$ I need only a single root. But I have no idea how to instruct solve to exit after finding a single solution? $\endgroup$
    – caesar
    Jan 19, 2017 at 9:46
  • 1
    $\begingroup$ "I need only a single root." - you really ought to have mentioned that in your question to begin with. :| $\endgroup$
    – J. M.'s eventual burnout ♦
    Jan 19, 2017 at 9:48
  • $\begingroup$ 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? $\endgroup$
    – caesar
    Jan 19, 2017 at 22:54
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.