good chemistry is complicated,and a little bit messy -LW PerlMonks

### Re: Re: Re: Factors

by gumby (Scribe)
 on Jun 13, 2002 at 15:43 UTC ( #174242=note: print w/replies, xml ) Need Help??

in reply to Re: Re: Factors

If you want to find a square root of a (mod p) use the Shanks-Tonelli algorithm.
```# Calculate the least non-negative 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;
}

# Shanks-Tonelli 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 non-residue mod p

\$s = \$q;
\$e = 0;
while (\$s % 2 == 0) {
\$s = \$s/2;
\$e++;
}

# p-1 = 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 non-residue mod p
\$g = mpow(\$n, \$s, \$p);
\$r = \$e;

OUT:  {
do {
\$y = \$b;
\$m = 0;
IN:   while (\$m <= \$r-1) {
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;
}
}

Replies are listed 'Best First'.
Re: Re: Re: Re: Factors
by gumby (Scribe) on Jun 15, 2002 at 19:21 UTC
Here's a more 'reference' style implementation. Basically, this is the algorithm you'll find in the paper 'Square Roots, From 1; 24, 51, 10 To Dan Shanks' by Ezra Brown, but at the cost of using Math::BigInt.
```#!/usr/bin/perl -w

use Math::BigInt ':constant';
use strict;

sub tonelli {
my (\$s, \$e, \$n, \$x, \$b, \$g, \$r, \$m, \$t);
my (\$a, \$p) = @_;

die "\$a has no square roots (mod \$p)" if \$a**((\$p-1)/2) % \$p == -1;

\$s = \$p-1;
\$e = 0;
while (\$s % 2 == 0) {
\$s = \$s / 2;
\$e++;
}

for (\$n = 2; \$n >= 2; \$n++) {
last if \$n**((\$p-1)/2) % \$p == \$p-1;
}

\$x = \$a**((\$s+1)/2) % \$p;
\$b = \$a**\$s % \$p;
\$g = \$n**\$s % \$p;
\$r = \$e;

\$t = \$b;
while (1) {
for (\$m = 0; \$m <= \$r-1; \$m++) {
last if \$t % \$p == 1;
\$t = \$t**2;
}

return \$x if \$m == 0;

\$x = \$x * \$g**(2*(\$r-\$m-1)) % \$p;
\$b = \$b * \$g**(2*(\$r-\$m-1)) % \$p;
\$g = \$g**(2*(\$r-\$m)) % \$p;
\$r = \$m;
}
}

Create A New User
Node Status?
node history
Node Type: note [id://174242]
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others exploiting the Monastery: (5)
As of 2017-08-19 20:49 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
Who is your favorite scientist and why?

Results (312 votes). Check out past polls.

Notices?