Beefy Boxes and Bandwidth Generously Provided by pair Networks
The stupid question is the question not asked
 
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
in thread Vampire Numbers

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; } }


Comment on Re: Re: Re: Factors
Download Code
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; } }

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://174242]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others perusing the Monastery: (12)
As of 2015-07-31 18:42 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The top three priorities of my open tasks are (in descending order of likelihood to be worked on) ...









    Results (280 votes), past polls