Math::Primality - Advanced Primality Algorithms using GMP

    use Math::Primality qw/:all/;

    my $t1 = is_pseudoprime($x,$base);
    my $t2 = is_strong_pseudoprime($x);

    print "Prime!" if is_prime($outrageously_large_prime);

    my $t3 = next_prime($x); 

    if ( is_pseudoprime($n,$b) ) {
        # it's a pseudoprime
    } else {
        # not a psuedoprime
    }

    if ( is_strong_pseudoprime($n,$b) ) {
        # it's a strong pseudoprime
    } else {
        # not a strong psuedoprime
    }

    if ( is_strong_lucas_pseudoprime($n) ) {
        # it's a strong Lucas-Selfridge pseudoprime
    } else {
        # not a strong Lucas-Selfridge psuedoprime
        # i.e. definitely composite
    }

    if ( is_prime($n) ) {
        # it's a prime
    } else {
        # definitely composite
    }

    my $q = next_prime($n);

    my $q = prev_prime($n);

    my $count = prime_count(1000);          # $count = 168
    my $bigger_count = prime_count(10000);  # $bigger_count = 1229

    perldoc Math::Primality

package Math::Primality;
use warnings;
use strict;
use Data::Dumper;
use Math::GMPz qw/:mpz/;
use base 'Exporter';
use Carp qw/croak/;
my %small_primes = (
    2   => 2,
    3   => 2,
    5   => 2,
    7   => 2,
    11  => 2,
    13  => 2,
    17  => 2,
    19  => 2,
    23  => 2,
    29  => 2,
    31  => 2,
    37  => 2,
    41  => 2,
    43  => 2,
    47  => 2,
    53  => 2,
    59  => 2,
    61  => 2,
    67  => 2,
    71  => 2,
    73  => 2,
    79  => 2,
    83  => 2,
    89  => 2,
    97  => 2,
    101 => 2,
    103 => 2,
    107 => 2,
    109 => 2,
    113 => 2,
    127 => 2,
    131 => 2,
    137 => 2,
    139 => 2,
    149 => 2,
    151 => 2,
    157 => 2,
    163 => 2,
    167 => 2,
    173 => 2,
    179 => 2,
    181 => 2,
    191 => 2,
    193 => 2,
    197 => 2,
    199 => 2,
    211 => 2,
    223 => 2,
    227 => 2,
    229 => 2,
    233 => 2,
    239 => 2,
    241 => 2,
    251 => 2,
    257 => 2,
);

use constant
	DEBUG => 0
;

use constant GMP => 'Math::GMPz';

our $VERSION = '0.04';
$VERSION = eval $VERSION;


our @EXPORT_OK = qw/is_pseudoprime is_strong_pseudoprime is_strong_lucas_pseudoprime is_prime next_prime prev_prime prime_count/;

our %EXPORT_TAGS = ( all => \@EXPORT_OK );

sub is_pseudoprime($;$)
{
    my ($n, $base) = @_;
    return 0 unless $n;
    $base ||= 2;
    # we should check if we are passed a GMPz object
    $base   = GMP->new($base);
    $n      = GMP->new($n);

    my $m    = GMP->new();
    Rmpz_sub_ui($m, $n, 1);              # $m = $n - 1

    my $mod = GMP->new();
    Rmpz_powm($mod, $base, $m, $n );     # $mod = ($base ^ $m) mod $n 
    return ! Rmpz_cmp_ui($mod, 1);       # pseudoprime if $mod = 1
}

# checks if $n is in %small_primes
# private functions expect a Math::GMPz object
sub _is_small_prime
{
    my $n = shift;
    $n = Rmpz_get_ui($n);
    return $small_primes{$n} ? 2 : 0;

}

sub debug {
    if ( DEBUG ) {
      warn $_[0];
    }
}

sub is_strong_pseudoprime($;$)
{
    my ($n, $base) = @_;

    $base ||= 2;
    # we should check if we are passed a GMPz object
    $base   = GMP->new($base);
    $n      = GMP->new($n);

    # unnecessary but faster if $n is even
    my $cmp = _check_two_and_even($n);
    return $cmp if $cmp != 2;

    my $m   = GMP->new();
    Rmpz_sub_ui($m,$n,1);              # $m = $n - 1

    my ($s,$d) = _find_s_d($m);
    debug "m=$m, s=$s, d=$d" if DEBUG;

    my $residue = GMP->new();
    Rmpz_powm($residue, $base,$d, $n); # $residue = ($base ^ $d) mod $n
    debug "$base^$d % $n = $residue" if DEBUG;

    # if $base^$d = +-1 (mod $n) , $n is a strong pseudoprime

    if ( Rmpz_cmp_ui($residue,1) == 0 ) {
        debug "found $n as spsp since $base^$d % $n == $residue == 1\n" if DEBUG;
        return 1;
    }

    if ( Rmpz_cmp($residue,$m) == 0 ) {
        debug "found $n as spsp since $base^$d % $n == $residue == $m\n" if DEBUG;
        return 1;
    }

    map {
        # successively square $residue, $n is a strong psuedoprime
        # if any of these are congruent to -1 (mod $n)
        Rmpz_mul($residue,$residue,$residue);   # $residue = $residue * $residue
        debug "$_: r=$residue" if DEBUG;

        my $mod = GMP->new();
        Rmpz_mod($mod, $residue, $n);           # $mod = $residue mod $n
        debug "$_:$residue % $n = $mod " if DEBUG;
        $mod = Rmpz_cmp($mod, $m);

        if ($mod == 0) {
            debug "$_:$mod == $m => spsp!" if DEBUG;
            return 1;
        }
    } ( 1 .. $s-1 );

    return 0;
}

# given an odd number N find (s, d) such that N = d * 2^s + 1
# private functions expect a Math::GMPz object
sub _find_s_d($)
{
    my $m   = $_[0];
    my $s   = Rmpz_scan1($m,1);
    my $d   = GMP->new();
    Rmpz_tdiv_q_2exp($d,$m,$s);
    return ($s,$d);
}

sub is_strong_lucas_pseudoprime($)
{
    my ($n) = @_;
    $n      = GMP->new($n);
    # we also need to handle all N < 3 and all even N 
    my $cmp = _check_two_and_even($n);
    return $cmp if $cmp != 2;
    # handle all perfect squares
    if ( Rmpz_perfect_square_p($n) ) {
        return 0;
    }
    # determine Selfridge parameters D, P and Q
    my ($D, $P, $Q) = _find_dpq_selfridge($n);
    if ($D == 0) {  #_find_dpq_selfridge found a factor of N
      return 0;
    }
    my $m = GMP->new();
    Rmpz_add_ui($m, $n, 1);  # $m = $n + 1

    # determine $s and $d such that $m = $d * 2^$s + 1
    my ($s,$d) = _find_s_d($m);
    # compute U_d and V_d
    # initalize $U, $V, $U_2m, $V_2m
    my $U = GMP->new(1);     # $U = U_1 = 1
    my $V = GMP->new($P);    # $V = V_1 = P
    my $U_2m = GMP->new(1);  # $U_2m = U_1
    my $V_2m = GMP->new($P); # $V_2m = P
    # initalize Q values (eventually need to calculate Q^d, which will be used in later stages of test)
    my $Q_m = GMP->new($Q);
    my $Q2_m = GMP->new(2 * $Q);  # Really 2Q_m, but perl will barf with a variable named like that
    my $Qkd = GMP->new($Q);
    # start doubling the indicies!
    my $dbits = Rmpz_sizeinbase($d,2);
    for (my $i = 1; $i < $dbits; $i++) {  #since d is odd, the zeroth bit is on so we skip it
      # U_2m = U_m * V_m (mod N)
      Rmpz_mul($U_2m, $U_2m, $V_2m);  # U_2m = U_m * V_m
      Rmpz_mod($U_2m, $U_2m, $n);     # U_2m = U_2m mod N
      # V_2m = V_m * V_m - 2 * Q^m (mod N)
      Rmpz_mul($V_2m, $V_2m, $V_2m);  # V_2m = V_2m * V_2m
      Rmpz_sub($V_2m, $V_2m, $Q2_m);  # V_2m = V_2m - 2Q_m
      Rmpz_mod($V_2m, $V_2m, $n);     # V_2m = V_2m mod N
      # calculate powers of Q for V_2m and Q^d (used later)
      # 2Q_m = 2 * Q_m * Q_m (mod N)
      Rmpz_mul($Q_m, $Q_m, $Q_m);     # Q_m = Q_m * Q_m
      Rmpz_mod($Q_m, $Q_m, $n);       # Q_m = Q_m mod N
      Rmpz_mul_2exp($Q2_m, $Q_m, 1);  # 2Q_m = Q_m * 2
      if (Rmpz_tstbit($d, $i)) {      # if bit i of d is set
        # add some indicies
        # initalize some temporary variables
        my $T1 = GMP->new();
        my $T2 = GMP->new();
        my $T3 = GMP->new();
        my $T4 = GMP->new();
        # this is how we do it
        # U_(m+n) = (U_m * V_n + U_n * V_m) / 2
        # V_(m+n) = (V_m * V_n + D * U_m * U_n) / 2
        Rmpz_mul($T1, $U_2m, $V);     # T1 = U_2m * V
        Rmpz_mul($T2, $U, $V_2m);     # T2 = U * V_2m
        Rmpz_mul($T3, $V_2m, $V);     # T3 = V_2m * V
        Rmpz_mul($T4, $U_2m, $U);     # T4 = U_2m * U
        Rmpz_mul_si($T4, $T4, $D);    # T4 = T4 * D = U_2m * U * D
        Rmpz_add($U, $T1, $T2);       # U = T1 + T2 = U_2m * V - U * V_2m
        if (Rmpz_odd_p($U)) {         # if U is odd
          Rmpz_add($U, $U, $n);       # U = U + n
        }
        Rmpz_fdiv_q_2exp($U, $U, 1);  # U = floor(U / 2)
        Rmpz_add($V, $T3, $T4);       # V = T3 + T4 = V_2m * V + U_2m * U * D
        if (Rmpz_odd_p($V)) {         # if V is odd
          Rmpz_add($V, $V, $n);       # V = V + n 
        }
        Rmpz_fdiv_q_2exp($V, $V, 1);  # V = floor(V / 2)
        Rmpz_mod($U, $U, $n);         # U = U mod N
        Rmpz_mod($V, $V, $n);         # V = V mod N
        # Get our Q^d calculating on (to be used later)
        Rmpz_mul($Qkd, $Qkd, $Q_m);   # Qkd = Qkd * Q_m
        Rmpz_mod($Qkd, $Qkd, $n);     # Qkd = Qkd mod N
      }
    }
    # if U_d or V_d = 0 mod N, then N is prime or a strong Lucas pseudoprime
    if(Rmpz_sgn($U) == 0 || Rmpz_sgn($V) == 0) {
      return 1;
    }
    # ok, if we're still here, we have to compute V_2d, V_4d, V_8d, ..., V_{2^(s-1)*d}
    # initalize 2Qkd
    my $Q2kd = GMP->new;
    Rmpz_mul_2exp($Q2kd, $Qkd, 1);    # 2Qkd = 2 * Qkd
    # V_2m = V_m * V_m - 2 * Q^m (mod N)
    for (my $r = 1; $r < $s; $r++) {
      Rmpz_mul($V, $V, $V);      # V = V * V;
      Rmpz_sub($V, $V, $Q2kd);   # V = V - 2Qkd
      Rmpz_mod($V, $V, $n);      # V = V mod N
      # if V = 0 mod N then N is a prime or a strong Lucas pseudoprime
      if(Rmpz_sgn($V) == 0) {
        return 1;
      }
      # calculate Q ^(d * 2^r) for next r (unless on final iteration)
      if ($r < ($s - 1)) {
        Rmpz_mul($Qkd, $Qkd, $Qkd);     # Qkd = Qkd * Qkd
        Rmpz_mod($Qkd, $Qkd, $n);       # Qkd = Qkd mod N
        Rmpz_mul_2exp($Q2kd, $Qkd, 1);  # 2Qkd = 2 * Qkd
      } 
    }
    # otherwise N is definitely composite 
    return 0;
}

# selfridge's method for finding the tuple (D,P,Q) for is_strong_lucas_pseudoprime
# private functions expect a Math::GMPz object
sub _find_dpq_selfridge($) {
  my $n = $_[0];
  my ($d,$sign,$wd) = (5,1,0);
  my $gcd = GMP->new;

  # determine D
  while (1) {
    $wd = $d * $sign;

    Rmpz_gcd_ui($gcd, $n, abs $wd);
    if ($gcd > 1 && Rmpz_cmp($n, $gcd) > 0) {
      debug "1 < $gcd < $n => $n is composite with factor $wd" if DEBUG;
      return 0;
    }
    my $j = Rmpz_jacobi(GMP->new($wd), $n);
    if ($j == -1) {
      debug "Rmpz_jacobi($wd, $n) == -1 => found D" if DEBUG;
      last; 
    }
    # didn't find D, increment and swap sign
    $d += 2;
    $sign = -$sign;
  }
  # P = 1
  my ($p,$q) = (1,0);
  {
      use integer;
      # Q = (1 - D) / 4
      $q = (1 - $wd) / 4;
  }
  debug "found P and Q: ($p, $q)" if DEBUG;
  return ($wd, $p, $q);
}

# method returns 0 if N < two or even, returns 1 if N == 2
# returns 2 if N > 2 and odd
# private functions expect a Math::GMPz object
sub _check_two_and_even($) {
  my $n = $_[0];

  my $cmp = Rmpz_cmp_ui($n, 2);
  return 1 if $cmp == 0;
  return 0 if $cmp < 0;
  return 0 if Rmpz_even_p($n);
  return 2;
}

sub is_prime($) {
    my $n = shift;
    $n = GMP->new($n);

    if (Rmpz_cmp_ui($n, 2) == -1) {
        return 0;
    }
    if (Rmpz_cmp_ui($n, 257) == -1) {
        return _is_small_prime($n);
    } elsif ( Rmpz_cmp_ui($n, 9_080_191) == -1 ) {
        return 0 unless is_strong_pseudoprime($n,31);
        return 0 unless is_strong_pseudoprime($n,73);
        return 2;
    } elsif ( Rmpz_cmp_ui($n, 4_759_123_141) == -1 ) {
        return 0 unless is_strong_pseudoprime($n,2);
        return 0 unless is_strong_pseudoprime($n,7);
        return 0 unless is_strong_pseudoprime($n,61);
        return 2;
    }
    # the strong_pseudoprime test is quicker, do it first
    return is_strong_pseudoprime($n,2) && is_strong_lucas_pseudoprime($n);
}

sub next_prime($) {
  my $n = GMP->new($_[0]);
  my $cmp = Rmpz_cmp_ui($n, 2 ); #check if $n < 2
  if ($cmp < 0) {
    return GMP->new(2);
  }
  if (Rmpz_odd_p($n)) {         # if N is odd
    Rmpz_add_ui($n, $n, 2);     # N = N + 2
  } else {
    Rmpz_add_ui($n, $n, 1);     # N = N + 1
  }
  # N is now the next odd number
  while (1) {
    return $n if is_prime($n);  # check primality of that number, return if prime
    Rmpz_add_ui($n, $n, 2);     # N = N + 2
  }
}

sub prev_prime($) {
  my $n = GMP->new($_[0]);
  my $cmp = Rmpz_cmp_ui($n, 3);   # compare N with 3
  if ($cmp == 0) {                # N = 3
    return GMP->new(2);
  } elsif ($cmp < 0) {            # N < 3
    return undef;
  } else {
    if (Rmpz_odd_p($n)) {         # if N is odd
      Rmpz_sub_ui($n, $n, 2);     # N = N - 2
    } else {
      Rmpz_sub_ui($n, $n, 1);     # N = N - 1
    }
    # N is now the previous odd number
    while (1) {
      return $n if is_prime($n);  # check primality of that number, return if prime
      Rmpz_sub_ui($n, $n, 2);     # N = N - 2
    } 
  }
}

sub prime_count($) {
  my $n = GMP->new($_[0]); # check if $n needs to be a Math::GMPz object
  my $primes = 0;
  return 0 if $n <= 1;

  for (my $i = GMP->new(0); Rmpz_cmp($i, $n) <= 0; Rmpz_add_ui($i, $i, 1)) {
    $primes++ if is_prime($i);
  }
  return $primes;
}

exp(0); # End of Math::Primality