Congruence Equation Partial Solutions

Question

I want to evaluate $x^2 = a\, mod\, n$ for large modulus. But I need only the first few solutions. How can I get it?

I tried following

Reduce[x^2 == 7946761, x, Modulus -> 130356633908760178920]

I got this output :

ReplaceAll::reps: {ToRules[SystemException[RootsListSizeExceeded,100000]]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

ReplaceAll::reps: {x==ToRules[SystemException[RootsListSizeExceeded,100000]]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

I also tried this

PowerModList[7946761, 1/2, 130356633908760178920]

with the output:

SystemException["PowerModListSizeExceeded", 1048576]

I don't need the complete set of solutions. Just first 10-15 are sufficient. How can I access them? Help please.

Answer 1

Using v10.4.1, Reduce returns 2^17=131072 solutions beginning with these 5.

With[{a = 7946761, m = 130356633908760178920},
      x /. {ToRules[Reduce[x^2 == a, x, Modulus -> m][[Range[5]]]]}]

(* {2819, 639260064875759, 2448281311653851, 3275200757643631, \
    3967030855412441} *)

Your modulus m is composite, with prime factorization containing the first 16 primes.

FactorInteger[130356633908760178920]
(* {{2, 3}, {3, 1}, {5, 1}, {7, 1}, {11, 1}, {13, 1}, {17, 1}, {19, 1}, 
    {23, 1}, {29, 1}, {31, 1}, {37, 1}, {41, 1}, {43, 1}, {47, 1}, {53, 1}} *)

The strategy is to use ChineseRemainder[r,k] for all moduli k equal to the prime powers in m=130356633908760178920.

Use PowerModList[a,1/2,k] to find the remainders r for each individual prime-power modulus k.

With[{a = 7946761, m = 130356633908760178920},
     Table[PowerModList[a, 1/2, k], {k, Apply[Power, FactorInteger[m], 1]}]]

The full solution.

CRTsolution[a_, m_] :=
   Block[{f = FactorInteger[m], k, r},
         k = Apply[Power, f, 1];
         r = Flatten[
                Outer[List, Apply[Sequence, Map[PowerModList[a, 1/2, #] &, k]]],
                Length[f] - 1];
         Map[ChineseRemainder[#, k] &, r]]

The answers from CRTsolution agree with Reduce. The hope is that, by breaking the modulus apart, CRTsolution will work with your version of Mathematica.

With[{a = 7946761, m = 130356633908760178920},
     Total[
        Sort[CRTsolution[a, m]] - (x /. {ToRules[Reduce[x^2==a, x, Modulus->m]]})]]
(* 0 *)