Cipolla's algorithm: Difference between revisions

=={{header|Mathematica}}/{{header|Wolfram Language}}==

<lang Mathematica>ClearAll[Cipolla]

Cipolla[n_, p_] := Module[{ls, omega2, nn, a, r, s},

ls = JacobiSymbol[n, p];

If[ls != 1,

{0, 0, False}

,

a = 0;

While[True,

omega2 = Mod[a^2 + p - n, p];

If[JacobiSymbol[omega2, p] == -1, Break[]];

a++

];

ClearAll[Mul];

Mul[{ax_, ay_}, {bx_, by_}] :=

Mod[{ax bx + ay by omega2, ax by + bx ay}, p];

r = {1, 0};

s = {a, 1};

nn = Mod[BitShiftRight[p + 1, 1], p];

While[nn > 0,

If[BitAnd[nn, 1] == 1,

r = Mul[r, s]

];

s = Mul[s, s];

nn = BitShiftRight[nn, 1];

];

If[r[[2]] != 0, Return[{0, 0, False}]];

If[Mod[r[[1]]^2, p] != n, Return[{0, 0, False}]];

Return[{r[[1]], p - r[[1]], True}]

]

]

Cipolla[10, 13]

Cipolla[56, 101]

Cipolla[8218, 10007]

Cipolla[8219, 10007]

Cipolla[331575, 1000003]

Cipolla[665165880, 1000000007]

Cipolla[881398088036, 1000000000039]

Cipolla[34035243914635549601583369544560650254325084643201, 10^50 + 151]</lang>

{{out}}

<pre>{6,7,True}

{37,64,True}

{9872,135,True}

{0,0,False}

{855842,144161,True}

{475131702,524868305,True}

{791399408049,208600591990,True}

{82563118828090362261378993957450213573687113690751,17436881171909637738621006042549786426312886309400,True}</pre>