################################################
#
# Paul Garrett copyright 1995-2001
#
# Most recent additions "push down" on top
#
#####################################################################
######################################################################
######################################################################
package Num;
######################################################################
sub p_minus_one_verbose { ## Pollard's p-1: Returns process string
my $n = shift;
my $a = 3;
my $a_custom = shift;
if (defined $a_custom) {
$a = $a_custom;
}
my $rval = "";
$rval .= " Initialize \$ a=$a \$, and take list of primes (2,3).";
$rval .= " We'll look for primes \$ p \$ with \$ p-1 \$ most likely to be
a product of 2,3's. ";
my @p = (2,3); ##__EDIT__
my($e,$gcd);
for $p (@p) {
$e = int(log($n)/log($p));
$rval .= " Compute $e=floor(log($n)/log($p)). ";
$rval .= " Replace \$a=$a\$ by \$ a=$a^\{$p^\{$e\}\} \\red $n\$. ";
$a = &fastpow($a, $p ** $e, $n);
$rval .= " This gives \$ a = $a \$. ";
$gcd = &gcd($n, $a-1);
$rval .= " Compute \$\\gcd($n,$a-1)=$gcd\$.";
if ($gcd>1 && $gcd < $n) {
$rval .= " Thus, we find the proper factor $gcd. ";
return $rval;
}
}
$rval .= " No proper factor \$ p \$ with \$ p-1 \$ smooth was found. ";
return $rval;
}
######################################################################
sub p_minus_one_semiverbose { ## Pollard's p-1: Returns process string
my $n = shift;
my $a = 3;
my $a_custom = shift;
if (defined $a_custom) {
$a = $a_custom;
}
my $rval = "";
$rval .= " Initialize a=$a, and take list of primes (2,3).";
$rval .= " We'll look for primes p with p-1 most likely to be
a product of 2,3's. ";
my @p = (2,3); ##__EDIT__
my($e,$gcd);
for $p (@p) {
$e = int(log($n)/log($p));
$rval .= " Compute $e=floor(log($n)/log($p)). ";
$rval .= "\n";
$rval .= " Replace a=$a by a=$a^($p^$e) % $n. ";
$rval .= "\n";
$a = &fastpow($a, $p ** $e, $n);
$rval .= " This gives a = $a. ";
$rval .= "\n";
$gcd = &gcd($n, $a-1);
$rval .= " Compute gcd($n,$a-1)=$gcd.";
if ($gcd>1 && $gcd < $n) {
$rval .= " Thus, we find the proper factor $gcd. ";
$rval .= "\n";
return $rval;
}
}
$rval .= " No proper factor p with p-1 smooth was found. ";
return $rval;
}
######################################################################
sub p_minus_one { ## Pollard's p-1: Returns a factor, or 0
my $n = shift;
my $a = 3;
my $a_custom = shift;
if (defined $a_custom) {
$a = $a_custom;
}
my @p = (2,3); ##__EDIT__
my($e,$gcd);
for $p (@p) {
$e = int(log($n)/log($p));
$a = &fastpow($a, $p ** $e, $n);
$gcd = &gcd($n, $a-1);
if ($gcd>1 && $gcd < $n) {
return $gcd;
}
}
return 0;
}
######################################################################
sub rho_verbose { ## Pollard's rho. returns TeX string
my $n = shift;
if ($n >= 46341) { die "Need input < 46341 \n"; }
my $y = 2; ########## seed = 2
my $x = ($y*$y+1) % $n;
my $count = 0;
my $gcd;
my $rval = "";
$rval .= ' Initialize $y=2$ and $x=y^2+1 \red n$: ';
while ($count < $n) { ## stupid safety bound
$rval .= "\$ (x,y) = ($x,$y) \$. Compute ";
$gcd = &gcd($x-$y, $n);
$rval .= " \$\\gcd($x-$y,$n)=$gcd.\$";
if ($gcd > 1 && $gcd < $n) {
$rval .= " And $gcd is a proper factor of $n, so
we are done. ";
return $rval;
}
$rval .= " Replace \$y\$ by \$y^2+1 \\red $n \$, and ";
$rval .= " \$x\$ by \$(x^2+1)^2+1 \\red $n \$. ";
$y = ($y*$y+1) % $n;
$x = ($x*$x+1) % $n;
$x = ($x*$x+1) % $n;
$count++;
}
$rval .= " Didn't find a proper factor in $count tries. ";
return $rval;
}
######################################################################
sub rho_semiverbose { ## Pollard's rho. returns plaintext string
my $n = shift;
if ($n >= 46341) { die "Need input < 46341 \n"; }
my $y = 2; ########## seed = 2
my $x = ($y*$y+1) % $n; ########## x's doublespeed
my $count = 0;
my $gcd;
my $rval = "";
$rval .= " Initialize y=2 and x=y^2+1 % n: \n";
my($x_next,$y_next);
while ($count < $n) { ## stupid safety bound
$rval .= "(x,y) = ($x,$y) . Compute ";
$gcd = &gcd($x-$y, $n);
$rval .= " gcd($x-$y,$n)=$gcd. \n";
if ($gcd > 1 && $gcd < $n) {
$rval .= "And $gcd is a proper factor of $n, so
we are done. \n";
return $rval;
}
$y_next = ($y*$y+1) % $n;
$x_next = ($x*$x+1) % $n;
$x_next = ($x_next*$x_next+1) % $n;
$rval .= "Replace y=$y by $y^2+1 % $n = $y_next, and ";
$rval .= " x=$x by (x^2+1)^2+1 % $n = $x_next. \n";
$y = $y_next;
$x = $x_next;
$count++;
}
$rval .= " Didn't find a proper factor in $count tries. ";
return $rval;
}
######################################################################
sub rho { ## Pollard's rho. returns a factor
my $n = shift;
if ($n >= 46341) { die "Need input < 46341 \n"; }
my $y = 2; ########## seed = 2
my $x = ($y*$y+1) % $n;
my $count = 0;
my $gcd;
while ($count < $n) { ## stupid safety bound
$gcd = &gcd($x-$y, $n);
if ($gcd > 1 && $gcd < $n) {
return $gcd;
}
$y = ($y*$y+1) % $n;
$x = ($x*$x+1) % $n;
$x = ($x*$x+1) % $n;
$count++;
## print "(x,y) = ($x,$y)\n";
}
return 0;
}
######################################################################
sub is_primitive_root_verbose { ## Usage &is_primitive_root($g,$p)
my $g = shift;
my $p = shift;
my $expd;
my $query = 1;
my $out = "";
my @pfactors = &prime_factors($p-1);
while ( $q = pop( @pfactors )) {
$expd = &fastpow($g,($p-1)/$q,$p);
if ( $expd == 1 ) {
$out .= "\$ $g\^\{($p-1)/$q\} \\red $p = 1 \$";
$out .= " so $g fails the test to be a primitive root mod $p. ";
}
else {
$out .= "\$ $g\^\{($p-1)/$q\} \\red $p = $expd\\not= 1 \$";
if ( @pfactors ) {
$out .= " so we continue. ";
}
else {
$out .= " so $g passes the test to be a primitive root mod $p. ";
}
}
}
return $out;
"$g is a primitive root mod $p";
}
sub is_primitive_root { ## Usage &is_primitive_root($g,$p)
my $g = shift;
my $p = shift;
my $d;
for $q (&prime_factors($p-1)) {
if ( &fastpow($g,($p-1)/$q,$p) == 1 ) {
return 0;
}
}
return 1;
}
sub prime_factors {
my $n = shift;
my @factors = ();
if ($n % 2 == 0) {
push(@factors, 2);
}
my $d = 3;
while ($d <= sqrt($n)) {
if ($n % $d == 0) {
push(@factors, $d) if &is_prime( $d );
}
$d += 2;
}
return @factors;
}
sub is_prime { # usage &is_prime( $n )
my $n = shift;
my $d=1;
if ($n == 2 || $n == 3) {
return 1;
}
if ($n % 2 == 0) {
return 0;
}
my $sqrt = sqrt( $n );
while ($d <= $sqrt) {
$d += 2;
return 0 if $n % $d == 0;
}
return 1;
}
######################################################################
sub square_root_verbose { # probabilistic test, finds non-square first
my($rt);
my($b) = shift;
my($p) = shift;
my($g) = nonsquare( $p );
my($s,$m) = odd_part( $p-1);
print " We first factor $p-1 as \$ $p-1 = 2^s\\cdot m = 2^\{$s\}\\cdot $m\$. ";
my($B,$G,$E) = ( fastpow($b,$m,$p), fastpow($g,$m,$p), 0);
my($i);
$s_minus_1 = $s-1;
print " We look at indices i in the range 2..s,
which here is 2..$s";
foreach $i (2..$s) {
print "Compute \$\$ (B \\cdot G^E)^\{2^\{s-i\}\} ";
print " =($B \\cdot $G^\{$E\})^\{2^\{$s-$i\}\} = ";
print fastpow( ($B* fastpow($G,$E,$p)), 2**($s-$i), $p ), "\\mod $p\$\$";
if ( fastpow( ($B* fastpow($G,$E,$p)), 2**($s-$i), $p ) != 1 ) {
$E += 2**($i-1);
}
print "So E is now ", $E, ".\n ";
}
my($F) = $p-1-$E;
print " The final E is $E, and \$\$ F=p-1-E= $p-1-$E = $F \$\$";
print " Plug into the formula: the square root is \$\$ ";
print " (b \\cdot g^E)^\{(m+1)/2} \\cdot g^\{F/2\} = ";
print " ($b \\cdot $g^\{$E\})^\{($m+1)/2} \\cdot $g^\{$F/2\} = ";
$rt = (fastpow( $b * fastpow($g,$E,$p), ($m+1)/2, $p)
* fastpow($g,$F/2,$p)) % $p;
print " $rt \$\$ Thus, $rt is the square root. Check by computing";
print " \$ $rt\\cdot $rt= $b\\mod $p \$";
}
######################################################################
sub square_root { # probabilistic test, finds non-square first
my($rt);
my($b) = shift;
my($p) = shift;
my($g) = nonsquare( $p );
my($s,$m) = odd_part( $p-1);
my($B,$G,$E) = ( fastpow($b,$m,$p), fastpow($g,$m,$p), 0);
my($i);
foreach $i (2..$s) {
if ( fastpow( ($B* fastpow($G,$E,$p)), 2**($s-$i), $p ) != 1 ) {
$E += 2**($i-1);
}
}
my($F) = $p-1-$E;
$rt = (fastpow( $b * fastpow($g,$E,$p), ($m+1)/2, $p)
* fastpow($g,$F/2,$p)) % $p;
return $rt;
}
######################################################################
sub odd_part { # finds odd factor and power of 2
my($n) = shift;
my($s) = 0;
while ($n % 2 == 0 ){
$n /= 2;
$s++;
}
return ($s,$n);
}
######################################################################
sub nonsquare { # finds non-square mod prime p
my($p) = shift;
my($g) = 2;
while ( fastpow($g, ($p-1)/2, $p) == 1 ) {
$g++;
}
return $g;
}
######################################################################
sub distinct_factors { # produces list of _distinct_ factors
my($x) = shift;
my(@factors) = ();
my($d) = 2;
until ($x == 1) {
if ($x % $d == 0) {
push(@factors, $d);
$x /= $d;
while ( $x % $d == 0 ) {
$x /= $d;
}
}
else {
$d += 1;
}
}
return @factors;
}
######################################################################
sub tex_half_euclid {
my($x,$y) = @_;
my $eu = &print_half_euclid( $x, $y );
my(@eu) = split(/\n+/,$eu);
my $line;
for $line (@eu) {
$line =~ s/\./\\cdot /g;
$line .= '\hfill\cr';
$line =~ s/=/\&=\&/; # just once on each line
}
$eu = join("\n", @eu);
$eu = '$$\matrix{'."\n".$eu."\n".'}$$';
return $eu; ## modified to _return_ rather than _print_!!!
}
######################################################################
sub tex_full_euclid {
my($x,$y) = @_;
my $eu = &print_full_euclid( $x, $y );
my(@eu) = split(/\n+/,$eu);
my $line;
for $line (@eu) {
$line =~ s/\./\\cdot /g;
$line .= '\hfill\cr';
$line =~ s/=/\&=\&/; # just once on each line
}
$eu = join("\n", @eu);
$eu = '$$\matrix{'."\n".$eu."\n".'}$$';
return $eu;
}
######################################################################
sub prime {
my $n = shift;
$n = eval $n;
($n < 2147483647) || die "$num too big for this implementation\n";
if ($n == 2) {
return 1;
}
if ($n % 2 == 0) {
return 0;
}
my $div = 3;
my $rt = int(sqrt($n));
while ($div <= $rt) {
if ($n % $div == 0) {
return 0;
}
else {
$div += 2;
}
}
return 1; # for 'true': it is prime
}
######################################################################
sub miller_rabin_verbose {
## tests $n for psprimality base $b
my($n, $b) = @_; # changed order or args 11-27-98
# changed again, 01-09-99
my($s,$m) = (0,$n-1);
while ($m % 2 == 0) {
$s++;
$m /= 2;
}
my($c) = &fastpow($b,$m,$n);
my $rval = "
Apply the Miller-Rabin test with a single auxiliary number, the base \$ b = $b\$.
Factor \$ $n-1 = 2^s\\cdot m \$ with \$ m \$ odd: here \$ m = $m \$ and
\$ s = $s \$. Then (using the fast exponentiation algorithm) compute
\$ c = b^m \\mod $n = $b^\{$m\} = $c \\mod $n\$.";
if ($c == 1 || $c == $n-1) {
$rval .= " Since \$ $c = \\pm 1 \\mod $n \$ we conclude that ";
$rval .= " \$ $n \$ is a strong pseudoprime base $b.\n";
return $rval;
}
else {
$rval .= " Since \$ $c \\not= \\pm 1 \\mod $n \$ we continue, ";
$rval .= " entering the squaring loop. ";
}
my($t) = 0;
while ($t < $s-1) {
$t++;
$rval .= " Square $c mod $n (and increment \$t\$ to $t), obtaining ";
$c = ($c * $c) % $n;
$rval .= "$c. ";
if ($c == $n-1) {
$rval .= "Since this is -1 mod $n, (and still \$ $t=t < s=$s \$) we conclude that ";
# print "we won't find an anomalous square root of 1 here, so ";
$rval .= " \$ $n \$ is a strong pseudoprime base $b.\n";
return $rval;
}
elsif ($c == 1) {
# $rval .= " Since we have obtained +1 prior to obtaining -1,";
# $rval .= " we have found a square root of 1 other than \$ \\pm 1\$,";
$rval .= " so we conclude that ";
$rval .= " \$ $n \$ is definitely composite ";
$rval .= "(since it is not a strong pseudoprime base $b.\n";
return $rval;
}
}
if ($c != 1) {
# $rval .= "That is, $n is not even a Fermat pseudoprime, so ";
$rval .= " It is not a strong pseudoprime.\n";
return $rval;
}
# $rval .= "Since we have obtained +1 prior to obtaining -1,";
# $rval .= " we have found a square root of 1 other than \$ \\pm 1\$,";
$rval .= " so we conclude that ";
$rval .= " \$ $n \$ is not a strong pseudoprime base $b.\n";
return $rval;
}
######################################################################
######################################################################
sub char_to_int { ## lowercases and puts into range 0-25
my($char) = @_;
return unpack("C",lc($char)) - 97;
}
sub int_to_char { ## range 0-25 to a-z (lowercase)
my($int) = @_;
return pack("C",$int+97);
}
sub word_to_ints {
my($in) = @_; ## a string
my(@in) = split(//,$in);
my($i);
foreach $i (0..$#in) {
$in[$i] = char_to_int($in[$i]);
}
return @in; ################## edited 07-13-98
}
sub ints_to_word {
my(@in) = @_;
my($i);
foreach $i (0..$#in) {
$in[$i] = int_to_char($in[$i]);
}
return join('',@in);
}
####################################################################
sub miller_rabin { # does single Miller-Rabin test for strong pseudoprimes
# tests n via auxiliary b
my($n) = shift;
my($b) = shift; #### changed order of args 11-27-98
my($m) = $n-1; #### changed again 11-25-2000
my($power_of_two) = 1;
my($s) = 0; ############### traditional but unsuggestive notation...
while ($m % 2 == 0) {
$m /= 2;
$power_of_two *= 2;
$s++;
} ############################ the above need not be re-computed for
############################ other $b's
my($c) = &fastpow($b,$m,$n);
if ($c == 1 || $c == $n-1) {
return 'prime'; # getting +1/-1 immediately
# surely won't find a bogus sqrt(1) this way
}
my($t) = 0;
while ($t < $s-1) {
$t++;
$c = ($c * $c) % $n;
if ($c == $n-1) { # get -1 mod $n
return 'prime';
}
elsif ($c == 1) { # get +1
return 'composite1';
}
}
# by now, $c is $b**($n-1), but has failed Fermat pseudoprime test
# or else satisfies Fermat, but previous number was _not_ -1, so we
# have a "bad" sqrt(1) in that case...
return 'composite2';
}
######################################################################
sub miller_rabin_semiverbose { ## in plain text, not TeX...
## tests $n for psprimality base $b
my($n, $b) = @_;
print "Is $n a strong pseudoprime base $b? \n";
my($s,$m) = (0,$n-1);
while ($m % 2 == 0) {
$s++;
$m /= 2;
}
print "$n-1 = 2^s * m with s=$s and m=$m \n";
my $c = Num::fastpow($b,$m,$n);
print "b^m % $n = $b^$m = $c mod $n. \n";
if ($c == 1 || $c == $n-1) {
print "Since $c=+/- 1 mod $n, $n is a strong pseudoprime base $b. \n";
return;
}
else {
print "Since $c is not +1/-1 mod $n, continue: \n";
print "Initialize t=0. \n";
}
my $t = 0;
while ($t < $s-1) {
my $t_old = $t;
$t++;
print "Replace $c by its square ";
$c = ($c * $c) % $n;
print "$c mod $n (and increment t=$t_old to $t) \n";
if ($c == $n-1) {
print "Since this is -1 mod $n, $n is a strong pseudoprime base $b.\n";
return;
}
elsif ($c == 1) {
print "So $n is definitely composite. \n";
my @factors = Num::factors($n);
my $smallest = @factors -> [0];
print "(Indeed, $n has smallest factor $smallest.) \n";
return;
}
}
## by now $t == $s !?
if ($c != 1) {
print "t = s-1 with c value not -1 so $n is definitely composite. \n";
my @factors = Num::factors($n);
my $smallest = @factors -> [0];
print "(Indeed, $n has smallest factor $smallest.) \n";
}
}
######################################################################
sub factors { ## produces array of factors, in ascending order
my($x) = @_;
my(@factors) = ();
my($d) = 2;
until ($d > sqrt($x)) {
if ($x % $d == 0) {
push(@factors, $d);
$x = $x/$d;
}
else {
$d += 1;
}
}
push(@factors,$x);
return @factors;
}
sub fastpow { ## fastpow(x,e,m) is x^e % m
my($x,$e,$m) = @_;
$out = 1;
while ($e > 0) {
if ($e % 2 == 1) {
$out = ($out * $x) % $m;
$e -= 1;
} ## boy am I clever.......
$x = ($x * $x) % $m;
$e /= 2;
}
return $out;
}
sub fastpow_verbose { ## fastpow(x,e,m) is x^e % m
my($x,$e,$m) = @_;
my $out = 1;
my $rvalue = "";
$rvalue .= "\$ ($x, $e, $out) \$, ";
while ($e > 0) {
if ($e % 2 == 1) {
$out = ($out * $x) % $m;
$e -= 1;
$rvalue .= "\$ ($x, $e, $out) \$";
}
else {
$x = ($x * $x) % $m;
$e /= 2;
$rvalue .= "\$ ($x, $e, $out) \$";
}
$rvalue .= ", " unless $e == 0;
}
return $rvalue;
}
sub fermat_psp { #__EDIT__: could exceed implementation bound...
my($n, $b) = @_;
my $pow = &fastpow($b,$n,$n);
if ($pow == $b) {
return "$n is Fermat pseudoprime base $b";
}
else {
return "$n is NOT Fermat pseudoprime base $b";
}
}
sub carmichael { # test primes... less than sqrt...
my($n) = shift;
my($b) = 2;
my($sqrt) = int(sqrt($n));
if (&fastpow($b,$n,$n) != $b) {
return "$n is not Fermat pseudoprime base $b\n";
}
my $nn = &prime($n); #__EDIT__eh?__
if ($nn == $n) { # only test this if base 2 Fermat pseudoprime already
return "$n is an actual prime\n";
}
$b++;
while ($b <= $sqrt) {
if (&fastpow($b,$n,$n) != $b) {
return "$n is not Fermat pseudoprime base $b\n";
}
$b += 2;
}
return "$n is a Carmichael number\n";
}
sub list_carmichael {
my $lower = shift;
my $upper = shift;
if ($lower % 2 == 0) {
$lower++;
}
my $n;
my $ret;
for ($n = $lower; $n < $upper; $n += 2) {
$ret = &carmichael($n);
print $ret if ($ret =~ /Carmichael/);
}
}
######################################################################
sub order { ## of x mod m ## Nov 21, 97
my($x,$m) = @_;
my($e) = 1;
my($y) = $x;
if (gcd($x,$m) == 1) {
until ($y == 1) {
$y = ($y * $x) % $m;
$e++;
}
}
return $e;
}
sub factorial { ## only for small...?
my($n) = @_;
my($out) = 1;
for ($i=1; $i<=$n; $i++) {
$out *= $i;
}
return $out;
}
sub binomial_coefficient {
my($n,$k) = @_;
my($out) = 1;
for ($i=$n; $i>$n-$k; $i--) {
$out *= $i;
$out /= $n-$i+1;
}
return $out;
}
sub slow_inverse { ##__EDIT__: very slow...
my($in, $modulus) = @_; ## presume $in is odd, not 13, etc...
my($out) = 0;
if (&gcd($in,$modulus) == 1) {
until (($out * $in) % $modulus == 1) {
$out++;
}
}
else {
$out = 0; ## error state
}
return $out;
}
sub full_euclid { ## use Euclidean Algorithm # revised 18 Oct 98: signs...
my($y) = shift;
my($x) = shift;
my($r,$q);
my($y_save) = $y;
my($x_save) = $x;
$x = abs($x); # revised 18 Oct 98
$y = abs($y);
my($out) = 0; ## with actual return of 0 as error msg
my($a,$b,$c,$d,$atmp,$btmp);
$a = 1; $d = 1; $b = 0; $c = 0;
while ($y != 0) {
my($q) = int($x/$y);
$atmp = $a; $btmp = $b;
$a = $c;
$b = $d;
$c = $atmp - $q * $c;
$d = $btmp - $q * $d;
$r = $x % $y;
$x = $y;
$y = $r;
}
$a = -$a if $x_save < 0;
$b = -$b if $y_save < 0;
#
# first input is $y, second is $x
#
return ($b,$y_save,$a,$x_save,$x); # form "by + ax = g"
}
sub gcd { ## repeats some of previous, but ok...
my($x,$y) = @_;
my($remainder) = $y;
until ($remainder == 0) {
$remainder = $x % $y;
$x = $y;
$y = $remainder;
}
return abs($x);
}
sub inverse {
my($x,$m) = @_;
my @eu = &full_euclid($x,$m);
if ((@eu -> [4]) == 1) {
(@eu -> [0] + $m) % $m;
}
else {
return 0; # error...
}
}
######################################################################
sub print_full_euclid {
my($x, $y) = @_;
$out = "";
@x = ($x);
@y = ($y);
(@r, @q) = ();
my $q;
my $r = $y;
until ($y == 0) {
$q = int($x/$y);
push(@q, $q);
$r = $x % $y;
push(@r, $r);
$out .= "$x - $q.$y = $r\n";
($x,$y) = ($y,$r);
push(@x,$x) unless ($r ==0);
push(@y,$y) unless ($r ==0);
}
pop @x;
pop @y;
pop @q;
pop @r;
$out .= "\n";
($x,$y,$q,$r) = (pop @x, pop @y, pop @q, pop @r);
my $c1 = 1;
my $c2 = -$q;
$out .= "$r = ($c1)$x + ($c2)$y ";
while (@q) {
($x,$y,$q,$r) = (pop @x, pop @y, pop @q, pop @r);
$out .= " = ($c1)$y + ($c2)($x - $q.$y) \n ";
($c1, $c2) = ($c2, $c1 - $c2*$q);
$out .= " = ($c1)$x + ($c2)$y ";
}
$out .= "\n";
return $out;
}
######################################################################
sub print_half_euclid {
my($x, $y) = @_;
$out = "";
my $q;
my $r = $y;
until ($y == 0) {
$q = int($x/$y);
$r = $x % $y;
push(@r, $r);
$out .= "$x - $q.$y = $r\n";
($x,$y) = ($y,$r);
}
return $out;
}
##################################################
sub prime_or_factor { # can pass evaluable expressions...
my $n = shift;
$n = eval $n;
($n < 2147483647) || die "$num too big for this implementation\n";
if ($n % 2 == 0) {
return 2;
}
my $div = 3;
my $rt = int(sqrt($n));
my $known_factor = 'no';
while ($div <= $rt) {
if ($n % $div == 0) {
return $div;
}
else {
$div += 2;
}
}
return $n;
}
#################################################3
## Nov 7, 1997
sub print_cycle {
my($hash_name) = @_; ## ref to name of hash
my($out) = "" ;
delete %{hash_name} -> {0} if (defined %{hash_name} -> {0});
my(@keys) = keys %{$hash_name}; ## these should be among 1,2,3,...,N
if ($#keys >= 0) {
my($first_key) = $keys[0];
$out .= "$first_key ";
my($next_key) = %{$hash_name} -> {$first_key};
until ($next_key eq $first_key) {
$out .= "$next_key ";
$next_key = %{$hash_name} -> {$next_key};
}
}
return $out;
}
##########################################################
sub decompose_into_cycles { ## this is deprecated form
my(@perm) = @_; ## array 0, img(1), img(2), ...
my(%perm) = ();
foreach $i (1..7) {
%perm -> {"$i"} = @perm -> [$i];
}
my(@keys) = sort keys %perm;
my($index) = 0;
my(@cycles);
my($key);
until ($#keys < 0) {
$key = @keys -> [0]; ## effectively at random...
$init_key = $key;
@cycles -> [$index] -> {$key} = %perm -> {$key};
$next_key = %perm -> {$key};
until ($next_key eq $init_key) {
@cycles -> [$index] -> {$next_key} = %perm -> {$next_key};
$key = $next_key;
$next_key = %perm -> {$key};
}
foreach $key (keys %{@cycles -> [$index]}) {
delete %perm -> {$key};
}
## then repeat:
@keys = keys %perm;
$index++;
}
my($out) = ""; ## string expressing the cycles...
foreach $i (0..$#cycles) {
$out .= '(';
$out .= &print_cycle(\%{@cycles->[$i]});
$out .= ')';
}
return $out;
}
######################################################################
sub cross_product {
my($x,$y,$z,$X,$Y,$Z) = @_;
my @out;
@out -> [0] = $y * $Z - $Y * $z ;
@out -> [1] = - $x * $Z + $X * $z ;
@out -> [2] = $x * $Y - $X * $y ;
return @out;
}
#####################################################################
sub qr_main { # returns (-1)^(x-1)(y-1)/4
my ($x,$y) = @_;
if (($x % 4 == 1) || ($y % 4 == 1)) {
return 1;
}
else {
return -1;
}
}
sub qr_two { # returns (2/x)_2
my $x = shift;
if (($x % 8 == 1) || ($x % 8 == 7)) {
return 1;
}
else {
return -1;
}
}
sub qr_minus_one { # returns (-1/x)_2
my $x = shift;
if ($x % 4 == 1) {
return 1;
}
else {
return -1;
}
}
sub twos { # returns ord_2(x)
my $x = shift;
my $log_twos = 0;
while ($x % 2 == 0) {
$x /= 2;
$log_twos++;
}
return $log_twos;
}
sub minus_one_pow { # returns (-1)^x
my $x = shift;
if ($x % 2 == 0) {
return 1;
}
else {
return -1;
}
}
######################################## QR stuff above...
sub quadratic_symbol_verbose {
my($y,$m) = @_; # computes (y/m)_2
my $sgn = 1;
my $tmp_sgn;
my $twos;
print "\$\$\\eqalign\{ ($y/$m)_2 & = \\cr\n"; # top
while ($y > 1) {
print " & = "; # beginning of line
print " ($sgn)\\cdot " if ($sgn == -1);
$twos = &twos($y,$m); # number of factors of 2
$y /= 2**$twos;
if ($twos > 0) { # there are some 2's spit out
if ($twos == 1) {
print " (2/$m)_2 \\cdot ($y/$m)_2 ";
}
else {
print " (2/$m)_2^$twos \\cdot ($y/$m)_2 ";
}
if ($twos % 2 == 0) {
# do nothing
}
else {
$sgn *= &qr_two($m);
}
print " = ($sgn) \\cdot ($y/$m)_2";
print " \\hbox\{\\hskip30pt (by two-part of Quad.Rec.) \}";
} # end of handling of twos, all on the same line
else { # no twos spit out
($y,$m) = ($m,$y);
$tmp_sgn = qr_main($y,$m);
print " ($tmp_sgn)\\cdot($y/$m)_2 ";
$y %= $m;
$sgn *= $tmp_sgn;
print " = ";
print " ($sgn) \\, " if ($sgn == -1);
print " ($y/$m)_2 \\hbox\{\\hskip30pt (by main part of Quad. Rec.) \}";
}
print " \\cr\n ";
} # end of $y > 1 loop
# by now $y == 1 :
print " & = $sgn \}\$\$ "; # end of equalign
}
######################################################################
sub two_digit_primes {
return (11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97);
}
sub three_digit_primes {
return (101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293,
307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,
383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457,
461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547,
557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619,
631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709,
719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883,
887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983,
991, 997);
}
sub small_primes {
return (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191,
193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263,
269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347,
349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421,
431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499,
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593,
599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661,
673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757,
761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853,
857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
947, 953, 967, 971, 977, 983, 991, 997);
}
######################################################################
1;
###############################################
##
## The End
##
################################################