File Coverage

blib/lib/Math/Prime/Util.pm

Criterion	Covered	Total	%
statement	384	508	75.5
branch	208	366	56.8
condition	69	159	43.4
subroutine	58	62	93.5
pod	27	27	100.0
total	746	1122	66.4

package Math::Prime::Util;
use strict;
use warnings;
use Carp qw/croak confess carp/;

BEGIN {
  $Math::Prime::Util::AUTHORITY = 'cpan:DANAJ';
  $Math::Prime::Util::VERSION = '0.68';
}

# parent is cleaner, and in the Perl 5.10.1 / 5.12.0 core, but not earlier.
# use parent qw( Exporter );
use base qw( Exporter );
our @EXPORT_OK =
  qw( prime_get_config prime_set_config
      prime_precalc prime_memfree
      is_prime is_prob_prime is_provable_prime is_provable_prime_with_cert
      prime_certificate verify_prime
      is_pseudoprime is_euler_pseudoprime is_strong_pseudoprime
      is_euler_plumb_pseudoprime
      is_lucas_pseudoprime
      is_strong_lucas_pseudoprime
      is_extra_strong_lucas_pseudoprime
      is_almost_extra_strong_lucas_pseudoprime
      is_frobenius_pseudoprime
      is_frobenius_underwood_pseudoprime is_frobenius_khashin_pseudoprime
      is_perrin_pseudoprime is_catalan_pseudoprime
      is_aks_prime is_bpsw_prime is_ramanujan_prime is_mersenne_prime
      is_power is_prime_power is_pillai is_semiprime is_square is_polygonal
      is_square_free is_primitive_root is_carmichael is_quasi_carmichael
      is_fundamental
      sqrtint rootint logint
      miller_rabin_random
      lucas_sequence lucasu lucasv
      primes twin_primes ramanujan_primes sieve_prime_cluster sieve_range
      forprimes forcomposites foroddcomposites fordivisors
      forpart forcomp forcomb forperm forderange formultiperm
      lastfor
      numtoperm permtonum randperm shuffle
      prime_iterator prime_iterator_object
      next_prime  prev_prime
      prime_count
      prime_count_lower prime_count_upper prime_count_approx
      nth_prime nth_prime_lower nth_prime_upper nth_prime_approx inverse_li
      twin_prime_count twin_prime_count_approx
      nth_twin_prime nth_twin_prime_approx
      ramanujan_prime_count ramanujan_prime_count_approx
      ramanujan_prime_count_lower ramanujan_prime_count_upper
      nth_ramanujan_prime nth_ramanujan_prime_approx
      nth_ramanujan_prime_lower nth_ramanujan_prime_upper
      sum_primes print_primes
      random_prime random_ndigit_prime random_nbit_prime random_strong_prime
      random_proven_prime random_proven_prime_with_cert
      random_maurer_prime random_maurer_prime_with_cert
      random_shawe_taylor_prime random_shawe_taylor_prime_with_cert
      random_semiprime random_unrestricted_semiprime
      primorial pn_primorial consecutive_integer_lcm gcdext chinese
      gcd lcm factor factor_exp divisors valuation hammingweight
      todigits fromdigits todigitstring sumdigits
      invmod sqrtmod addmod mulmod divmod powmod
      vecsum vecmin vecmax vecprod vecreduce vecextract
      vecany vecall vecnotall vecnone vecfirst vecfirstidx
      moebius mertens euler_phi jordan_totient exp_mangoldt liouville
      partitions bernfrac bernreal harmfrac harmreal
      chebyshev_theta chebyshev_psi
      divisor_sum carmichael_lambda kronecker hclassno
      ramanujan_tau ramanujan_sum
      binomial stirling znorder znprimroot znlog legendre_phi
      factorial factorialmod
      ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR LambertW Pi
      irand irand64 drand urandomb urandomm csrand random_bytes entropy_bytes
  );
our %EXPORT_TAGS = (all  => [ @EXPORT_OK ],
                    rand => [qw/srand rand irand irand64/],
                   );

# These are only exported if specifically asked for
push @EXPORT_OK, (qw/trial_factor fermat_factor holf_factor lehman_factor squfof_factor prho_factor pbrent_factor pminus1_factor pplus1_factor ecm_factor rand srand/);

my %_Config;
my %_GMPfunc;  # List of available MPU::GMP functions

# Similar to how boolean handles its option
sub import {
    if ($] < 5.020) {   # Prevent "used only once" warnings
      my $pkg = caller;
      no strict 'refs';  ## no critic(strict)
      ${"${pkg}::a"} = ${"${pkg}::a"};
      ${"${pkg}::b"} = ${"${pkg}::b"};
    }
    foreach my $opt (qw/nobigint secure/) {
      my @options = grep $_ ne "-$opt", @_;
      $_Config{$opt} = 1 if @options != @_;
      @_ = @options;
    }
    _XS_set_secure() if $_Config{'xs'} && $_Config{'secure'};
    goto &Exporter::import;
}

#############################################################################


BEGIN {

  # Separate lines to keep compatible with default from 5.6.2.
  # We could alternately use Config's $Config{uvsize} for MAXBITS
  use constant OLD_PERL_VERSION=> $] < 5.008;
  use constant MPU_MAXBITS     => (~0 == 4294967295) ? 32 : 64;
  use constant MPU_32BIT       => MPU_MAXBITS == 32;
  use constant MPU_MAXPARAM    => MPU_32BIT ? 4294967295 : 18446744073709551615;
  use constant MPU_MAXDIGITS   => MPU_32BIT ?         10 : 20;
  use constant MPU_MAXPRIME    => MPU_32BIT ? 4294967291 : 18446744073709551557;
  use constant MPU_MAXPRIMEIDX => MPU_32BIT ?  203280221 :   425656284035217743;
  use constant UVPACKLET       => MPU_32BIT ?        'L' : 'Q';
  use constant INTMAX          => (!OLD_PERL_VERSION || MPU_32BIT) ? ~0 : 562949953421312;

  eval {
    return 0 if defined $ENV{MPU_NO_XS} && $ENV{MPU_NO_XS} == 1;
    require XSLoader;
    XSLoader::load(__PACKAGE__, $Math::Prime::Util::VERSION);
    prime_precalc(0);
    $_Config{'maxbits'} = _XS_prime_maxbits();
    $_Config{'xs'} = 1;
    1;
  } or do {
    carp "Using Pure Perl implementation: $@"
      unless defined $ENV{MPU_NO_XS} && $ENV{MPU_NO_XS} == 1;

    $_Config{'xs'} = 0;
    $_Config{'maxbits'} = MPU_MAXBITS;

    # Load PP front end code
    require Math::Prime::Util::PPFE;

    # Init rand
    Math::Prime::Util::csrand();

    *prime_count   = \&Math::Prime::Util::_generic_prime_count;
    *factor        = \&Math::Prime::Util::_generic_factor;
    *factor_exp    = \&Math::Prime::Util::_generic_factor_exp;
  };

  $_Config{'secure'} = 0;
  $_Config{'nobigint'} = 0;
  $_Config{'gmp'} = 0;
  # See if they have the GMP module and haven't requested it not to be used.
  if (!defined $ENV{MPU_NO_GMP} || $ENV{MPU_NO_GMP} != 1) {
    if (eval { require Math::Prime::Util::GMP;
               Math::Prime::Util::GMP->import();
               1; }) {
      $_Config{'gmp'} = int(100*$Math::Prime::Util::GMP::VERSION);
    }
    for my $e (@Math::Prime::Util::GMP::EXPORT_OK) {
      $Math::Prime::Util::_GMPfunc{"$e"} = $_Config{'gmp'};
    }
    # If we have GMP, it is not seeded properly but we are, seed with our data.
    if (   $_Config{'gmp'} >= 42
        && !Math::Prime::Util::GMP::is_csprng_well_seeded()
        && Math::Prime::Util::_is_csprng_well_seeded()) {
      Math::Prime::Util::GMP::seed_csprng(256, random_bytes(256));
    }
  }
}

croak "Perl and XS don't agree on bit size"
      if $_Config{'xs'} && MPU_MAXBITS != _XS_prime_maxbits();

$_Config{'maxparam'}    = MPU_MAXPARAM;
$_Config{'maxdigits'}   = MPU_MAXDIGITS;
$_Config{'maxprime'}    = MPU_MAXPRIME;
$_Config{'maxprimeidx'} = MPU_MAXPRIMEIDX;
$_Config{'assume_rh'}   = 0;
$_Config{'verbose'}     = 0;

# used for code like:
#    return _XS_foo($n)  if $n <= $_XS_MAXVAL
# which builds into one scalar whether XS is available and if we can call it.
my $_XS_MAXVAL = $_Config{'xs'}  ?  MPU_MAXPARAM  :  -1;
my $_HAVE_GMP = $_Config{'gmp'};
_XS_set_callgmp($_HAVE_GMP) if $_Config{'xs'};

# Infinity in Perl is rather O/S specific.
our $_Infinity = 0+'inf';
$_Infinity = 20**20**20 if 65535 > $_Infinity;   # E.g. Windows
our $_Neg_Infinity = -$_Infinity;

sub prime_get_config {
  my %config = %_Config;

  $config{'precalc_to'} = ($_Config{'xs'})
                        ? _get_prime_cache_size()
                        : Math::Prime::Util::PP::_get_prime_cache_size();

  return \%config;
}

# Note: You can cause yourself pain if you turn on xs or gmp when they're not
# loaded.  Your calls will probably die horribly.
sub prime_set_config {
  my %params = (@_);  # no defaults
  foreach my $param (keys %params) {
    my $value = $params{$param};
    $param = lc $param;
    # dispatch table should go here.
    if      ($param eq 'xs') {
      $_Config{'xs'} = ($value) ? 1 : 0;
      $_XS_MAXVAL = $_Config{'xs'}  ?  MPU_MAXPARAM  :  -1;
    } elsif ($param eq 'gmp') {
      $_HAVE_GMP = ($value) ? int(100*$Math::Prime::Util::GMP::VERSION) : 0;
      $_Config{'gmp'} = $_HAVE_GMP;
      $Math::Prime::Util::_GMPfunc{$_} = $_HAVE_GMP
        for keys %Math::Prime::Util::_GMPfunc;
      _XS_set_callgmp($_HAVE_GMP) if $_Config{'xs'};
    } elsif ($param eq 'nobigint') {
      $_Config{'nobigint'} = ($value) ? 1 : 0;
    } elsif ($param eq 'secure') {
      croak "Cannot disable secure once set" if !$value && $_Config{'secure'};
      if ($value) {
        $_Config{'secure'} = 1;
        _XS_set_secure() if $_Config{'xs'};
      }
    } elsif ($param eq 'irand') {
      carp "ntheory irand option is deprecated";
    } elsif ($param eq 'use_primeinc') {
      carp "ntheory use_primeinc option is deprecated";
    } elsif ($param =~ /^(assume[_ ]?)?[ge]?rh$/ || $param =~ /riemann\s*h/) {
      $_Config{'assume_rh'} = ($value) ? 1 : 0;
    } elsif ($param eq 'verbose') {
      if    ($value =~ /^\d+$/) { }
      elsif ($value =~ /^[ty]/i) { $value = 1; }
      elsif ($value =~ /^[fn]/i) { $value = 0; }
      else { croak("Invalid setting for verbose.  0, 1, 2, etc."); }
      $_Config{'verbose'} = $value;
      _XS_set_verbose($value) if $_Config{'xs'};
      Math::Prime::Util::GMP::_GMP_set_verbose($value) if $_Config{'gmp'};
    } else {
      croak "Unknown or invalid configuration setting: $param\n";
    }
  }
  1;
}

sub _bigint_to_int {
  return (OLD_PERL_VERSION) ? unpack(UVPACKLET,pack(UVPACKLET,$_[0]->bstr))
                            : int($_[0]->bstr);
}
sub _to_bigint {
  do { require Math::BigInt;  Math::BigInt->import(try=>"GMP,Pari"); }
    unless defined $Math::BigInt::VERSION;
  return (!defined($_[0]) || ref($_[0]) eq 'Math::BigInt') ? $_[0] : Math::BigInt->new("$_[0]");
}
sub _to_gmpz {
  do { require Math::GMPz; } unless defined $Math::GMPz::VERSION;
  return (ref($_[0]) eq 'Math::GMPz') ? $_[0] : Math::GMPz->new($_[0]);
}
sub _to_gmp {
  do { require Math::GMP; } unless defined $Math::GMP::VERSION;
  return (ref($_[0]) eq 'Math::GMP') ? $_[0] : Math::GMP->new($_[0]);
}
sub _reftyped {
  return unless defined $_[1];
  my $ref0 = ref($_[0]);
  if ($ref0) {
    return  ($ref0 eq ref($_[1])) ?  $_[1]  :  $ref0->new("$_[1]");
  }
  if (OLD_PERL_VERSION) {
    # Perl 5.6 truncates arguments to doubles if you look at them funny
    return "$_[1]" if "$_[1]" <= INTMAX;
  } elsif ($_[1] >= 0) {
    # TODO: This wasn't working right in 5.20.0-RC1, verify correct
    return $_[1] if $_[1] <= ~0;
  } else {
    return $_[1] if ''.int($_[1]) >= -(~0>>1);
  }
  do { require Math::BigInt;  Math::BigInt->import(try=>"GMP,Pari"); }
    unless defined $Math::BigInt::VERSION;
  return Math::BigInt->new("$_[1]");
}


#*_validate_positive_integer = \&Math::Prime::Util::PP::_validate_positive_integer;
sub _validate_positive_integer {
  my($n, $min, $max) = @_;
  croak "Parameter must be defined" if !defined $n;
  if (ref($n) eq 'CODE') {
    $_[0] = $_[0]->();
    $n = $_[0];
  }
  if (ref($n) eq 'Math::BigInt') {
    croak "Parameter '$n' must be a positive integer"
      if $n->sign() ne '+' || !$n->is_int();
    $_[0] = _bigint_to_int($_[0]) if $n <= '' . INTMAX;
  } elsif (ref($n) eq 'Math::GMPz') {
    croak "Parameter '$n' must be a positive integer" if Math::GMPz::Rmpz_sgn($n) < 0;
    $_[0] = _bigint_to_int($_[0]) if $n <= INTMAX;
  } else {
    my $strn = "$n";
    croak "Parameter '$strn' must be a positive integer"
      if $strn eq '' || ($strn =~ tr/0123456789//c && $strn !~ /^\+?\d+$/);
    if ($n <= INTMAX) {
      $_[0] = $strn if ref($n);
    } else {
      #$_[0] = Math::BigInt->new($strn)
      $_[0] = _to_bigint($strn);
    }
  }
  $_[0]->upgrade(undef) if ref($_[0]) eq 'Math::BigInt' && $_[0]->upgrade();
  croak "Parameter '$_[0]' must be >= $min" if defined $min && $_[0] < $min;
  croak "Parameter '$_[0]' must be <= $max" if defined $max && $_[0] > $max;
  1;
}

#############################################################################

sub primes {
  my($low,$high) = @_;
  if (scalar @_ > 1) {
    _validate_num($low) || _validate_positive_integer($low);
  } else {
    ($low,$high) = (2, $low);
  }
  _validate_num($high) || _validate_positive_integer($high);

  my $sref = [];
  return $sref if ($low > $high) || ($high < 2);

  if ($high > $_XS_MAXVAL) {
    if ($_HAVE_GMP) {
      $sref = Math::Prime::Util::GMP::primes($low,$high);
      if ($high > ~0) {
        # Convert the returned strings into BigInts
        @$sref = map { _reftyped($_[0],$_) } @$sref;
      } else {
        @$sref = map { int($_) } @$sref;
      }
      return $sref;
    }
    require Math::Prime::Util::PP;
    return Math::Prime::Util::PP::primes($low,$high);
  }

  # Decide the method to use.  We have four to choose from:
  #  1. Trial     No memory, no overhead, but more time per prime.
  #  2. Sieve     Monolithic cached sieve.
  #  3. Erat      Monolithic uncached sieve.
  #  4. Segment   Segment sieve.  Never a bad decision.

  if (($low+1) >= $high ||                      # Tiny range, or
      $high > 10**14 && ($high-$low) < 50000) { # Small relative range

      $sref = trial_primes($low, $high);

  } elsif ($high <= (65536*30) ||                # Very small, or
           $high <= _get_prime_cache_size()) {   # already in the main cache.

      $sref = sieve_primes($low, $high);

  } else {

      $sref = segment_primes($low, $high);

  }

  # We could return an array ref in scalar context, array in list context with:
  #   return (wantarray) ? @{$sref} : $sref;
  # but I think the dual interface could be confusing, albeit often handy.
  return $sref;
}

# Shortcut for primes returning an array instead of array reference.
# sub aprimes { @{primes(@_)}; }

sub twin_primes {
  my($low,$high) = @_;
  if (scalar @_ > 1) {
    _validate_num($low) || _validate_positive_integer($low);
  } else {
    ($low,$high) = (2, $low);
  }
  _validate_num($high) || _validate_positive_integer($high);

  return [] if ($low > $high) || ($high < 2);

  if ($high > $_XS_MAXVAL) {
    my @tp;
    if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::sieve_twin_primes && $low > 2**31) {
      @tp = map { _reftyped($_[0],$_) } Math::Prime::Util::GMP::sieve_twin_primes($low, $high);
    } else {
      require Math::Prime::Util::PP;
      @tp = map { _reftyped($_[0],$_) }  Math::Prime::Util::PP::sieve_prime_cluster($low,$high, 2);
    }
    return \@tp;
  }

  return segment_twin_primes($low, $high);
}

sub ramanujan_primes {
  my($low,$high) = @_;
  if (scalar @_ > 1) {
    _validate_num($low) || _validate_positive_integer($low);
  } else {
    ($low,$high) = (2, $low);
  }
  _validate_num($high) || _validate_positive_integer($high);

  return [] if ($low > $high) || ($high < 2);

  if ($high > $_XS_MAXVAL) {
    require Math::Prime::Util::PP;
    return Math::Prime::Util::PP::_ramanujan_primes($low,$high);
  }
  return _ramanujan_primes($low, $high);
}

#############################################################################
# Random primes.  These are front end functions that do input validation,
# load the RandomPrimes module, and call its function.

sub random_maurer_prime_with_cert {
  my($bits) = @_;
  _validate_num($bits, 2) || _validate_positive_integer($bits, 2);

  if ($Math::Prime::Util::_GMPfunc{"random_maurer_prime_with_cert"}) {
    my($n,$cert) = Math::Prime::Util::GMP::random_maurer_prime_with_cert($bits);
    return (Math::Prime::Util::_reftyped($_[0],$n), $cert);
  }

  require Math::Prime::Util::RandomPrimes;
  return Math::Prime::Util::RandomPrimes::random_maurer_prime_with_cert($bits);
}

sub random_shawe_taylor_prime_with_cert {
  my($bits) = @_;
  _validate_num($bits, 2) || _validate_positive_integer($bits, 2);

  if ($Math::Prime::Util::_GMPfunc{"random_shawe_taylor_prime_with_cert"}) {
    my($n,$cert) =Math::Prime::Util::GMP::random_shawe_taylor_prime_with_cert($bits);
    return (Math::Prime::Util::_reftyped($_[0],$n), $cert);
  }

  require Math::Prime::Util::RandomPrimes;
  return Math::Prime::Util::RandomPrimes::random_shawe_taylor_prime_with_cert($bits);
}

sub random_proven_prime_with_cert {
  my($bits) = @_;
  _validate_num($bits, 2) || _validate_positive_integer($bits, 2);

  # Go to Maurer with GMP
  if ($Math::Prime::Util::_GMPfunc{"random_maurer_prime_with_cert"}) {
    my($n,$cert) = Math::Prime::Util::GMP::random_maurer_prime_with_cert($bits);
    return (Math::Prime::Util::_reftyped($_[0],$n), $cert);
  }

  require Math::Prime::Util::RandomPrimes;
  return Math::Prime::Util::RandomPrimes::random_proven_prime_with_cert($bits);
}


#############################################################################
# These functions almost always return bigints, so there is no XS
# implementation.  Try to run the GMP version, and if it isn't available,
# load PP and call it.

sub primorial {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);
  return (1,1,2,6,6,30,30,210,210,210)[$n] if $n < 10;

  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::primorial) {
    return _reftyped($_[0], Math::Prime::Util::GMP::primorial($n));
  }
  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::primorial($n);
}

sub pn_primorial {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);
  return (1,2,6,30,210,2310,30030,510510,9699690,223092870)[$n] if $n < 10;

  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::pn_primorial) {
    return _reftyped($_[0], Math::Prime::Util::GMP::pn_primorial($n));
  }

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::primorial(nth_prime($n));
}

sub consecutive_integer_lcm {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);
  return 0 if $n < 1;

  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::consecutive_integer_lcm) {
    return _reftyped($_[0],Math::Prime::Util::GMP::consecutive_integer_lcm($n));
  }
  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::consecutive_integer_lcm($n);
}

# See 2011+ FLINT and Fredrik Johansson's work for state of the art.
# Very crude timing estimates (ignores growth rates).
#   Perl-comb   partitions(10^5)  ~ 300 seconds  ~200,000x slower than Pari
#   GMP-comb    partitions(10^6)  ~ 120 seconds    ~1,000x slower than Pari
#   Pari        partitions(10^8)  ~ 100 seconds
#   Bober 0.6   partitions(10^9)  ~  20 seconds       ~50x faster than Pari
#   Johansson   partitions(10^12) ~  10 seconds     >1000x faster than Pari
sub partitions {
  my($n) = @_;
  return 1 if defined $n && $n <= 0;
  _validate_num($n) || _validate_positive_integer($n);

  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::partitions) {
    return _reftyped($_[0],Math::Prime::Util::GMP::partitions($n));
  }

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::partitions($n);
}

#############################################################################
# forprimes, forcomposites, fordivisors.
# These are used when the XS code can't handle it.

sub _generic_forprimes {
  my($sub, $beg, $end) = @_;
  if (!defined $end) { $end = $beg; $beg = 2; }
  _validate_positive_integer($beg);
  _validate_positive_integer($end);
  $beg = 2 if $beg < 2;
  my $oldforexit = Math::Prime::Util::_start_for_loop();
  {
    my $pp;
    local *_ = \$pp;
    for (my $p = next_prime($beg-1);  $p <= $end;  $p = next_prime($p)) {
      $pp = $p;
      $sub->();
      last if Math::Prime::Util::_get_forexit();
    }
  }
  Math::Prime::Util::_end_for_loop($oldforexit);
}

sub _generic_forcomposites {
  my($sub, $beg, $end) = @_;
  if (!defined $end) { $end = $beg; $beg = 4; }
  _validate_positive_integer($beg);
  _validate_positive_integer($end);
  $beg = 4 if $beg < 4;
  $end = Math::BigInt->new(''.~0) if ref($end) ne 'Math::BigInt' && $end == ~0;
  my $oldforexit = Math::Prime::Util::_start_for_loop();
  {
    my $pp;
    local *_ = \$pp;
    for (my $p = next_prime($beg-1);  $beg <= $end;  $p = next_prime($p)) {
      for ( ; $beg < $p && $beg <= $end ; $beg++ ) {
        $pp = $beg;
        $sub->();
        last if Math::Prime::Util::_get_forexit();
      }
      $beg++;
      last if Math::Prime::Util::_get_forexit();
    }
  }
  Math::Prime::Util::_end_for_loop($oldforexit);
}

sub _generic_foroddcomposites {
  my($sub, $beg, $end) = @_;
  if (!defined $end) { $end = $beg; $beg = 9; }
  _validate_positive_integer($beg);
  _validate_positive_integer($end);
  $beg = 9 if $beg < 9;
  $beg++ unless $beg & 1;
  $end = Math::BigInt->new(''.~0) if ref($end) ne 'Math::BigInt' && $end == ~0;
  my $oldforexit = Math::Prime::Util::_start_for_loop();
  {
    my $pp;
    local *_ = \$pp;
    for (my $p = next_prime($beg-1);  $beg <= $end;  $p = next_prime($p)) {
      for ( ; $beg < $p && $beg <= $end ; $beg += 2 ) {
        $pp = $beg;
        $sub->();
        last if Math::Prime::Util::_get_forexit();
      }
      $beg += 2;
      last if Math::Prime::Util::_get_forexit();
    }
  }
  Math::Prime::Util::_end_for_loop($oldforexit);
}

sub _generic_fordivisors {
  my($sub, $n) = @_;
  _validate_positive_integer($n);
  my @divisors = divisors($n);
  my $oldforexit = Math::Prime::Util::_start_for_loop();
  {
    my $pp;
    local *_ = \$pp;
    foreach my $d (@divisors) {
      $pp = $d;
      $sub->();
      last if Math::Prime::Util::_get_forexit();
    }
  }
  Math::Prime::Util::_end_for_loop($oldforexit);
}

sub formultiperm (&$) {    ## no critic qw(ProhibitSubroutinePrototypes)
  my($sub, $iref) = @_;
  croak("formultiperm first argument must be an array reference") unless ref($iref) eq 'ARRAY';

  my($sum, %h, @n) = (0);
  $h{$_}++ for @$iref;
  @n = map { [$_, $h{$_}] } sort(keys(%h));
  $sum += $_->[1] for @n;

  require Math::Prime::Util::PP;
  my $oldforexit = Math::Prime::Util::_start_for_loop();
  Math::Prime::Util::PP::_multiset_permutations( $sub, [], \@n, $sum );
  Math::Prime::Util::_end_for_loop($oldforexit);
}

#############################################################################
# Iterators

sub prime_iterator {
  my($start) = @_;
  $start = 0 unless defined $start;
  _validate_num($start) || _validate_positive_integer($start);
  my $p = ($start > 0) ? $start-1 : 0;
  # This works fine:
  #   return sub { $p = next_prime($p); return $p; };
  # but we can optimize a little
  if (ref($p) ne 'Math::BigInt' && $p <= $_XS_MAXVAL) {
    # This is simple and low memory, but slower than segments:
    #     return sub { $p = next_prime($p); return $p; };
    my $pr = [];
    return sub {
      if (scalar(@$pr) == 0) {
        # Once we're into bigints, just use next_prime
        return $p=next_prime($p) if $p >= MPU_MAXPRIME;
        # Get about 10k primes
        my $segment = ($p <= 1e4) ? 10_000 : int(10000*log($p)+1);
        $segment = ~0-$p if $p+$segment > ~0 && $p+1 < ~0;
        $pr = primes($p+1, $p+$segment);
      }
      return $p = shift(@$pr);
    };
  } elsif ($_HAVE_GMP) {
    return sub { $p = $p-$p+Math::Prime::Util::GMP::next_prime($p); return $p;};
  } else {
    require Math::Prime::Util::PP;
    return sub { $p = Math::Prime::Util::PP::next_prime($p); return $p; }
  }
}

sub prime_iterator_object {
  my($start) = @_;
  require Math::Prime::Util::PrimeIterator;
  return Math::Prime::Util::PrimeIterator->new($start);
}

#############################################################################
# Front ends to functions.
#
# These will do input validation, then call the appropriate internal function
# based on the input (XS, GMP, PP).
#############################################################################

sub _generic_prime_count {
  my($low,$high) = @_;
  if (defined $high) {
    _validate_num($low) || _validate_positive_integer($low);
    _validate_num($high) || _validate_positive_integer($high);
  } else {
    ($low,$high) = (2, $low);
    _validate_num($high) || _validate_positive_integer($high);
  }
  return 0 if $high < 2  ||  $low > $high;

  # We can relax these constraints if MPU::GMP gets a fast implementation.
  return Math::Prime::Util::GMP::prime_count($low,$high) if $_HAVE_GMP
                       && defined &Math::Prime::Util::GMP::prime_count
                       && (   (ref($high) eq 'Math::BigInt')
                           || (($high-$low) < int($low/1_000_000))
                          );
  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::prime_count($low,$high);
}

sub _generic_factor {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);

  if ($_HAVE_GMP) {
    my @factors;
    if ($n != 1) {
      @factors = Math::Prime::Util::GMP::factor($n);
      #if (ref($_[0]) eq 'Math::BigInt') {
      #  @factors = map { ($_ > ~0) ? Math::BigInt->new(''.$_) : $_ } @factors;
      #}
      if (ref($_[0])) {
        @factors = map { ($_ > ~0) ? ref($_[0])->new(''.$_) : $_ } @factors;
      }
    }
    return @factors;
  }

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::factor($n);
}

sub _generic_factor_exp {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);

  my %exponents;
  my @factors = grep { !$exponents{$_}++ } factor($n);
  return scalar @factors unless wantarray;
  return (map { [$_, $exponents{$_}] } @factors);
}

#############################################################################

sub _is_gaussian_prime {
  my($a,$b) = @_;
  return ((($b % 4) == 3) ? is_prime($b) : 0) if $a == 0;
  return ((($a % 4) == 3) ? is_prime($a) : 0) if $b == 0;
  is_prime( vecsum( vecprod($a,$a), vecprod($b,$b) ) );
}

#############################################################################

# Return just the cert portion.
sub prime_certificate {
  my($n) = @_;
  my ($is_prime, $cert) = is_provable_prime_with_cert($n);
  return $cert;
}


sub is_provable_prime_with_cert {
  my($n) = @_;
  return 0 if defined $n && $n < 2;
  _validate_num($n) || _validate_positive_integer($n);
  my $header = "[MPU - Primality Certificate]\nVersion 1.0\n\nProof for:\nN $n\n\n";

  if ($n <= $_XS_MAXVAL) {
    my $isp = is_prime($n);
    return ($isp, '') unless $isp == 2;
    return (2, $header . "Type Small\nN $n\n");
  }

  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::is_provable_prime_with_cert) {
    my ($isp, $cert) = Math::Prime::Util::GMP::is_provable_prime_with_cert($n);
    # New version that returns string format.
    #return ($isp, $cert) if ref($cert) ne 'ARRAY';
    if (ref($cert) ne 'ARRAY') {
      # Fix silly 0.13 mistake (TODO: deprecate this)
      $cert =~ s/^Type Small\n(\d+)/Type Small\nN $1/smg;
      return ($isp, $cert);
    }
    # Old version.  Convert.
    require Math::Prime::Util::PrimalityProving;
    return ($isp, Math::Prime::Util::PrimalityProving::convert_array_cert_to_string($cert));
  }

  {
    my $isp = is_prob_prime($n);
    return ($isp, '') if $isp == 0;
    return (2, $header . "Type Small\nN $n\n") if $isp == 2;
  }

  # Choice of methods for proof:
  #   ECPP         needs a fair bit of programming work
  #   APRCL        needs a lot of programming work
  #   BLS75 combo  Corollary 11 of BLS75.  Trial factor n-1 and n+1 to B, find
  #                factors F1 of n-1 and F2 of n+1.  Quit when:
  #                B > (N/(F1*F1*(F2/2)))^1/3 or B > (N/((F1/2)*F2*F2))^1/3
  #   BLS75 n+1    Requires factoring n+1 to (n/2)^1/3 (theorem 19)
  #   BLS75 n-1    Requires factoring n-1 to (n/2)^1/3 (theorem 5 or 7)
  #   Pocklington  Requires factoring n-1 to n^1/2 (BLS75 theorem 4)
  #   Lucas        Easy, requires factoring of n-1 (BLS75 theorem 1)
  #   AKS          horribly slow
  # See http://primes.utm.edu/prove/merged.html or other sources.

  require Math::Prime::Util::PrimalityProving;
  #my ($isp, $pref) = Math::Prime::Util::PrimalityProving::primality_proof_lucas($n);
  my ($isp, $pref) = Math::Prime::Util::PrimalityProving::primality_proof_bls75($n);
  carp "proved $n is not prime\n" if !$isp;
  return ($isp, $pref);
}


sub verify_prime {
  require Math::Prime::Util::PrimalityProving;
  return Math::Prime::Util::PrimalityProving::verify_cert(@_);
}

#############################################################################

sub RiemannZeta {
  my($n) = @_;
  croak("Invalid input to ReimannZeta:  x must be > 0") if $n <= 0;

  return $n-$n if $n > 10_000_000;   # Over 3M leading zeros

  return _XS_RiemannZeta($n)
  if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $_Config{'xs'};

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::RiemannZeta($n);
}

sub RiemannR {
  my($n) = @_;
  croak("Invalid input to ReimannR:  x must be > 0") if $n <= 0;

  return _XS_RiemannR($n)
  if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $_Config{'xs'};

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::RiemannR($n);
}

sub ExponentialIntegral {
  my($n) = @_;
  return $_Neg_Infinity if $n == 0;
  return 0              if $n == $_Neg_Infinity;
  return $_Infinity     if $n == $_Infinity;

  return _XS_ExponentialIntegral($n)
  if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $_Config{'xs'};

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::ExponentialIntegral($n);
}

sub LogarithmicIntegral {
  my($n) = @_;
  return 0              if $n == 0;
  return $_Neg_Infinity if $n == 1;
  return $_Infinity     if $n == $_Infinity;

  croak("Invalid input to LogarithmicIntegral:  x must be >= 0") if $n <= 0;

  if (!defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $_Config{'xs'}) {
    return 1.045163780117492784844588889194613136522615578151 if $n == 2;
    return _XS_LogarithmicIntegral($n);
  }

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::LogarithmicIntegral(@_);
}

sub LambertW {
  my($x) = @_;

  return _XS_LambertW($x)
  if !defined $bignum::VERSION && ref($x) ne 'Math::BigFloat' && $_Config{'xs'};

  # TODO: Call GMP function here directly

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::LambertW($x);
}

sub bernfrac {
  my($n) = @_;
  return map { _to_bigint($_) } (0,1) if defined $n && $n < 0;
  _validate_num($n) || _validate_positive_integer($n);
  return map { _to_bigint($_) } (0,1) if $n > 1 && ($n & 1);

  if ($Math::Prime::Util::_GMPfunc{"bernfrac"}) {
    return map { _to_bigint($_) } Math::Prime::Util::GMP::bernfrac($n);
  }

  require Math::Prime::Util::PP;
  return Math::Prime::Util::PP::bernfrac($n);
}
sub bernreal {
  my($n, $precision) = @_;
  do { require Math::BigFloat; Math::BigFloat->import(); } unless defined $Math::BigFloat::VERSION;

  if ($Math::Prime::Util::_GMPfunc{"bernreal"}) {
    return Math::BigFloat->new(Math::Prime::Util::GMP::bernreal($n)) if !defined $precision;
    return Math::BigFloat->new(Math::Prime::Util::GMP::bernreal($n,$precision),$precision);
  }

  my($num,$den) = bernfrac($n);
  return Math::BigFloat->bzero if $num->is_zero;
  scalar Math::BigFloat->new($num)->bdiv($den, $precision);
}

sub harmfrac {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);
  return map { _to_bigint($_) } (0,1) if $n <= 0;

  if ($Math::Prime::Util::_GMPfunc{"harmfrac"}) {
    return map { _to_bigint($_) } Math::Prime::Util::GMP::harmfrac($n);
  }

  require Math::Prime::Util::PP;
  Math::Prime::Util::PP::harmfrac($n);
}
sub harmreal {
  my($n, $precision) = @_;
  _validate_num($n) || _validate_positive_integer($n);
  do { require Math::BigFloat; Math::BigFloat->import(); } unless defined $Math::BigFloat::VERSION;
  return Math::BigFloat->bzero if $n <= 0;

  if ($Math::Prime::Util::_GMPfunc{"harmreal"}) {
    return Math::BigFloat->new(Math::Prime::Util::GMP::harmreal($n)) if !defined $precision;
    return Math::BigFloat->new(Math::Prime::Util::GMP::harmreal($n,$precision),$precision);
  }

  # If low enough precision, use native floating point.  Fast.
  if (defined $precision && $precision <= 13) {
    return Math::BigFloat->new(
      ($n < 80) ? do { my $h = 0; $h += 1/$_ for 1..$n; $h; }
                : log($n) + 0.57721566490153286060651209 + 1/(2*$n) - 1/(12*$n*$n) + 1/(120*$n*$n*$n*$n)
      ,$precision
    );
  }

  if ($Math::Prime::Util::_GMPfunc{"harmfrac"}) {
    my($num,$den) = map { _to_bigint($_) } Math::Prime::Util::GMP::harmfrac($n);
    return scalar Math::BigFloat->new($num)->bdiv($den, $precision);
  }

  require Math::Prime::Util::PP;
  Math::Prime::Util::PP::harmreal($n, $precision);
}

#############################################################################

use Math::Prime::Util::MemFree;

1;

__END__