File Coverage

blib/lib/Math/Prime/Util.pm

Criterion	Covered	Total	%
statement	260	370	70.2
branch	142	260	54.6
condition	48	99	48.4
subroutine	47	53	88.6
pod	18	18	100.0
total	515	800	64.3

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.74';
}

# 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_delicate_prime is_chen_prime
      is_odd is_even is_divisible is_congruent
      is_power is_prime_power is_perfect_power is_square
      is_square_free is_powerfree
      is_pillai is_polygonal is_congruent_number is_perfect_number
      is_semiprime is_almost_prime is_omega_prime
      is_primitive_root is_carmichael is_quasi_carmichael is_cyclic
      is_fundamental is_totient is_gaussian_prime is_sum_of_squares
      is_smooth is_rough is_powerful is_practical is_lucky is_happy
      sqrtint rootint logint lshiftint rshiftint rashiftint absint negint
      signint cmpint addint subint add1int sub1int mulint powint
      divint modint cdivint divrem fdivrem cdivrem tdivrem
      miller_rabin_random
      lucas_sequence
      lucasu lucasv lucasuv lucasumod lucasvmod lucasuvmod pisano_period
      primes twin_primes semi_primes almost_primes omega_primes ramanujan_primes
      sieve_prime_cluster sieve_range prime_powers lucky_numbers
      forprimes forcomposites foroddcomposites forsemiprimes foralmostprimes
      forpart forcomp forcomb forperm forderange formultiperm forsetproduct
      fordivisors forfactored forsquarefree forsquarefreeint
      lastfor
      numtoperm permtonum randperm shuffle vecsample
      prime_iterator prime_iterator_object
      next_prime prev_prime
      next_prime_power prev_prime_power
      next_perfect_power prev_perfect_power
      next_chen_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 inverse_li_nv
      twin_prime_count twin_prime_count_approx
      nth_twin_prime nth_twin_prime_approx
      semiprime_count semiprime_count_approx
      nth_semiprime nth_semiprime_approx
      almost_prime_count almost_prime_count_approx
      almost_prime_count_lower almost_prime_count_upper
      nth_almost_prime nth_almost_prime_approx
      nth_almost_prime_lower nth_almost_prime_upper
      omega_prime_count nth_omega_prime
      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
      powerful_count nth_powerful sumpowerful powerful_numbers
      prime_power_count prime_power_count_approx
      prime_power_count_lower prime_power_count_upper
      nth_prime_power nth_prime_power_approx
      nth_prime_power_lower nth_prime_power_upper
      perfect_power_count perfect_power_count_approx
      perfect_power_count_lower perfect_power_count_upper
      nth_perfect_power nth_perfect_power_approx
      nth_perfect_power_lower nth_perfect_power_upper
      nth_powerfree powerfree_count powerfree_sum squarefree_kernel
      powerfree_part powerfree_part_sum
      smooth_count rough_count powersum
      lucky_count lucky_count_approx lucky_count_lower lucky_count_upper
      nth_lucky nth_lucky_approx nth_lucky_lower nth_lucky_upper
      minimal_goldbach_pair goldbach_pairs goldbach_pair_count
      sum_primes print_primes
      random_prime random_ndigit_prime
      random_nbit_prime random_safe_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
      random_factored_integer
      primorial pn_primorial consecutive_integer_lcm gcdext chinese chinese2
      gcd lcm factor factor_exp divisors valuation hammingweight
      frobenius_number
      todigits fromdigits todigitstring sumdigits
      tozeckendorf fromzeckendorf
      sqrtmod allsqrtmod rootmod allrootmod cornacchia
      negmod invmod addmod submod mulmod divmod powmod muladdmod mulsubmod
      vecsum vecmin vecmax vecprod vecreduce vecextract vecequal vecuniq
      vecany vecall vecnotall vecnone vecfirst vecfirstidx vecmex vecpmex
      vecsort vecsorti vecfreq vecsingleton vecslide
      setbinop sumset toset
      setunion setintersect setminus setdelta
      setcontains setcontainsany setinsert setremove setinvert
      is_sidon_set is_sumfree_set
      set_is_disjoint set_is_equal set_is_proper_intersection
      set_is_subset set_is_proper_subset set_is_superset set_is_proper_superset
      moebius mertens liouville sumliouville prime_omega prime_bigomega
      euler_phi jordan_totient exp_mangoldt sumtotient
      partitions bernfrac bernreal harmfrac harmreal
      chebyshev_theta chebyshev_psi
      divisor_sum carmichael_lambda hclassno inverse_totient
      kronecker is_qr qnr
      ramanujan_tau ramanujan_sum
      stirling fubini znorder znprimroot znlog legendre_phi
      factorial factorialmod subfactorial binomial binomialmod
      falling_factorial rising_factorial
      contfrac from_contfrac
      next_calkin_wilf next_stern_brocot
      calkin_wilf_n stern_brocot_n
      nth_calkin_wilf nth_stern_brocot
      nth_stern_diatomic
      farey next_farey farey_rank
      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 cheb_factor ecm_factor rand srand/);

my %_Config;
my %_GMPfunc;  # 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;
  use constant INTMIN          => -(INTMAX >> 1) - 1;

  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();
  };

  $_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 + 1e-6);
    }
    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));
    }
  }

  # Alias PP and GMP if requested.  Very convenient but a big hammer.
  if (defined $ENV{MPU_DEVNAMES} && $ENV{MPU_DEVNAMES} == 1) {
    no strict 'refs';  ## no critic(strict)
    *MPU:: = \*Math::Prime::Util::;
    *PP:: = \*Math::Prime::Util::PP:: if eval { require Math::Prime::Util::PP; Math::Prime::Util::PP->import(); 1; };
    *GMP:: = \*Math::Prime::Util::GMP:: if $_Config{'gmp'};
  }
}

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;
$_Config{'bigintclass'} = undef;

# 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'};

our $_BIGINT = $_Config{'bigintclass'};

# 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 'bigint' || $param eq 'trybigint') {
      my $class = _load_bigint_class($value);
      if (defined $class) {
        $_BIGINT = $_Config{'bigintclass'} = $class;
      } else {
        carp "ntheory could not load bigint class from '$value'"
          unless $param =~ /try/;
      }
    } 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;
}

# Input:  object, or comma separated list of class names
# Output: class name or undef
sub _load_bigint_class {
  my($val) = @_;

  my $class = undef;
  if (ref($val)) {      # We are given an object, e.g. a Math::GMPz number
    $class = ref($val);
  } else {              # Comma separated list of class names
    for my $name (split /,/, $val) {
      $name =~ s/^\s+|\s+$//g;
      (my $cfname="$name.pm")=~s|::|/|g; # Foo::Bar::Baz => Foo/Bar/Baz.pm
      if ($INC{$cfname} || eval { require $cfname; $name->import(); 1; }) {
        $class = $name;
        last;
      }
    }
  }
  if ($class) {         # Check we can make a number with it
    $class = undef unless eval { $class->new(1) == 1 };
  }
  return $class;
}

# This is for loading the default bigint class the very first time.
sub _load_bigint {
  return $_BIGINT if defined $_BIGINT;

  # TODO: turn this on for next release
  #prime_set_config( trybigint => 'Math::GMPz,Math::GMP' );
  #return $_BIGINT if defined $_BIGINT;

  do { require Math::BigInt;  Math::BigInt->import(try=>"GMP,GMPz,LTM,Pari"); } unless defined $Math::BigInt::VERSION;
  $_BIGINT = $_Config{'bigintclass'} = 'Math::BigInt';
}

sub _bigint_to_int {
  return (OLD_PERL_VERSION && $_[0] >= 0) ? unpack(UVPACKLET,pack(UVPACKLET,"$_[0]"))
                            : int("$_[0]");
}
sub _to_bigint {
  return undef unless defined($_[0]);
  _load_bigint() unless defined $_BIGINT;
  # We don't do any validation other than that the class is happy.
  my $n;
  if (ref($_[0]) eq $_BIGINT) {
    $n = $_[0];
  } elsif (ref($_[0]) eq 'Math::BigFloat' && !$_[0]->is_int()) {
    $n = Math::BigInt->bnan;
  } elsif ($_BIGINT eq 'Math::Pari' && $_[0] =~ /^0[bx]/) {
    # Pari added support for this in 2.8, so not in Math::Pari
    do { require Math::BigInt;  Math::BigInt->import(try=>"GMP,GMPz,LTM,Pari"); } unless defined $Math::BigInt::VERSION;
    $n = Math::BigInt->new("$_[0]");
    $n = $_BIGINT->new("$n");
  } else {
    $n = $_BIGINT->new("$_[0]");
  }
  croak "Parameter '$_[0]' must be an integer" unless $_BIGINT ne 'Math::BigInt' || $n->is_int();
  $n;
}
sub _to_bigint_nonneg {
  return undef unless defined($_[0]);
  _load_bigint() unless defined $_BIGINT;
  my $n;
  if (ref($_[0]) eq $_BIGINT) {
    $n = $_[0];
  } elsif (ref($_[0]) eq 'Math::BigFloat' && !$_[0]->is_int()) {
    $n = Math::BigInt->bnan;
  } else {
    $n = $_BIGINT->new("$_[0]");
  }
  croak "Parameter '$_[0]' must be a non-negative integer" unless ($_BIGINT ne 'Math::BigInt' || $n->is_int()) && $n >= 0;
  $n;
}
sub _to_bigint_abs {
  return undef unless defined($_[0]);
  my $n = _to_bigint($_[0]);
  return ($n < 0) ? -$n : $n;
}
sub _to_bigint_if_needed {
  return $_[0] if !defined $_[0] || ref($_[0]);
  if ($_[0] >= INTMAX || $_[0] <= INTMIN) {    # Probably a bigint
    my $n = _to_bigint($_[0]);
    return $n if $n > INTMAX || $n < INTMIN;   # Definitely a bigint
  }
  $_[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 undef unless defined $_[1];
  my $ref0 = ref($_[0]);
  if (OLD_PERL_VERSION) {
    # Perl 5.6 truncates arguments to doubles if you look at them funny
    return "$_[1]" if "$_[1]" <= INTMAX && "$_[1]" >= INTMIN;
  } elsif ($_[1] >= 0) {
    return $_[1] if $_[1] < ~0;
  } else {
    return $_[1] if $_[1] > -(~0>>1);
  }
  my $bign;
  if ($ref0) {
    $bign = $ref0->new("$_[1]");
  } else {
    _load_bigint() unless defined $_BIGINT;
    $bign = $_BIGINT->new("$_[1]");
  }
  return $bign if $bign > INTMAX || $bign < INTMIN;
  0+"$_[1]";
}

sub _maybe_bigint_allargs {
  _load_bigint() unless defined $_BIGINT;
  for my $i (0..$#_) {
    next if !defined $_[$i] || ref($_[$i]);
    next if $_[$i] < INTMAX && $_[$i] > INTMIN;
    my $n = $_BIGINT->new("$_[$i]");
    $_[$i] = $n if $n > INTMAX || $n < INTMIN;
  }
  @_;
}

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

# These are called by the XS code to keep the GMP CSPRNG in sync with us.

sub _srand_p {
  my($seedval) = @_;
  return unless $_Config{'gmp'} >= 42;
  $seedval = unpack("L",entropy_bytes(4)) unless defined $seedval;
  Math::Prime::Util::GMP::seed_csprng(4, pack("L",$seedval));
  $seedval;
}

sub _csrand_p {
  my($str) = @_;
  return unless $_Config{'gmp'} >= 42;
  $str = entropy_bytes(256) unless defined $str;
  Math::Prime::Util::GMP::seed_csprng(length($str), $str);
}

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


#############################################################################
# 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_integer_nonneg($bits);
  croak "random_maurer_prime bits must be >= 2" unless $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_integer_nonneg($bits);
  croak "random_shawe_taylor_prime bits must be >= 2" unless $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_integer_nonneg($bits);
  croak "random_proven_prime bits must be >= 2" unless $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);
}

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

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_integer_nonneg($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) && $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 = addint(0,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).
#############################################################################

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

# 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) = @_;
  _validate_integer($n);
  return 0 if $n < 2;
  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 RiemannZeta:  x must be > 0") if $n <= 0;

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

  return _XS_RiemannZeta($n) if !ref($n) && $_Config{'xs'};

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

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

  return _XS_RiemannR($n) if !ref($n) && $_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 !ref($n) && $_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 (!ref($n) && $_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 !ref($x) && $_Config{'xs'};

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

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) = map { _to_bigint($_) } bernfrac($n);
  return Math::BigFloat->bzero if $num == 0;
  scalar Math::BigFloat->new($num)->bdiv($den, $precision);
}

sub harmreal {
  my($n, $precision) = @_;
  _validate_integer_nonneg($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);
}



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

1;

__END__