Cipolla's algorithm: Difference between revisions
m
→{{header|Wren}}: Minor tidy
Thundergnat (talk | contribs) (→{{header|Perl 6}}: Add Perl 6) |
m (→{{header|Wren}}: Minor tidy) |
||
| (48 intermediate revisions by 21 users not shown) | |||
Line 75:
* [[Tonelli-Shanks algorithm]]
<br><br>
=={{header|11l}}==
{{trans|Python}}
<syntaxhighlight lang="11l">F convertToBase(n, b)
I (n < 2)
R [n]
V temp = n
V ans = [Int]()
L (temp != 0)
ans = [temp % b] [+] ans
temp I/= b
R ans
F cipolla(=n, p)
n %= p
I (n == 0 | n == 1)
R [n, (p - n) % p]
V phi = p - 1
I (pow(BigInt(n), phi I/ 2, p) != 1)
R [Int]()
I (p % 4 == 3)
V ans = Int(pow(BigInt(n), (p + 1) I/ 4, p))
R [ans, (p - ans) % p]
V aa = 0
L(i) 1 .< p
V temp = pow(BigInt((i * i - n) % p), phi I/ 2, p)
I (temp == phi)
aa = i
L.break
V exponent = convertToBase((p + 1) I/ 2, 2)
F cipollaMult(ab, cd, w, p)
V (a, b) = ab
V (c, d) = cd
R ((a * c + b * d * w) % p, (a * d + b * c) % p)
V x1 = (aa, 1)
V x2 = cipollaMult(x1, x1, aa * aa - n, p)
L(i) 1 .< exponent.len
I (exponent[i] == 0)
x2 = cipollaMult(x2, x1, aa * aa - n, p)
x1 = cipollaMult(x1, x1, aa * aa - n, p)
E
x1 = cipollaMult(x1, x2, aa * aa - n, p)
x2 = cipollaMult(x2, x2, aa * aa - n, p)
R [x1[0], (p - x1[0]) % p]
print(‘Roots of 2 mod 7: ’cipolla(2, 7))
print(‘Roots of 8218 mod 10007: ’cipolla(8218, 10007))
print(‘Roots of 56 mod 101: ’cipolla(56, 101))
print(‘Roots of 1 mod 11: ’cipolla(1, 11))
print(‘Roots of 8219 mod 10007: ’cipolla(8219, 10007))</syntaxhighlight>
{{out}}
<pre>
Roots of 2 mod 7: [4, 3]
Roots of 8218 mod 10007: [9872, 135]
Roots of 56 mod 101: [37, 64]
Roots of 1 mod 11: [1, 10]
Roots of 8219 mod 10007: []
</pre>
=={{header|C}}==
{{trans|FreeBasic}}
<syntaxhighlight lang="c">#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
struct fp2 {
int64_t x, y;
};
uint64_t randULong(uint64_t min, uint64_t max) {
uint64_t t = (uint64_t)rand();
return min + t % (max - min);
}
// returns a * b mod modulus
uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t modulus) {
uint64_t x = 0, y = a % modulus;
while (b > 0) {
if ((b & 1) == 1) {
x = (x + y) % modulus;
}
y = (y << 1) % modulus;
b = b >> 1;
}
return x;
}
//returns b ^^ power mod modulus
uint64_t pow_mod(uint64_t b, uint64_t power, uint64_t modulus) {
uint64_t x = 1;
while (power > 0) {
if ((power & 1) == 1) {
x = mul_mod(x, b, modulus);
}
b = mul_mod(b, b, modulus);
power = power >> 1;
}
return x;
}
// miller-rabin prime test
bool isPrime(uint64_t n, int64_t k) {
uint64_t a, x, n_one = n - 1, d = n_one;
uint32_t s = 0;
uint32_t r;
if (n < 2) {
return false;
}
// limit 2^63, pow_mod/mul_mod can't handle bigger numbers
if (n > 9223372036854775808ull) {
printf("The number is too big, program will end.\n");
exit(1);
}
if ((n % 2) == 0) {
return n == 2;
}
while ((d & 1) == 0) {
d = d >> 1;
s = s + 1;
}
while (k > 0) {
k = k - 1;
a = randULong(2, n);
x = pow_mod(a, d, n);
if (x == 1 || x == n_one) {
continue;
}
for (r = 1; r < s; r++) {
x = pow_mod(x, 2, n);
if (x == 1) return false;
if (x == n_one) goto continue_while;
}
if (x != n_one) {
return false;
}
continue_while: {}
}
return true;
}
int64_t legendre_symbol(int64_t a, int64_t p) {
int64_t x = pow_mod(a, (p - 1) / 2, p);
if ((p - 1) == x) {
return x - p;
} else {
return x;
}
}
struct fp2 fp2mul(struct fp2 a, struct fp2 b, int64_t p, int64_t w2) {
struct fp2 answer;
uint64_t tmp1, tmp2;
tmp1 = mul_mod(a.x, b.x, p);
tmp2 = mul_mod(a.y, b.y, p);
tmp2 = mul_mod(tmp2, w2, p);
answer.x = (tmp1 + tmp2) % p;
tmp1 = mul_mod(a.x, b.y, p);
tmp2 = mul_mod(a.y, b.x, p);
answer.y = (tmp1 + tmp2) % p;
return answer;
}
struct fp2 fp2square(struct fp2 a, int64_t p, int64_t w2) {
return fp2mul(a, a, p, w2);
}
struct fp2 fp2pow(struct fp2 a, int64_t n, int64_t p, int64_t w2) {
struct fp2 ret;
if (n == 0) {
ret.x = 1;
ret.y = 0;
return ret;
}
if (n == 1) {
return a;
}
if ((n & 1) == 0) {
return fp2square(fp2pow(a, n / 2, p, w2), p, w2);
} else {
return fp2mul(a, fp2pow(a, n - 1, p, w2), p, w2);
}
}
void test(int64_t n, int64_t p) {
int64_t a, w2;
int64_t x1, x2;
struct fp2 answer;
printf("Find solution for n = %lld and p = %lld\n", n, p);
if (p == 2 || !isPrime(p, 15)) {
printf("No solution, p is not an odd prime.\n\n");
return;
}
//p is checked and is a odd prime
if (legendre_symbol(n, p) != 1) {
printf(" %lld is not a square in F%lld\n\n", n, p);
return;
}
while (true) {
do {
a = randULong(2, p);
w2 = a * a - n;
} while (legendre_symbol(w2, p) != -1);
answer.x = a;
answer.y = 1;
answer = fp2pow(answer, (p + 1) / 2, p, w2);
if (answer.y != 0) {
continue;
}
x1 = answer.x;
x2 = p - x1;
if (mul_mod(x1, x1, p) == n && mul_mod(x2, x2, p) == n) {
printf("Solution found: x1 = %lld, x2 = %lld\n\n", x1, x2);
return;
}
}
}
int main() {
srand((size_t)time(0));
test(10, 13);
test(56, 101);
test(8218, 10007);
test(8219, 10007);
test(331575, 1000003);
test(665165880, 1000000007);
//test(881398088036, 1000000000039);
return 0;
}</syntaxhighlight>
{{out}}
<pre>Find solution for n = 10 and p = 13
Solution found: x1 = 7, x2 = 6
Find solution for n = 56 and p = 101
Solution found: x1 = 37, x2 = 64
Find solution for n = 8218 and p = 10007
Solution found: x1 = 9872, x2 = 135
Find solution for n = 8219 and p = 10007
8219 is not a square in F10007
Find solution for n = 331575 and p = 1000003
Solution found: x1 = 144161, x2 = 855842
Find solution for n = 665165880 and p = 1000000007
Solution found: x1 = 524868305, x2 = 475131702</pre>
=={{header|C sharp|C#}}==
<syntaxhighlight lang="csharp">using System;
using System.Numerics;
namespace CipollaAlgorithm {
class Program {
static readonly BigInteger BIG = BigInteger.Pow(10, 50) + 151;
private static Tuple<BigInteger, BigInteger, bool> C(string ns, string ps) {
BigInteger n = BigInteger.Parse(ns);
BigInteger p = ps.Length > 0 ? BigInteger.Parse(ps) : BIG;
// Legendre symbol. Returns 1, 0, or p-1
BigInteger ls(BigInteger a0) => BigInteger.ModPow(a0, (p - 1) / 2, p);
// Step 0: validate arguments
if (ls(n) != 1) {
return new Tuple<BigInteger, BigInteger, bool>(0, 0, false);
}
// Step 1: Find a, omega2
BigInteger a = 0;
BigInteger omega2;
while (true) {
omega2 = (a * a + p - n) % p;
if (ls(omega2) == p - 1) {
break;
}
a += 1;
}
// Multiplication in Fp2
BigInteger finalOmega = omega2;
Tuple<BigInteger, BigInteger> mul(Tuple<BigInteger, BigInteger> aa, Tuple<BigInteger, BigInteger> bb) {
return new Tuple<BigInteger, BigInteger>(
(aa.Item1 * bb.Item1 + aa.Item2 * bb.Item2 * finalOmega) % p,
(aa.Item1 * bb.Item2 + bb.Item1 * aa.Item2) % p
);
}
// Step 2: Compute power
Tuple<BigInteger, BigInteger> r = new Tuple<BigInteger, BigInteger>(1, 0);
Tuple<BigInteger, BigInteger> s = new Tuple<BigInteger, BigInteger>(a, 1);
BigInteger nn = ((p + 1) >> 1) % p;
while (nn > 0) {
if ((nn & 1) == 1) {
r = mul(r, s);
}
s = mul(s, s);
nn >>= 1;
}
// Step 3: Check x in Fp
if (r.Item2 != 0) {
return new Tuple<BigInteger, BigInteger, bool>(0, 0, false);
}
// Step 5: Check x * x = n
if (r.Item1 * r.Item1 % p != n) {
return new Tuple<BigInteger, BigInteger, bool>(0, 0, false);
}
// Step 4: Solutions
return new Tuple<BigInteger, BigInteger, bool>(r.Item1, p - r.Item1, true);
}
static void Main(string[] args) {
Console.WriteLine(C("10", "13"));
Console.WriteLine(C("56", "101"));
Console.WriteLine(C("8218", "10007"));
Console.WriteLine(C("8219", "10007"));
Console.WriteLine(C("331575", "1000003"));
Console.WriteLine(C("665165880", "1000000007"));
Console.WriteLine(C("881398088036", "1000000000039"));
Console.WriteLine(C("34035243914635549601583369544560650254325084643201", ""));
}
}
}</syntaxhighlight>
{{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>
=={{header|D}}==
{{trans|Kotlin}}
<syntaxhighlight lang="d">import std.bigint;
import std.stdio;
import std.typecons;
enum BIGZERO = BigInt(0);
/// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
BigInt modPow(BigInt b, BigInt e, BigInt n) {
if (n == 1) return BIGZERO;
BigInt result = 1;
b = b % n;
while (e > 0) {
if (e % 2 == 1) {
result = (result * b) % n;
}
e >>= 1;
b = (b*b) % n;
}
return result;
}
alias Point = Tuple!(BigInt, "x", BigInt, "y");
alias Triple = Tuple!(BigInt, "x", BigInt, "y", bool, "b");
Triple c(string ns, string ps) {
auto n = BigInt(ns);
BigInt p;
if (ps.length > 0) {
p = BigInt(ps);
} else {
p = BigInt(10)^^50 + 151;
}
// Legendre symbol, returns 1, 0 or p - 1
auto ls = (BigInt a) => modPow(a, (p-1)/2, p);
// Step 0, validate arguments
if (ls(n) != 1) return Triple(BIGZERO, BIGZERO, false);
// Step 1, find a, omega2
auto a = BIGZERO;
BigInt omega2;
while (true) {
omega2 = (a * a + p - n) % p;
if (ls(omega2) == p-1) break;
a++;
}
// multiplication in Fp2
auto mul = (Point aa, Point bb) => Point(
(aa.x * bb.x + aa.y * bb.y * omega2) % p,
(aa.x * bb.y + bb.x * aa.y) % p
);
// Step 2, compute power
auto r = Point(BigInt(1), BIGZERO);
auto s = Point(a, BigInt(1));
auto nn = ((p+1) >> 1) % p;
while (nn > 0) {
if ((nn & 1) == 1) r = mul(r, s);
s = mul(s, s);
nn >>= 1;
}
// Step 3, check x in Fp
if (r.y != 0) return Triple(BIGZERO, BIGZERO, false);
// Step 5, check x * x = n
if (r.x*r.x%p!=n) return Triple(BIGZERO, BIGZERO, false);
// Step 4, solutions
return Triple(r.x, p-r.x, true);
}
void main() {
writeln(c("10", "13"));
writeln(c("56", "101"));
writeln(c("8218", "10007"));
writeln(c("8219", "10007"));
writeln(c("331575", "1000003"));
writeln(c("665165880", "1000000007"));
writeln(c("881398088036", "1000000000039"));
writeln(c("34035243914635549601583369544560650254325084643201", ""));
}</syntaxhighlight>
{{out}}
<pre>Tuple!(BigInt, "x", BigInt, "y", bool, "b")(6, 7, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(37, 64, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(9872, 135, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(0, 0, false)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(855842, 144161, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(475131702, 524868305, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(791399408049, 208600591990, true)
Tuple!(BigInt, "x", BigInt, "y", bool, "b")(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400, true)</pre>
=={{header|EchoLisp}}==
<
(lib 'struct)
(lib 'types)
Line 129 ⟶ 579:
(assert (mod= n (* x x) p)) ;; checking the result
(printf "Roots of %d are (%d,%d) (mod %d)" n x (% (- p x) p) p))
</syntaxhighlight>
{{out}}
<pre>
Line 150 ⟶ 600:
</pre>
=={{header|F_Sharp|F#}}==
===The function===
This task uses [http://www.rosettacode.org/wiki/Extensible_prime_generator#The_function Extensible Prime Generator (F#)]
<syntaxhighlight lang="fsharp">
// Cipolla's algorithm. Nigel Galloway: June 16th., 2019
let Cipolla n g =
let rec fN i g e l=match e with n when n=0I->i |_ when e%2I=1I->fN ((i*g)%l) ((g*g)%l) (e/2I) l |_-> fN i ((g*g)%l) (e/2I) l
let rec fG g=match (n/g+g)>>>1 with n when bigint.Abs(g-n)>>>1<2I->n+1I |g->fG g
let a,b=let rec fI i=let q=i*i-n in if fN 1I q ((g-1I)/2I) g>1I then (i,q) else fI (i+1I) in fI(fG (bigint(sqrt(double n))))
let fE=Seq.unfold(fun(n,i)->Some((n,i),((n*n+i*i*b)%g,(2I*n*i)%g)))(a,1I)|>Seq.cache
let rec fL Πn Πi α β=match 2I**α with
Ω when Ω<β->fL Πn Πi (α+1) β
|Ω when Ω>β->let n,i=Seq.item (α-1) fE in fL ((Πn*n+Πi*i*b)%g) ((Πn*i+Πi*n)%g) 0 (β-Ω/2I)
|_->let n,i=Seq.item α fE in ((Πn*n+Πi*i*b)%g)
if fN 1I n ((g-1I)/2I) g<>1I then None else Some(fL 1I 0I 0 ((g+1I)/2I))
</syntaxhighlight>
===The Task===
<syntaxhighlight lang="fsharp">
let test=[(10I,13I);(56I,101I);(8218I,10007I);(8219I,10007I);(331575I,1000003I);(665165880I,1000000007I);(881398088036I,1000000000039I);(34035243914635549601583369544560650254325084643201I,10I**50+151I)]
test|>List.iter(fun(n,g)->match Cipolla n g with Some r->printfn "Cipolla %A %A -> %A (%A) check %A" n g r (g-r) ((r*r)%g) |_->printfn "Cipolla %A %A -> has no result" n g)
</syntaxhighlight>
{{out}}
<pre>
Cipolla 10 13 -> 7 (6) check 10
Cipolla 56 101 -> 64 (37) check 56
Cipolla 8218 10007 -> 135 (9872) check 8218
Cipolla 8219 10007 -> has no result
Cipolla 331575 1000003 -> 144161 (855842) check 331575
Cipolla 665165880 1000000007 -> 475131702 (524868305) check 665165880
Cipolla 881398088036 1000000000039 -> 208600591990 (791399408049) check 881398088036
Cipolla 34035243914635549601583369544560650254325084643201 100000000000000000000000000000000000000000000000151 -> 17436881171909637738621006042549786426312886309400 (82563118828090362261378993957450213573687113690751) check 34035243914635549601583369544560650254325084643201
Real: 00:00:00.089, CPU: 00:00:00.090, GC gen0: 2, gen1: 0
</pre>
=={{header|Factor}}==
{{trans|Sidef}}
{{works with|Factor|0.99 2020-08-14}}
<syntaxhighlight lang="factor">USING: accessors assocs interpolate io kernel literals locals
math math.extras math.functions ;
TUPLE: point x y ;
C: <point> point
:: (cipolla) ( n p -- m )
0 0 :> ( a! ω2! )
[ ω2 p legendere -1 = ]
[ a sq n - p rem ω2! a 1 + a! ] do until a 1 - a!
[| a b |
a x>> b x>> * a y>> b y>> ω2 * * + p mod
a x>> b y>> * b x>> a y>> * + p mod <point>
] :> [mul]
1 0 <point> :> r!
a 1 <point> :> s!
p 1 + -1 shift :> n!
[ n 0 > ] [
n odd? [ r s [mul] call r! ] when
s s [mul] call s!
n -1 shift n!
] while
r y>> zero? r x>> f ? ;
: cipolla ( n p -- m/f )
2dup legendere 1 = [ (cipolla) ] [ 2drop f ] if ;
! Task
{
{ 10 13 }
{ 56 101 }
{ 8218 10007 }
{ 8219 10007 }
{ 331575 1000003 }
{ 665165880 1000000007 }
{ 881398088036 1000000000039 }
${
34035243914635549601583369544560650254325084643201
10 50 ^ 151 +
}
}
[
2dup cipolla
[ 2dup - [I Roots of ${3} are (${1} ${0}) mod ${2}I] ]
[ [I No solution for (${}, ${})I] ] if* nl
] assoc-each</syntaxhighlight>
{{out}}
<pre>
Roots of 10 are (6 7) mod 13
Roots of 56 are (37 64) mod 101
Roots of 8218 are (9872 135) mod 10007
No solution for (8219, 10007)
Roots of 331575 are (855842 144161) mod 1000003
Roots of 665165880 are (475131702 524868305) mod 1000000007
Roots of 881398088036 are (791399408049 208600591990) mod 1000000000039
Roots of 34035243914635549601583369544560650254325084643201 are (82563118828090362261378993957450213573687113690751 17436881171909637738621006042549786426312886309400) mod 100000000000000000000000000000000000000000000000151
</pre>
=={{header|FreeBASIC}}==
===LongInt version===
Had a close look at the EchoLisp code for step 2.
Used the FreeBASIC code from the Miller-Rabin task for prime testing.
<syntaxhighlight lang="freebasic">' version 08-04-2017
' compile with: fbc -s console
' maximum for p is 17 digits to be on the save side
' TRUE/FALSE are built-in constants since FreeBASIC 1.04
' But we have to define them for older versions.
#Ifndef TRUE
#Define FALSE 0
#Define TRUE Not FALSE
#EndIf
Type fp2
x As LongInt
y As LongInt
End Type
Function mul_mod(a As ULongInt, b As ULongInt, modulus As ULongInt) As ULongInt
' returns a * b mod modulus
Dim As ULongInt x, y = a mod modulus
While b > 0
If (b And 1) = 1 Then
x = (x + y) Mod modulus
End If
y = (y Shl 1) Mod modulus
b = b Shr 1
Wend
Return x
End Function
Function pow_mod(b As ULongInt, power As ULongInt, modulus As ULongInt) As ULongInt
' returns b ^ power mod modulus
Dim As ULongInt x = 1
While power > 0
If (power And 1) = 1 Then
' x = (x * b) Mod modulus
x = mul_mod(x, b, modulus)
End If
' b = (b * b) Mod modulus
b = mul_mod(b, b, modulus)
power = power Shr 1
Wend
Return x
End Function
Function Isprime(n As ULongInt, k As Long) As Long
' miller-rabin prime test
If n > 9223372036854775808ull Then ' limit 2^63, pow_mod/mul_mod can't handle bigger numbers
Print "number is to big, program will end"
Sleep
End
End If
' 2 is a prime, if n is smaller then 2 or n is even then n = composite
If n = 2 Then Return TRUE
If (n < 2) OrElse ((n And 1) = 0) Then Return FALSE
Dim As ULongInt a, x, n_one = n - 1, d = n_one
Dim As UInteger s
While (d And 1) = 0
d = d Shr 1
s = s + 1
Wend
While k > 0
k = k - 1
a = Int(Rnd * (n -2)) +2 ' 2 <= a < n
x = pow_mod(a, d, n)
If (x = 1) Or (x = n_one) Then Continue While
For r As Integer = 1 To s -1
x = pow_mod(x, 2, n)
If x = 1 Then Return FALSE
If x = n_one Then Continue While
Next
If x <> n_one Then Return FALSE
Wend
Return TRUE
End Function
Function legendre_symbol (a As LongInt, p As LongInt) As LongInt
Dim As LongInt x = pow_mod(a, ((p -1) \ 2), p)
If p -1 = x Then
Return x - p
Else
Return x
End If
End Function
Function fp2mul(a As fp2, b As fp2, p As LongInt, w2 As LongInt) As fp2
Dim As fp2 answer
Dim As ULongInt tmp1, tmp2
' needs to be broken down in smaller steps to avoid overflow
' answer.x = (a.x * b.x + a.y * b.y * w2) Mod p
' answer.y = (a.x * b.y + a.y * b.x) Mod p
tmp1 = mul_mod(a.x, b.x, p)
tmp2 = mul_mod(a.y, b.y, p)
tmp2 = mul_mod(tmp2, w2, p)
answer.x = (tmp1 + tmp2) Mod p
tmp1 = mul_mod(a.x, b.y, p)
tmp2 = mul_mod(a.y, b.x, p)
answer.y = (tmp1 + tmp2) Mod p
Return answer
End Function
Function fp2square(a As fp2, p As LongInt, w2 As LongInt) As fp2
Return fp2mul(a, a, p, w2)
End Function
Function fp2pow(a As fp2, n As LongInt, p As LongInt, w2 As LongInt) As fp2
If n = 0 Then Return Type (1, 0)
If n = 1 Then Return a
If n = 2 Then Return fp2square(a, p, w2)
If (n And 1) = 0 Then
Return fp2square(fp2pow(a, n \ 2, p, w2), p , w2)
Else
Return fp2mul(a, fp2pow(a, n -1, p, w2), p, w2)
End If
End Function
' ------=< MAIN >=------
Data 10, 13, 56, 101, 8218, 10007,8219, 10007
Data 331575, 1000003, 665165880, 1000000007
Data 881398088036, 1000000000039
Randomize Timer
Dim As LongInt n, p, a, w2
Dim As LongInt i, x1, x2
Dim As fp2 answer
For i = 1 To 7
Read n, p
Print
Print "Find solution for n =";n ; " and p =";p
If p = 2 OrElse Isprime(p,15) = FALSE Then
Print "No solution, p is not a odd prime"
Continue For
End If
' p is checked and is a odd prime
If legendre_symbol(n, p) <> 1 Then
Print n; " is not a square in F";Str(p)
Continue For
End If
Do
Do
a = Rnd * (p -2) +2
w2 = a * a - n
Loop Until legendre_symbol(w2, p) = -1
answer = Type(a, 1)
answer = fp2pow(answer, (p +1) \ 2, p, w2)
If answer.y <> 0 Then Continue Do
x1 = answer.x : x2 = p - x1
If mul_mod(x1, x1, p) = n AndAlso mul_mod(x2, x2, p) = n Then
Print "Solution found: x1 ="; x1; ", "; "x2 ="; x2
Exit Do
End If
Loop ' loop until solution is found
Next
' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End</syntaxhighlight>
{{out}}
<pre>Find solution for n = 10 and p = 13
Solution found: x1 = 7, x2 = 6
Find solution for n = 56 and p = 101
Solution found: x1 = 37, x2 = 64
Find solution for n = 8218 and p = 10007
Solution found: x1 = 9872, x2 = 135
Find solution for n = 8219 and p = 10007
8219 is not a square in F10007
Find solution for n = 331575 and p = 1000003
Solution found: x1 = 144161, x2 = 855842
Find solution for n = 665165880 and p = 1000000007
Solution found: x1 = 475131702, x2 = 524868305
Find solution for n = 881398088036 and p = 1000000000039
Solution found: x1 = 791399408049, x2 = 208600591990</pre>
===GMP version===
{{libheader|GMP}}
<syntaxhighlight lang="freebasic">' version 12-04-2017
' compile with: fbc -s console
#Include Once "gmp.bi"
Type fp2
x As Mpz_ptr
y As Mpz_ptr
End Type
Data "10", "13"
Data "56", "101"
Data "8218", "10007"
Data "8219", "10007"
Data "331575", "1000003"
Data "665165880", "1000000007"
Data "881398088036", "1000000000039"
Data "34035243914635549601583369544560650254325084643201" ', 10^50 + 151
Function fp2mul(a As fp2, b As fp2, p As Mpz_ptr, w2 As Mpz_ptr) As fp2
Dim As fp2 r
r.x = Allocate(Len(__Mpz_struct)) : Mpz_init(r.x)
r.y = Allocate(Len(__Mpz_struct)) : Mpz_init(r.y)
Mpz_mul (r.x, a.y, b.y)
Mpz_mul (r.x, r.x, w2)
Mpz_addmul(r.x, a.x, b.x)
Mpz_mod (r.x, r.x, p)
Mpz_mul (r.y, a.x, b.y)
Mpz_addmul(r.y, a.y, b.x)
Mpz_mod (r.y, r.y, p)
Return r
End Function
Function fp2square(a As fp2, p As Mpz_ptr, w2 As Mpz_ptr) As fp2
Return fp2mul(a, a, p, w2)
End Function
Function fp2pow(a As fp2, n As Mpz_ptr, p As Mpz_ptr, w2 As Mpz_ptr) As fp2
If Mpz_cmp_ui(n, 0) = 0 Then
Mpz_set_ui(a.x, 1)
Mpz_set_ui(a.y, 0)
Return a
End If
If Mpz_cmp_ui(n, 1) = 0 Then Return a
If Mpz_cmp_ui(n, 2) = 0 Then Return fp2square(a, p, w2)
If Mpz_tstbit(n, 0) = 0 Then
Mpz_fdiv_q_2exp(n, n, 1) ' even
Return fp2square(fp2pow(a, n, p, w2), p, w2)
Else
Mpz_sub_ui(n, n, 1) ' odd
Return fp2mul(a, fp2pow(a, n, p, w2), p, w2)
End If
End Function
' ------=< MAIN >=------
Dim As Long i
Dim As ZString Ptr zstr
Dim As String n_str, p_str
Dim As Mpz_ptr a, n, p, p2, w2, x1, x2
a = Allocate(Len(__Mpz_struct)) : Mpz_init(a)
n = Allocate(Len(__Mpz_struct)) : Mpz_init(n)
p = Allocate(Len(__Mpz_struct)) : Mpz_init(p)
p2 = Allocate(Len(__Mpz_struct)) : Mpz_init(p2)
w2 = Allocate(Len(__Mpz_struct)) : Mpz_init(w2)
x1 = Allocate(Len(__Mpz_struct)) : Mpz_init(x1)
x2 = Allocate(Len(__Mpz_struct)) : Mpz_init(x2)
Dim As fp2 answer
answer.x = Allocate(Len(__Mpz_struct)) : Mpz_init(answer.x)
answer.y = Allocate(Len(__Mpz_struct)) : Mpz_init(answer.y)
For i = 1 To 8
Read n_str
Mpz_set_str(n, n_str, 10)
If i < 8 Then
Read p_str
Mpz_set_str(p, p_str, 10)
Else
p_str = "10^50 + 151" ' set up last n
Mpz_set_str(p, "1" + String(50, "0"), 10)
Mpz_add_ui(p, p, 151)
End If
Print "Find solution for n = "; n_str; " and p = "; p_str
If Mpz_tstbit(p, 0) = 0 OrElse Mpz_probab_prime_p(p, 20) = 0 Then
Print p_str; "is not a odd prime"
Print
Continue For
End If
' p is checked and is a odd prime
' legendre symbol needs to be 1
If Mpz_legendre(n, p) <> 1 Then
Print n_str; " is not a square in F"; p_str
Print
Continue For
End If
Mpz_set_ui(a, 1)
Do
Do
Do
Mpz_add_ui(a, a, 1)
Mpz_mul(w2, a, a)
Mpz_sub(w2, w2, n)
Loop Until Mpz_legendre(w2, p) = -1
Mpz_set(answer.x, a)
Mpz_set_ui(answer.y, 1)
Mpz_add_ui(p2, p, 1) ' p2 = p + 1
Mpz_fdiv_q_2exp(p2, p2, 1) ' p2 = p2 \ 2 (p2 shr 1)
answer = fp2pow(answer, p2, p, w2)
Loop Until Mpz_cmp_ui(answer.y, 0) = 0
Mpz_set(x1, answer.x)
Mpz_sub(x2, p, x1)
Mpz_powm_ui(a, x1, 2, p)
Mpz_powm_ui(p2, x2, 2, p)
If Mpz_cmp(a, n) = 0 AndAlso Mpz_cmp(p2, n) = 0 Then Exit Do
Loop
zstr = Mpz_get_str(0, 10, x1)
Print "Solution found: x1 = "; *zstr;
zstr = Mpz_get_str(0, 10, x2)
Print ", x2 = "; *zstr
Print
Next
Mpz_clear(x1) : Mpz_clear(p2) : Mpz_clear(p) : Mpz_clear(a) : Mpz_clear(n)
Mpz_clear(x2) : Mpz_clear(w2) : Mpz_clear(answer.x) : Mpz_clear(answer.y)
' empty keyboard buffer
While Inkey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End</syntaxhighlight>
{{out}}
<pre>Find solution for n = 10 and p = 13
Solution found: x1 = 6, x2 = 7
Find solution for n = 56 and p = 101
Solution found: x1 = 37, x2 = 64
Find solution for n = 8218 and p = 10007
Solution found: x1 = 9872, x2 = 135
Find solution for n = 8219 and p = 10007
8219 is not a square in F10007
Find solution for n = 331575 and p = 1000003
Solution found: x1 = 855842, x2 = 144161
Find solution for n = 665165880 and p = 1000000007
Solution found: x1 = 524868305, x2 = 475131702
Find solution for n = 881398088036 and p = 1000000000039
Solution found: x1 = 208600591990, x2 = 791399408049
Find solution for n = 34035243914635549601583369544560650254325084643201 and p = 10^50 + 151
Solution found: x1 = 17436881171909637738621006042549786426312886309400, x2 = 82563118828090362261378993957450213573687113690751</pre>
=={{header|Go}}==
===int===
Implementation following the pseudocode in the task description.
<
import "fmt"
Line 216 ⟶ 1,149:
fmt.Println(c(8219, 10007))
fmt.Println(c(331575, 1000003))
}</
{{out}}
<pre>
Line 227 ⟶ 1,160:
===big.Int===
Extra credit:
<
import (
Line 286 ⟶ 1,219:
fmt.Println(&R1)
fmt.Println(&R2)
}</
{{out}}
<pre>
Line 294 ⟶ 1,227:
17436881171909637738621006042549786426312886309400
</pre>
=={{header|J}}==
Based on the echolisp implementation:
<
x (y&|)@^ (y-1)%2
)
Line 332 ⟶ 1,264:
assert. 'not a valid square root'
end.
)</
Task examples:
<
6 7
56 cipolla 101
Line 352 ⟶ 1,284:
208600591990 791399408049
34035243914635549601583369544560650254325084643201x cipolla (10^50x) + 151
17436881171909637738621006042549786426312886309400 82563118828090362261378993957450213573687113690751</
=={{header|Java}}==
{{trans|Kotlin}}
{{works with|Java|8}}
<syntaxhighlight lang="java">import java.math.BigInteger;
import java.util.function.BiFunction;
import java.util.function.Function;
public class CipollasAlgorithm {
private static final BigInteger BIG = BigInteger.TEN.pow(50).add(BigInteger.valueOf(151));
private static final BigInteger BIG_TWO = BigInteger.valueOf(2);
private static class Point {
BigInteger x;
BigInteger y;
Point(BigInteger x, BigInteger y) {
this.x = x;
this.y = y;
}
@Override
public String toString() {
return String.format("(%s, %s)", this.x, this.y);
}
}
private static class Triple {
BigInteger x;
BigInteger y;
boolean b;
Triple(BigInteger x, BigInteger y, boolean b) {
this.x = x;
this.y = y;
this.b = b;
}
@Override
public String toString() {
return String.format("(%s, %s, %s)", this.x, this.y, this.b);
}
}
private static Triple c(String ns, String ps) {
BigInteger n = new BigInteger(ns);
BigInteger p = !ps.isEmpty() ? new BigInteger(ps) : BIG;
// Legendre symbol, returns 1, 0 or p - 1
Function<BigInteger, BigInteger> ls = (BigInteger a)
-> a.modPow(p.subtract(BigInteger.ONE).divide(BIG_TWO), p);
// Step 0, validate arguments
if (!ls.apply(n).equals(BigInteger.ONE)) {
return new Triple(BigInteger.ZERO, BigInteger.ZERO, false);
}
// Step 1, find a, omega2
BigInteger a = BigInteger.ZERO;
BigInteger omega2;
while (true) {
omega2 = a.multiply(a).add(p).subtract(n).mod(p);
if (ls.apply(omega2).equals(p.subtract(BigInteger.ONE))) {
break;
}
a = a.add(BigInteger.ONE);
}
// multiplication in Fp2
BigInteger finalOmega = omega2;
BiFunction<Point, Point, Point> mul = (Point aa, Point bb) -> new Point(
aa.x.multiply(bb.x).add(aa.y.multiply(bb.y).multiply(finalOmega)).mod(p),
aa.x.multiply(bb.y).add(bb.x.multiply(aa.y)).mod(p)
);
// Step 2, compute power
Point r = new Point(BigInteger.ONE, BigInteger.ZERO);
Point s = new Point(a, BigInteger.ONE);
BigInteger nn = p.add(BigInteger.ONE).shiftRight(1).mod(p);
while (nn.compareTo(BigInteger.ZERO) > 0) {
if (nn.and(BigInteger.ONE).equals(BigInteger.ONE)) {
r = mul.apply(r, s);
}
s = mul.apply(s, s);
nn = nn.shiftRight(1);
}
// Step 3, check x in Fp
if (!r.y.equals(BigInteger.ZERO)) {
return new Triple(BigInteger.ZERO, BigInteger.ZERO, false);
}
// Step 5, check x * x = n
if (!r.x.multiply(r.x).mod(p).equals(n)) {
return new Triple(BigInteger.ZERO, BigInteger.ZERO, false);
}
// Step 4, solutions
return new Triple(r.x, p.subtract(r.x), true);
}
public static void main(String[] args) {
System.out.println(c("10", "13"));
System.out.println(c("56", "101"));
System.out.println(c("8218", "10007"));
System.out.println(c("8219", "10007"));
System.out.println(c("331575", "1000003"));
System.out.println(c("665165880", "1000000007"));
System.out.println(c("881398088036", "1000000000039"));
System.out.println(c("34035243914635549601583369544560650254325084643201", ""));
}
}</syntaxhighlight>
{{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>
=={{header|jq}}==
{{trans|Wren}}
'''Works with gojq, the Go implementation of jq'''
A Point is represented by a numeric array of length two (i.e. [x,y]).
<syntaxhighlight lang="jq"># To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
# If $j is 0, then an error condition is raised;
# otherwise, assuming infinite-precision integer arithmetic,
# if the input and $j are integers, then the result will be an integer.
def idivide($j):
. as $i
| ($i % $j) as $mod
| ($i - $mod) / $j ;
def bigBig: (10|power(50)) + 151;
# The remainder when $b raised to the power $e is divided by $m:
def modPow($b; $e; $m):
if ($m == 1) then 0
else {r: 1, b: ($b % $m), $e}
| until (.e <= 0;
if (.e % 2) == 1 then .r = (.r*.b) % $m else . end
| .e |= idivide(2)
| .b |= ((.*.) % $m) )
| .r
end;
def c($ns; $ps):
$ns as $n
| (if $ps != "" then $ps else bigBig end) as $p
# Legendre symbol, returns 1, 0 or p - 1
| def ls($a): modPow($a; ($p - 1) | idivide(2); $p);
# multiplication in Fp2 where .omega2 comes from .
def mul($aa; $bb):
[($aa[0] * $bb[0] + $aa[1] * $bb[1] * .omega2) % $p,
($aa[0] * $bb[1] + $bb[0] * $aa[1]) % $p] ;
# Step 0, validate arguments
if (ls($n) != 1) then [0, 0, false]
else
# Step 1, find a, omega2
{ a: 0, stop: false }
| until( .stop;
.omega2 = ((.a * .a + $p - $n) % $p)
| if ls(.omega2) == ($p - 1)
then .stop = true
else .a += 1
end )
# Step 2, compute power
| { r: [1, 0],
s: [.a, 1],
nn: (( ($p + 1) | idivide(2)) % $p),
omega2 }
| until (.nn <= 0;
if (.nn % 2) == 1 then .r = mul(.r; .s) else . end
| .s = mul(.s; .s)
| .nn |= idivide(2) )
# Step 3, check x in Fp; or that x * x = n
| if (.r[1] != 0) or ((.r[0] * .r[0]) % $p != $n) then [0, 0, false]
# Step 4, solutions
else [.r[0], $p - .r[0], true]
end
end;
def exercise:
c(10; 13),
c(56; 101),
c(8218; 10007),
c(8219; 10007),
c(331575; 1000003),
c(665165880; 1000000007),
c(881398088036; 1000000000039),
c(34035243914635549601583369544560650254325084643201; "")
;
exercise</syntaxhighlight>
{{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>
=={{header|Julia}}==
{{trans|Perl}}
<syntaxhighlight lang="julia">using Primes
function legendre(n, p)
if p != 2 && isprime(p)
x = powermod(BigInt(n), div(p - 1, 2), p)
return x == 0 ? 0 : x == 1 ? 1 : -1
end
return -1
end
function cipolla(n, p)
if legendre(n, p) != 1
return NaN
end
a, w2 = BigInt(0), BigInt(0)
while true
w2 = (a^2 + p - n) % p
if legendre(w2, p) < 0
break
end
a += 1
end
r, s, i = (1, 0), (a, 1), p + 1
while (i >>= 1) > 0
if isodd(i)
r = ((r[1] * s[1] + r[2] * s[2] * w2) % p, (r[1] * s[2] + s[1] * r[2]) % p)
end
s = ((s[1] * s[1] + s[2] * s[2] * w2) % p, (2 * s[1] * s[2]) % p)
end
return r[2] != 0 ? NaN : r[1]
end
const ctests = [(10, 13),
(56, 101),
(8218, 10007),
(8219, 10007),
(331575, 1000003),
(665165880, 1000000007),
(881398088036, 1000000000039),
(big"34035243914635549601583369544560650254325084643201",
big"100000000000000000000000000000000000000000000000151")]
for (n, p) in ctests
r = cipolla(n, p)
println(r > 0 ? "Roots of $n are ($r, $(p - r)) mod $p." : "No solution for ($n, $p)")
end
</syntaxhighlight>{{out}}
<pre>
Roots of 10 are (6, 7) mod 13.
Roots of 56 are (37, 64) mod 101.
Roots of 8218 are (9872, 135) mod 10007.
No solution for (8219, 10007)
Roots of 331575 are (855842, 144161) mod 1000003.
Roots of 665165880 are (475131702, 524868305) mod 1000000007.
Roots of 881398088036 are (791399408049, 208600591990) mod 1000000000039.
Roots of 34035243914635549601583369544560650254325084643201 are (82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400) mod 100000000000000000000000000000000000000000000000151.
</pre>
=={{header|Kotlin}}==
{{trans|Go}}
<syntaxhighlight lang="scala">// version 1.2.0
import java.math.BigInteger
class Point(val x: BigInteger, val y: BigInteger)
val bigZero = BigInteger.ZERO
val bigOne = BigInteger.ONE
val bigTwo = BigInteger.valueOf(2L)
val bigBig = BigInteger.TEN.pow(50) + BigInteger.valueOf(151L)
fun c(ns: String, ps: String): Triple<BigInteger, BigInteger, Boolean> {
val n = BigInteger(ns)
val p = if (!ps.isEmpty()) BigInteger(ps) else bigBig
// Legendre symbol, returns 1, 0 or p - 1
fun ls(a: BigInteger) = a.modPow((p - bigOne) / bigTwo, p)
// Step 0, validate arguments
if (ls(n) != bigOne) return Triple(bigZero, bigZero, false)
// Step 1, find a, omega2
var a = bigZero
var omega2: BigInteger
while (true) {
omega2 = (a * a + p - n) % p
if (ls(omega2) == p - bigOne) break
a++
}
// multiplication in Fp2
fun mul(aa: Point, bb: Point) =
Point(
(aa.x * bb.x + aa.y * bb.y * omega2) % p,
(aa.x * bb.y + bb.x * aa.y) % p
)
// Step 2, compute power
var r = Point(bigOne, bigZero)
var s = Point(a, bigOne)
var nn = ((p + bigOne) shr 1) % p
while (nn > bigZero) {
if ((nn and bigOne) == bigOne) r = mul(r, s)
s = mul(s, s)
nn = nn shr 1
}
// Step 3, check x in Fp
if (r.y != bigZero) return Triple(bigZero, bigZero, false)
// Step 5, check x * x = n
if (r.x * r.x % p != n) return Triple(bigZero, bigZero, false)
// Step 4, solutions
return Triple(r.x, p - r.x, true)
}
fun main(args: Array<String>) {
println(c("10", "13"))
println(c("56", "101"))
println(c("8218", "10007"))
println(c("8219", "10007"))
println(c("331575", "1000003"))
println(c("665165880", "1000000007"))
println(c("881398088036", "1000000000039"))
println(c("34035243914635549601583369544560650254325084643201", ""))
}</syntaxhighlight>
{{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>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight 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]</syntaxhighlight>
{{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>
=={{header|Nim}}==
{{trans|Kotlin}}
{{libheader|bignum}}
<syntaxhighlight lang="nim">import options
import bignum
let
Zero = newInt(0)
One = newInt(1)
BigBig = newInt(10)^50 + 15
type
Point = tuple[x, y: Int]
Solutions = (Int, Int)
proc exp(x, y, m: Int): Int =
## Missing function in "bignum" module.
if m == 1: return Zero
result = newInt(1)
var x = x mod m
var y = y.clone
while not y.isZero:
if y mod 2 == 1:
result = result * x mod m
y = y shr 1
x = x * x mod m
proc c(ns, ps: string): Option[Solutions] =
let n = newInt(ns)
let p = if ps.len != 0: newInt(ps) else: BigBig
# Legendre symbol: returns 1, 0 or p - 1.
proc ls(a: Int): Int = a.exp((p - 1) div 2, p)
# Step 0, validate arguments.
if ls(n) != One: return none(Solutions)
# Step 1, find a, omega2.
var a = newInt(0)
var omega2: Int
while true:
omega2 = (a * a + p - n) mod p
if ls(omega2) == p - 1: break
a += 1
# Multiplication in Fp2.
proc `*`(a, b: Point): Point =
((a.x * b.x + a.y * b.y * omega2) mod p,
(a.x * b.y + b.x * a.y) mod p)
# Step 2, compute power.
var
r: Point = (One, Zero)
s: Point = (a, One)
nn = ((p + 1) shr 1) mod p
while not nn.isZero:
if (nn and One) == One: r = r * s
s = s * s
nn = nn shr 1
# Step 3, check x in Fp.
if not r.y.isZero: return none(Solutions)
# Step 5, check x * x = n.
if r.x * r.x mod p != n: return none(Solutions)
# Step 4, solutions.
result = some((r.x, p - r.x))
when isMainModule:
const Values = [("10", "13"), ("56", "101"), ("8218", "10007"),
("8219", "10007"), ("331575", "1000003"),
("665165880", "1000000007"), ("881398088036", "1000000000039"),
("34035243914635549601583369544560650254325084643201", "")]
for (n, p) in Values:
let sols = c(n, p)
if sols.isSome: echo sols.get()
else: echo "No solutions."</syntaxhighlight>
{{out}}
<pre>(6, 7)
(37, 64)
(9872, 135)
No solutions.
(855842, 144161)
(475131702, 524868305)
(791399408049, 208600591990)
(82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400)</pre>
=={{header|Perl}}==
{{trans|Raku}}
{{libheader|ntheory}}
<syntaxhighlight lang="perl">use bigint;
use ntheory qw(is_prime);
sub Legendre {
my($n,$p) = @_;
return -1 unless $p != 2 && is_prime($p);
my $x = ($n->as_int())->bmodpow(int(($p-1)/2), $p); # $n coerced to BigInt
if ($x==0) { return 0 }
elsif ($x==1) { return 1 }
else { return -1 }
}
sub Cipolla {
my($n, $p) = @_;
return undef if Legendre($n,$p) != 1;
my $w2;
my $a = 0;
$a++ until Legendre(($w2 = ($a**2 - $n) % $p), $p) < 0;
my %r = ( x=> 1, y=> 0 );
my %s = ( x=> $a, y=> 1 );
my $i = $p + 1;
while (1 <= ($i >>= 1)) {
%r = ( x => (($r{x} * $s{x} + $r{y} * $s{y} * $w2) % $p),
y => (($r{x} * $s{y} + $s{x} * $r{y}) % $p)
) if $i % 2;
%s = ( x => (($s{x} * $s{x} + $s{y} * $s{y} * $w2) % $p),
y => (($s{x} * $s{y} + $s{x} * $s{y}) % $p)
)
}
$r{y} ? undef : $r{x}
}
my @tests = (
(10, 13),
(56, 101),
(8218, 10007),
(8219, 10007),
(331575, 1000003),
(665165880, 1000000007),
(881398088036, 1000000000039),
);
while (@tests) {
$n = shift @tests;
$p = shift @tests;
my $r = Cipolla($n, $p);
$r ? printf "Roots of %d are (%d, %d) mod %d\n", $n, $r, $p-$r, $p
: print "No solution for ($n, $p)\n"
}</syntaxhighlight>
{{out}}
<pre>Roots of 10 are (6, 7) mod 13
Roots of 56 are (37, 64) mod 101
Roots of 8218 are (9872, 135) mod 10007
No solution for (8219, 10007)
Roots of 331575 are (855842, 144161) mod 1000003
Roots of 665165880 are (475131702, 524868305) mod 1000000007
Roots of 881398088036 are (791399408049, 208600591990) mod 1000000000039</pre>
=={{header|Phix}}==
{{trans|Kotlin}}
{{libheader|Phix/mpfr}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpz</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- Legendre symbol, returns 1, 0 or p - 1 (in r)</span>
<span style="color: #7060A8;">mpz_sub_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_fdiv_q_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_powm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">mul_point</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">mpz</span> <span style="color: #000000;">omega2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- (modifies a)</span>
<span style="color: #004080;">mpz</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">xa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ya</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">xb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">yb</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">xaxb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">yayb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">xayb</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">xbya</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xaxb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xb</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">yayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ya</span><span style="color: #0000FF;">,</span><span style="color: #000000;">yb</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">yb</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xbya</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ya</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">yayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">yayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">omega2</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xaxb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xaxb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">yayb</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xa</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xaxb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- xa := mod(xaxb+yayb*omega2,p)</span>
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xbya</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ya</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- ya := mod(xayb+xbya,p)</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">xaxb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">yayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xbya</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_free</span><span style="color: #0000FF;">({</span><span style="color: #000000;">xaxb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">yayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xayb</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xbya</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">cipolla</span><span style="color: #0000FF;">(</span><span style="color: #004080;">object</span> <span style="color: #000000;">no</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">po</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">no</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">po</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
<span style="color: #000080;font-style:italic;">-- Step 0, validate arguments</span>
<span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"false"</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000080;font-style:italic;">-- Step 1, find a, omega2</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">omega2</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">pm1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">mpz_sub_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pm1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">mpz_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_add_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">a</span><span style="color: #0000FF;">*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">omega2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">legendre</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">omega2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pm1</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">a</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #000080;font-style:italic;">-- Step 2, compute power</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">r</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)},</span>
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)}</span>
<span style="color: #004080;">mpz</span> <span style="color: #000000;">nn</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">mpz_add_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_fdiv_q_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpz_mod</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">while</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)></span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_fdiv_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">mul_point</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #000000;">omega2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">mul_point</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #000000;">omega2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_fdiv_q_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nn</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #000080;font-style:italic;">-- Step 3, check x in Fp</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">],</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"false"</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000080;font-style:italic;">-- Step 5, check x * x = n</span>
<span style="color: #7060A8;">mpz_powm_ui</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpz_cmp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)!=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"0"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"false"</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000080;font-style:italic;">-- Step 4, solutions</span>
<span style="color: #7060A8;">mpz_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]),</span> <span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">),</span> <span style="color: #008000;">"true"</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">tests</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">13</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">56</span><span style="color: #0000FF;">,</span><span style="color: #000000;">101</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">8218</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10007</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">8219</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10007</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">331575</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1000003</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">665165880</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1000000007</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #008000;">"881398088036"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"1000000000039"</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #008000;">"34035243914635549601583369544560650254325084643201"</span><span style="color: #0000FF;">,</span>
<span style="color: #008000;">"100000000000000000000000000000000000000000000000151"</span><span style="color: #0000FF;">}}</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tests</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">object</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tests</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #0000FF;">?{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cipolla</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
Obviously were you to use that in anger, you would probably rip out a few mpz_get_str() and return NULL/false rather than "0"/"false", etc.
{{out}}
<pre>
{10,13,{"6","7","true"}}
{56,101,{"37","64","true"}}
{8218,10007,{"9872","135","true"}}
{8219,10007,{"0","0","false"}}
{331575,1000003,{"855842","144161","true"}}
{665165880,1000000007,{"475131702","524868305","true"}}
{"881398088036","1000000000039",{"791399408049","208600591990","true"}}
{"34035243914635549601583369544560650254325084643201","100000000000000000000000000000000000000000000000151",
{"82563118828090362261378993957450213573687113690751","17436881171909637738621006042549786426312886309400","true"}}
</pre>
=={{header|PicoLisp}}==
{{trans|Go}}
<syntaxhighlight lang="picolisp"># from @lib/rsa.l
(de **Mod (X Y N)
(let M 1
(loop
(when (bit? 1 Y)
(setq M (% (* M X) N)) )
(T (=0 (setq Y (>> 1 Y)))
M )
(setq X (% (* X X) N)) ) ) )
(de legendre (N P)
(**Mod N (/ (dec P) 2) P) )
(de mul ("A" B P W2)
(let (A (copy "A") B (copy B))
(set
"A"
(%
(+
(* (car A) (car B))
(* (cadr A) (cadr B) W2) )
P )
(cdr "A")
(%
(+
(* (car A) (cadr B))
(* (car B) (cadr A)) )
P ) ) ) )
(de ci (N P)
(and
(=1 (legendre N P))
(let
(A 0
W2 0
R NIL
S NIL )
(loop
(setq W2
(% (- (+ (* A A) P) N) P) )
(T (= (dec P) (legendre W2 P)))
(inc 'A) )
(setq R (list 1 0) S (list A 1))
(for
(N
(% (>> 1 (inc P)) P)
(> N 0)
(>> 1 N) )
(and (bit? 1 N) (mul R S P W2))
(mul S S P W2) )
(=0 (cadr R))
(=
N
(% (* (car R) (car R)) P) )
(list (car R) (- P (car R))) ) ) )
(println (ci 10 13))
(println (ci 56 101))
(println (ci 8218 10007))
(println (ci 8219 10007))
(println (ci 331575 1000003))
(println (ci 665165880 1000000007))
(println (ci 881398088036 1000000000039))
(println (ci 34035243914635549601583369544560650254325084643201 (+ (** 10 50) 151)))</syntaxhighlight>
{{out}}
<pre>
(6 7)
(37 64)
(9872 135)
NIL
(855842 144161)
(475131702 524868305)
(791399408049 208600591990)
(82563118828090362261378993957450213573687113690751 17436881171909637738621006042549786426312886309400)
</pre>
=={{header|Python}}==
<syntaxhighlight lang="python">
#Converts n to base b as a list of integers between 0 and b-1
#Most-significant digit on the left
def convertToBase(n, b):
if(n < 2):
return [n];
temp = n;
ans = [];
while(temp != 0):
ans = [temp % b]+ ans;
temp /= b;
return ans;
#Takes integer n and odd prime p
#Returns both square roots of n modulo p as a pair (a,b)
#Returns () if no root
def cipolla(n,p):
n %= p
if(n == 0 or n == 1):
return (n,-n%p)
phi = p - 1
if(pow(n, phi/2, p) != 1):
return ()
if(p%4 == 3):
ans = pow(n,(p+1)/4,p)
return (ans,-ans%p)
aa = 0
for i in xrange(1,p):
temp = pow((i*i-n)%p,phi/2,p)
if(temp == phi):
aa = i
break;
exponent = convertToBase((p+1)/2,2)
def cipollaMult((a,b),(c,d),w,p):
return ((a*c+b*d*w)%p,(a*d+b*c)%p)
x1 = (aa,1)
x2 = cipollaMult(x1,x1,aa*aa-n,p)
for i in xrange(1,len(exponent)):
if(exponent[i] == 0):
x2 = cipollaMult(x2,x1,aa*aa-n,p)
x1 = cipollaMult(x1,x1,aa*aa-n,p)
else:
x1 = cipollaMult(x1,x2,aa*aa-n,p)
x2 = cipollaMult(x2,x2,aa*aa-n,p)
return (x1[0],-x1[0]%p)
print(f"Roots of 2 mod 7: {cipolla(2, 7)}")
print(f"Roots of 8218 mod 10007: {cipolla(8218, 10007)}")
print(f"Roots of 56 mod 101: {cipolla(56, 101)}")
print(f"Roots of 1 mod 11: {cipolla(1, 11)}")
print(f"Roots of 8219 mod 10007: {cipolla(8219, 10007)}")
</syntaxhighlight>
{{out}}
<pre>
Roots of 2 mod 7: (4, 3)
Roots of 8218 mod 10007: (9872, 135)
Roots of 56 mod 101: (37, 64)
Roots of 1 mod 11: (1, 10)
Roots of 8219 mod 10007: ()
</pre>
=={{header|Racket}}==
{{trans|EchoLisp}}
<syntaxhighlight lang="racket">#lang racket
(require math/number-theory)
;; math/number-theory allows us to parameterize a "current-modulus"
;; which obviates the need for p to be passed around constantly
(define (Cipolla n p) (with-modulus p (mod-Cipolla n)))
(define (mod-Legendre a)
(modexpt a (/ (sub1 (current-modulus)) 2)))
;; Arithmetic in Fp²
(struct Fp² (x y))
(define-syntax-rule (Fp²-destruct* (a a.x a.y) ...)
(begin (match-define (Fp² a.x a.y) a) ...) )
;; a + b
(define (Fp²-add a b ω2)
(Fp²-destruct* (a a.x a.y) (b b.x b.y))
(Fp² (mod+ a.x b.x) (mod+ a.y b.y)))
;; a * b
(define (Fp²-mul a b ω2)
(Fp²-destruct* (a a.x a.y) (b b.x b.y))
(Fp² (mod+ (* a.x b.x) (* ω2 a.y b.y)) (mod+ (* a.x b.y) (* a.y b.x))))
;; a * a
(define (Fp²-square a ω2)
(Fp²-destruct* (a a.x a.y))
(Fp² (mod+ (sqr a.x) (* ω2 (sqr a.y))) (mod* 2 a.x a.y)))
;; a ^ n
(define (Fp²-pow a n ω2)
(Fp²-destruct* (a a.x a.y))
(cond
((= 0 n) (Fp² 1 0))
((= 1 n) a)
((= 2 n) (Fp²-mul a a ω2))
((even? n) (Fp²-square (Fp²-pow a (/ n 2) ω2) ω2))
(else (Fp²-mul a (Fp²-pow a (sub1 n) ω2) ω2))))
;; x^2 ≡ n (mod p) ?
(define (mod-Cipolla n)
;; check n is a square
(unless (= 1 (mod-Legendre n)) (error 'Cipolla "~a not a square (mod ~a)" n (current-modulus)))
;; iterate until suitable 'a' found
(define a (for/first ((t (in-range 2 (current-modulus))) ;; t = tentative a
#:when (= (sub1 (current-modulus))
(mod-Legendre (- (* t t) n))))
t))
(define ω2 (- (* a a) n))
;; (Fp² a 1) = a + ω
(define r (Fp²-pow (Fp² a 1) (/ (add1 (current-modulus)) 2) ω2))
(define x (Fp²-x r))
(unless (zero? (Fp²-y r)) (error 'Cipolla "ω has not vanished")) ;; hope that ω has vanished
(unless (mod= n (* x x)) (error 'Cipolla "result check failed")) ;; checking the result
(values x (mod- (current-modulus) x)))
(define (report-Cipolla n p)
(with-handlers ((exn:fail? (λ (x) (eprintf "Caught error: ~s~%" (exn-message x)))))
(define-values (r1 r2) (Cipolla n p))
(printf "Roots of ~a are (~a,~a) (mod ~a)~%" n r1 r2 p)))
(module+ test
(report-Cipolla 10 13)
(report-Cipolla 56 101)
(report-Cipolla 8218 10007)
(report-Cipolla 8219 10007)
(report-Cipolla 331575 1000003)
(report-Cipolla 665165880 1000000007)
(report-Cipolla 881398088036 1000000000039)
(report-Cipolla 34035243914635549601583369544560650254325084643201
100000000000000000000000000000000000000000000000151))</syntaxhighlight>
{{out}}
<pre>Roots of 10 are (6,7) (mod 13)
Roots of 56 are (37,64) (mod 101)
Roots of 8218 are (9872,135) (mod 10007)
Caught error: "Cipolla: 8219 not a square (mod 10007)"
Roots of 331575 are (855842,144161) (mod 1000003)
Roots of 665165880 are (524868305,475131702) (mod 1000000007)
Roots of 881398088036 are (208600591990,791399408049) (mod 1000000000039)
Roots of 34035243914635549601583369544560650254325084643201 are (17436881171909637738621006042549786426312886309400,82563118828090362261378993957450213573687113690751) (mod 100000000000000000000000000000000000000000000000151)</pre>
=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2016.10}}
{{trans|Sidef}}
<syntaxhighlight lang="raku"
sub infix:<│> (Int \𝑛, Int \𝑝 where
given 𝑛.expmod( (𝑝-1) div 2, 𝑝 ) {
when 0 { 0 }
Line 376 ⟶ 2,203:
note "Invalid parameters ({𝑛}, {𝑝})"
and return Nil if (𝑛│𝑝) != 1;
my
loop
last if ($ω2 = ($a² - 𝑛) % 𝑝)│𝑝 < 0;
}
# define a local multiply operator for Field coordinates
multi sub infix:<*> ( Fp $a, Fp $b ){
Fp.new: :x(($a.x * $b.x + $a.y * $b.y * $ω2) % 𝑝),
:
}
my $r = Fp.new
my $s = Fp.new
for (𝑝+1) +> 1, * +> 1 ...
$r *= $s if $_ % 2;
$s *= $s;
Line 416 ⟶ 2,241:
for @tests -> ($n, $p) {
my $r = cipolla($n, $p);
say $r ?? "Roots of $n are ($r, {$p-$r}) mod $p"
}
</syntaxhighlight>
{{out}}
<pre>Roots of 10 are (6, 7) mod 13
Line 433 ⟶ 2,256:
Roots of 34035243914635549601583369544560650254325084643201 are (82563118828090362261378993957450213573687113690751, 17436881171909637738621006042549786426312886309400) mod 100000000000000000000000000000000000000000000000151
</pre>
=={{header|Sage}}==
{{works with|Sage|7.
<syntaxhighlight lang="sage">
def eulerCriterion(a, p):
return -1 if pow(a, int((p-1)/2), p) == p-1 else 1
def cipollaMult(x1, y1, x2, y2, u, p):
return ((x1*x2 + y1*y2*u) % p), ((x1*y2 + x2*y1) % p)
def cipollaAlgorithm(n, p):
a =
out = []
if eulerCriterion(a, p) == -1:
print "❌ " + str(a) + " is not a quadratic residue modulo " + str(p)
return False
conglst = [] #congruence list
for k in list(factor(p)):
conglst.append(cipollaAlgorithm(a, f))
for i in Permutations([0, 1] * len(factors), len(factors)).list():
for j in range(len(factors)):
crtlst.append(int(conglst[ j ][ i[j] ]))
out.append(crt(crtlst, factors))
crtlst = []
return sorted(out)
if pow(p, 1, 4) == 3:
temp = pow(a, int((p+1)/4), p)
return [temp, p - temp]
t = randrange(2, p)
u = pow(t**2 - a, 1, p)
while (eulerCriterion(u, p) == 1):
t = randrange(2, p)
u = pow(t**2 - a, 1, p)
x0, y0 = t, 1
x, y = t, 1
for i in range(int((p + 1) / 2) - 1):
x, y = cipollaMult(x, y, x0, y0, u, p)
out.extend([x, p - x])
return sorted(out)
</syntaxhighlight>
{{out}}
<pre>
sage:
sage:
sage:
[135, 9872]
sage:
[144161, 855842]
sage:
❌ 8219 is not a
False
</pre>
=={{header|Scala}}==
===Imperative solution===
<syntaxhighlight lang="scala">object CipollasAlgorithm extends App {
private val BIG = BigInt(10).pow(50) + BigInt(151)
println(c("10", "13"))
println(c("56", "101"))
println(c("8218", "10007"))
println(c("8219", "10007"))
println(c("331575", "1000003"))
println(c("665165880", "1000000007"))
println(c("881398088036", "1000000000039"))
println(c("34035243914635549601583369544560650254325084643201", ""))
private def c(ns: String, ps: String): Triple = {
val (n, p) = (BigInt(ns), if (ps.isEmpty) BIG else BigInt(ps))
// Legendre symbol, returns 1, 0 or p - 1
def ls(a: BigInt) = a.modPow((p - BigInt(1)) / BigInt(2), p)
// multiplication in Fp2
def mul(aa: Point, bb: Point, omega2: BigInt) =
new Point((aa.x * bb.x + aa.y * bb.y * omega2) % p, (aa.x * bb.y + (bb.x * aa.y)) % p)
// Step 0, validate arguments
if (ls(n) != BigInt(1)) new Triple(0, 0, false)
else {
// Step 1, find a, omega2
var (a, flag, omega2) = (BigInt(0), true, BigInt(0))
while (flag) {
omega2 = (a * a + p - n) % p
if (ls(omega2) == p - BigInt(1)) flag = false else a = a + BigInt(1)
}
// Step 2, compute power
var (nn, r, s) = ((p + BigInt(1) >> 1) % p, new Point(BigInt(1), 0), new Point(a, BigInt(1)))
while (nn > 0) {
if ((nn & BigInt(1)) == BigInt(1)) r = mul(r, s, omega2)
s = mul(s, s, omega2)
nn = nn >> 1
}
// Step 3, check x in Fp
if (r.y != 0) new Triple(0, 0, false)
else // Step 5, check x * x = n
if ((r.x * r.x) % p != n) new Triple(0, 0, false)
else new Triple(r.x, p - r.x, true) // Step 4, solutions
}
}
private class Point(val x: BigInt, val y: BigInt)
private class Triple(val x: BigInt, val y: BigInt, val b: Boolean) {
override def toString: String = f"($x%s, $y%s, $b%s)"
}
}</syntaxhighlight>
{{Out}}See it running in your browser by [https://scalafiddle.io/sf/QQBsMza/3 ScalaFiddle (JavaScript, non JVM)] or by [https://scastie.scala-lang.org/NEP5hOWmSBqqpwmF30LpUA Scastie (JVM)].
=={{header|Sidef}}==
{{trans|Go}}
<
legendre(n, p) == 1 || return nil
Line 520 ⟶ 2,434:
say "No solution for (#{n}, #{p})"
}
}</
{{out}}
<pre>Roots of 10 are (6 7) mod 13
Roots of 56 are (37 64) mod 101
Roots of 8218 are (9872 135) mod 10007
Line 530 ⟶ 2,443:
Roots of 665165880 are (475131702 524868305) mod 1000000007
Roots of 881398088036 are (791399408049 208600591990) mod 1000000000039
Roots of 34035243914635549601583369544560650254325084643201 are (82563118828090362261378993957450213573687113690751 17436881171909637738621006042549786426312886309400) mod 100000000000000000000000000000000000000000000000151</pre>
Simpler implementation:
<syntaxhighlight lang="ruby">func cipolla(n, p) {
kronecker(n, p) == 1 || return nil
var (a, ω) = (
0..Inf -> lazy.map {|a|
[a, submod(a*a, n, p)]
}.first_by {|t|
kronecker(t[1], p) == -1
}...
)
var r = lift(Mod(Quadratic(a, 1, ω), p)**((p+1)>>1))
r.b == 0 ? r.a : nil
}</syntaxhighlight>
=={{header|Visual Basic .NET}}==
{{trans|C#}}
<syntaxhighlight lang="vbnet">Imports System.Numerics
Module Module1
ReadOnly BIG = BigInteger.Pow(10, 50) + 151
Function C(ns As String, ps As String) As Tuple(Of BigInteger, BigInteger, Boolean)
Dim n = BigInteger.Parse(ns)
Dim p = If(ps.Length > 0, BigInteger.Parse(ps), BIG)
' Legendre symbol. Returns 1, 0, or p-1
Dim ls = Function(a0 As BigInteger) BigInteger.ModPow(a0, (p - 1) / 2, p)
' Step 0: validate arguments
If ls(n) <> 1 Then
Return Tuple.Create(BigInteger.Zero, BigInteger.Zero, False)
End If
' Step 1: Find a, omega2
Dim a = BigInteger.Zero
Dim omega2 As BigInteger
Do
omega2 = (a * a + p - n) Mod p
If ls(omega2) = p - 1 Then
Exit Do
End If
a += 1
Loop
' Multiplication in Fp2
Dim mul = Function(aa As Tuple(Of BigInteger, BigInteger), bb As Tuple(Of BigInteger, BigInteger))
Return Tuple.Create((aa.Item1 * bb.Item1 + aa.Item2 * bb.Item2 * omega2) Mod p, (aa.Item1 * bb.Item2 + bb.Item1 * aa.Item2) Mod p)
End Function
' Step 2: Compute power
Dim r = Tuple.Create(BigInteger.One, BigInteger.Zero)
Dim s = Tuple.Create(a, BigInteger.One)
Dim nn = ((p + 1) >> 1) Mod p
While nn > 0
If nn Mod 2 = 1 Then
r = mul(r, s)
End If
s = mul(s, s)
nn >>= 1
End While
' Step 3: Check x in Fp
If r.Item2 <> 0 Then
Return Tuple.Create(BigInteger.Zero, BigInteger.Zero, False)
End If
' Step 5: Check x * x = n
If r.Item1 * r.Item1 Mod p <> n Then
Return Tuple.Create(BigInteger.Zero, BigInteger.Zero, False)
End If
' Step 4: Solutions
Return Tuple.Create(r.Item1, p - r.Item1, True)
End Function
Sub Main()
Console.WriteLine(C("10", "13"))
Console.WriteLine(C("56", "101"))
Console.WriteLine(C("8218", "10007"))
Console.WriteLine(C("8219", "10007"))
Console.WriteLine(C("331575", "1000003"))
Console.WriteLine(C("665165880", "1000000007"))
Console.WriteLine(C("881398088036", "1000000000039"))
Console.WriteLine(C("34035243914635549601583369544560650254325084643201", ""))
End Sub
End Module</syntaxhighlight>
{{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>
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-big}}
{{libheader|Wren-dynamic}}
<syntaxhighlight lang="wren">import "./big" for BigInt
import "./dynamic" for Tuple
var Point = Tuple.create("Point", ["x", "y"])
var bigBig = BigInt.ten.pow(50) + BigInt.new(151)
var c = Fn.new { |ns, ps|
var n = BigInt.new(ns)
var p = (ps != "") ? BigInt.new(ps) : bigBig
// Legendre symbol, returns 1, 0 or p - 1
var ls = Fn.new { |a| a.modPow((p - BigInt.one) / BigInt.two, p) }
// Step 0, validate arguments
if (ls.call(n) != BigInt.one) return [BigInt.zero, BigInt.zero, false]
// Step 1, find a, omega2
var a = BigInt.zero
var omega2
while (true) {
omega2 = (a * a + p - n) % p
if (ls.call(omega2) == p - BigInt.one) break
a = a.inc
}
// multiplication in Fp2
var mul = Fn.new { |aa, bb|
return Point.new((aa.x * bb.x + aa.y * bb.y * omega2) % p,
(aa.x * bb.y + bb.x * aa.y) % p)
}
// Step 2, compute power
var r = Point.new(BigInt.one, BigInt.zero)
var s = Point.new(a, BigInt.one)
var nn = ((p + BigInt.one) >> 1) % p
while (nn > BigInt.zero) {
if ((nn & BigInt.one) == BigInt.one) r = mul.call(r, s)
s = mul.call(s, s)
nn = nn >> 1
}
// Step 3, check x in Fp
if (r.y != BigInt.zero) return [BigInt.zero, BigInt.zero, false]
// Step 5, check x * x = n
if (r.x * r.x % p != n) return [BigInt.zero, BigInt.zero, false]
// Step 4, solutions
return [r.x, p - r.x, true]
}
System.print(c.call("10", "13"))
System.print(c.call("56", "101"))
System.print(c.call("8218", "10007"))
System.print(c.call("8219", "10007"))
System.print(c.call("331575", "1000003"))
System.print(c.call("665165880", "1000000007"))
System.print(c.call("881398088036", "1000000000039"))
System.print(c.call("34035243914635549601583369544560650254325084643201", ""))</syntaxhighlight>
{{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>
Line 536 ⟶ 2,626:
{{trans|EchoLisp}}
Uses lib GMP (GNU MP Bignum Library).
<
fcn modEq(a,b,p) { (a-b)%p==0 }
fcn Legendre(a,p){ a.powm((p - 1)/2,p) }
Line 566 ⟶ 2,656:
println("Roots of %d are (%d,%d) (mod %d)".fmt(n,x,(p-x)%p,p));
return(x,(p-x)%p);
}</
<
T(10,13),T(56,101),T(8218,10007),T(8219,10007),T(331575,1000003),
T(665165880,1000000007),T(881398088036,1000000000039),
Line 573 ⟶ 2,663:
BN(10).pow(50) + 151) )){
try{ Cipolla(n,p) }catch{ println(__exception) }
}</
{{out}}
<pre>
| |||