# Prime.pm

# $Id: Prime.pm,v 1.2 1999/12/30 19:35:33 friedman Exp $

# Commentary:
# Code:

package NF::Prime;

BEGIN
{
  use Exporter ();
  use vars qw ($VERSION @ISA @EXPORT_OK);

  # set the version for version checking
  $VERSION     = 1.00;

  @ISA         = qw (Exporter);
  @EXPORT_OK   = qw (primep next_prime);
}


# Compute the first prime following (or including) n.
#
# TODO: Consider whether it's worth caching computed primes in a table.
# For now this is fast enough.
sub next_prime
{
  my ($n) = @_;

  $n++ while (! primep ($n));
  return $n;
}

# Miller-Rabin probablistic test for primality.
sub primep
{
  use integer;
  my $n  = shift;
  my $p  = shift;

  my $n1 = $n - 1;
  my $one = $n - $n1; # 1, but ensure the right type of number.

  # find the largest power of two less than n-1.
  my $p2;
  my $p2index;
  if (defined ($p))
    {
      $p2index = $p;
      $p2 = (2 ** $p);
    }
  else
    {
      $p2 = $one;
      $p2index++, $p2 *= 2
        while $p2 < $n1;
    }

  # number of iterations: 5 for 260-bit number, go up to
  # 25 for much smaller numbers.
  my $last_witness = 5;
  $last_witness += (260 - $p2index)/13 if $p2index < 260;

  my $witness = $one * 100;
  for $witness_count ( 1..$last_witness )
    {
      $witness *= 1024;
      $witness += int(rand(1024));
      $witness = $witness % $n if $witness > $n;
      $witness = $one * 100, redo if $witness == 0;

      my $prod = $one;
      my $n1bits = $n1;
      my $p2next = $p2;

      # compute $witness ** ($n - 1).
      while (1)
        {
          # Is $prod, the power so far, a square root of 1?
          # (plus or minus 1)
          my $rootone = ($prod == 1 || $prod == $n1);

          $prod = ($prod * $prod) % $n;

          # An extra root of 1 disproves the primality.
          return 0 if ($prod == 1 && !$rootone);

          if ( $n1bits >= $p2next )
            {
              $prod = ($prod * $witness) % $n;
              $n1bits -= $p2next;
            }

          last if $p2next == 1;
          $p2next /= 2;
        }
      return 0 unless $prod == 1;
    }
  return 1;
}

1;

# eof