68 lines
1.4 KiB
Plaintext
68 lines
1.4 KiB
Plaintext
(* Legendre symbol *)
|
|
LegendreSymbol[a_, p_] := PowerMod[a, (p - 1)/2, p]
|
|
|
|
(* Tonelli-Shanks Algorithm *)
|
|
Tonelli[n_, p_] := Module[
|
|
{q, s, z, c, r, t, m, t2, i, b},
|
|
(* Check that n is a quadratic residue mod p *)
|
|
If[LegendreSymbol[n, p] =!= 1,
|
|
Message[Tonelli::notsquare, n, p]; Abort[]];
|
|
|
|
(* Step 1: Factor out powers of 2 from p - 1 *)
|
|
q = p - 1;
|
|
s = 0;
|
|
While[EvenQ[q], q = q/2; s++];
|
|
|
|
(* Special case when s == 1 *)
|
|
If[s == 1,
|
|
Return[PowerMod[n, (p + 1)/4, p]]
|
|
];
|
|
|
|
(* Find a non-residue z *)
|
|
For[z = 2, z < p, z++,
|
|
If[LegendreSymbol[z, p] == p - 1, Break[]]
|
|
];
|
|
|
|
(* Initialize variables *)
|
|
c = PowerMod[z, q, p];
|
|
r = PowerMod[n, (q + 1)/2, p];
|
|
t = PowerMod[n, q, p];
|
|
m = s;
|
|
|
|
(* Main loop *)
|
|
While[Mod[t - 1, p] =!= 0,
|
|
t2 = Mod[t^2, p];
|
|
For[i = 1, i < m, i++,
|
|
If[Mod[t2 - 1, p] === 0, Break[]];
|
|
t2 = Mod[t2^2, p];
|
|
];
|
|
b = PowerMod[c, 2^(m - i - 1), p];
|
|
r = Mod[r*b, p];
|
|
c = Mod[b^2, p];
|
|
t = Mod[t*c, p];
|
|
m = i;
|
|
];
|
|
r
|
|
]
|
|
|
|
(* Test cases *)
|
|
ttest = {
|
|
{10, 13},
|
|
{56, 101},
|
|
{1030, 10009},
|
|
{44402, 100049},
|
|
{665820697, 1000000009},
|
|
{881398088036, 1000000000039},
|
|
{41660815127637347468140745042827704103445750172002, 10^50 + 577}
|
|
};
|
|
|
|
Do[
|
|
n = ttest[[i, 1]];
|
|
p = ttest[[i, 2]];
|
|
r = Tonelli[n, p];
|
|
Assert[Mod[r^2 - n, p] == 0];
|
|
Print["n = ", n, ", p = ", p];
|
|
Print["\troots : ", r, " ", p - r],
|
|
{i, Length[ttest]}
|
|
]
|