If you want to find a square root of a (mod p) use the ShanksTonelli algorithm.
# Calculate the least nonnegative remainder when an integer a
# is divided by a positive integer b.
sub mod {
my ($a, $b) = @_;
my $c = $a % $b;
return $c if ($a >= 0);
return 0 if ($c == 0);
return ($c + $b);
}
# Calculate a^b (mod c), where a,b and c are integers and a,b>=0, c>=1
+
sub mpow {
my ($a, $b, $c) = @_;
my ($x, $y, $z) = ($a, $b, 1);
while ($y > 0) {
while ($y % 2 == 0) {
$y = $y / 2;
$x = $x**2 % $c;
}
$y;
$z = mod($z * $x, $c);
}
return $z;
}
# ShanksTonelli algorithm to calculate y^2 = a (mod p) for p an odd p
+rime
sub tonelli {
my ($a, $p) = @_;
my ($b, $e, $g, $h, $i, $m, $n, $q, $r, $s, $t, $x, $y, $z);
$q = $p  1;
$t = mpow($a, $q/2, $p);
return 0 if ($t = $q); # a is a quadratic nonresidue mod p
$s = $q;
$e = 0;
while ($s % 2 == 0) {
$s = $s/2;
$e++;
}
# p1 = s * 2^e;
$x = mpow($a, ++$s/2, $p);
$b = mpow($a, $s, $p);
return $x if ($b == 1);
$n = 2;
RES: while ($n >= 2) {
$t = mpow($n, $q/2, $p);
last RES if ($t == $q);
} continue {
$n++;
}
# n is the least quadratic nonresidue mod p
$g = mpow($n, $s, $p);
$r = $e;
OUT: {
do {
$y = $b;
$m = 0;
IN: while ($m <= $r1) {
last IN if ($y == 1);
$y = $y**2 % $p;
} continue {
$m++;
}
return $x if ($m == 0);
$z = $r  $m  1;
$h = $g;
for ($i = 1; $i < $z; $i++) {$h = $h**2 % $p;}
$x = ($x * $h) % $p;
$b = ($b * $h**2) % $p;
$g = $h**2 % $p;
$r = $m;
} while 1;
}
}
