Congruence Equation Partial Solutions

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.

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]]]}]

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


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

FactorInteger
(* {{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 *)