# primep($x) returns whether $x is a prime
sub primep_rm {
 my($cand, $base, $p, $iter, $exp, $odd, $t, $x, $xprev);
 $cand = 0 + $_[0];
 1 < $cand && $cand == int($cand) or
  die "error: invalid number to primep_rm: $cand";
 for $p (@smallprimes) {
  0 != $cand % $p or
   return $p == $cand;
 }
 ($exp, $odd) = (0, $cand - 1);
 while (0 == $odd % 2) {
  ($exp, $odd) = ($exp + 1, $odd >> 1);
 }
 for $iter (1 .. 52) {
  $base = 1 + int(rand($cand - 1));
  $x = powmod($base, $odd, $cand);
  SQUARING: for $t (1 .. $exp) {
   ($xprev, $x) = ($x, mulmod($x, $x, $cand));
   $x == 1 and do {
    $xprev == 1 || $xprev == $cand - 1 or
     return;
    last;
   }
  }
  # KEY STEP: $x now contains $base**($cand - 1) % $cand
  # which should be 1 because of the Fermat theorem.
  $x == 1 or
   return;
 }
 return 1;
}