File Coverage

blib/lib/Math/Prime/Util/PP.pm

Criterion	Covered	Total	%
statement	2864	4401	65.0
branch	1570	3170	49.5
condition	497	1200	41.4
subroutine	180	254	70.8
pod	4	154	2.6
total	5115	9179	55.7

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

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

BEGIN {
  do { require Math::BigInt;  Math::BigInt->import(try=>"GMP,Pari"); }
    unless defined $Math::BigInt::VERSION;
}

# The Pure Perl versions of all the Math::Prime::Util routines.
#
# Some of these will be relatively similar in performance, some will be
# very slow in comparison.
#
# Most of these are pretty simple.  Also, you really should look at the C
# code for more detailed comments, including references to papers.

BEGIN {
  use constant OLD_PERL_VERSION=> $] < 5.008;
  use constant MPU_MAXBITS     => (~0 == 4294967295) ? 32 : 64;
  use constant MPU_64BIT       => MPU_MAXBITS == 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 MPU_HALFWORD    => MPU_32BIT ? 65536 : OLD_PERL_VERSION ? 33554432 : 4294967296;
  use constant UVPACKLET       => MPU_32BIT ? 'L' : 'Q';
  use constant MPU_INFINITY    => (65535 > 0+'inf') ? 20**20**20 : 0+'inf';
  use constant CONST_EULER     => '0.577215664901532860606512090082402431042159335939923598805767';
  use constant CONST_LI2       => '1.04516378011749278484458888919461313652261557815120157583290914407501320521';
  use constant BZERO           => Math::BigInt->bzero;
  use constant BONE            => Math::BigInt->bone;
  use constant BTWO            => Math::BigInt->new(2);
  use constant INTMAX          => (!OLD_PERL_VERSION || MPU_32BIT) ? ~0 : 562949953421312;
  use constant BMAX            => Math::BigInt->new('' . INTMAX);
  use constant B_PRIM767       => Math::BigInt->new("261944051702675568529303");
  use constant B_PRIM235       => Math::BigInt->new("30");
  use constant PI_TIMES_8      => 25.13274122871834590770114707;
}

{
  my $_have_MPFR = -1;
  sub _MPFR_available {
    if ($_have_MPFR < 0) {
      $_have_MPFR = 0;
      $_have_MPFR = 1 if (!defined $ENV{MPU_NO_MPFR} || $ENV{MPU_NO_MPFR} != 1)
                      && eval { require Math::MPFR; $Math::MPFR::VERSION>=2.03; };
      # Minimum MPFR library version is 3.0 (2010).
      $_have_MPFR = 0 if $_have_MPFR && Math::MPFR::MPFR_VERSION_MAJOR() < 3;
    }
    if ($_have_MPFR && scalar(@_) == 2) {
      my($major,$minor) = @_;
      return 0 if Math::MPFR::MPFR_VERSION_MAJOR() < $major;
      return 0 if Math::MPFR::MPFR_VERSION_MAJOR() == $major && Math::MPFR::MPFR_VERSION_MINOR() < $minor;
    }
    return $_have_MPFR;
  }
}

my $_precalc_size = 0;
sub prime_precalc {
  my($n) = @_;
  croak "Parameter '$n' must be a positive integer" unless _is_positive_int($n);
  $_precalc_size = $n if $n > $_precalc_size;
}
sub prime_memfree {
  $_precalc_size = 0;
  Math::MPFR::Rmpfr_free_cache() if defined $Math::MPFR::VERSION;
}
sub _get_prime_cache_size { $_precalc_size }
sub _prime_memfreeall { prime_memfree; }


sub _is_positive_int {
  ((defined $_[0]) && $_[0] ne '' && ($_[0] !~ tr/0123456789//c));
}

sub _bigint_to_int {
  #if (OLD_PERL_VERSION) {
  #  my $pack = ($_[0] < 0) ? lc(UVPACKLET) : UVPACKLET;
  #  return unpack($pack,pack($pack,"$_[0]"));
  #}
  int("$_[0]");
}

sub _upgrade_to_float {
  do { require Math::BigFloat; Math::BigFloat->import(); }
    if !defined $Math::BigFloat::VERSION;
  Math::BigFloat->new(@_);
}

# Get the accuracy of variable x, or the max default from BigInt/BigFloat
# One might think to use ref($x)->accuracy() but numbers get upgraded and
# downgraded willy-nilly, and it will do the wrong thing from the user's
# perspective.
sub _find_big_acc {
  my($x) = @_;
  my $b;

  $b = $x->accuracy() if ref($x) =~ /^Math::Big/;
  return $b if defined $b;

  my ($i,$f) = (Math::BigInt->accuracy(), Math::BigFloat->accuracy());
  return (($i > $f) ? $i : $f)   if defined $i && defined $f;
  return $i if defined $i;
  return $f if defined $f;

  ($i,$f) = (Math::BigInt->div_scale(), Math::BigFloat->div_scale());
  return (($i > $f) ? $i : $f)   if defined $i && defined $f;
  return $i if defined $i;
  return $f if defined $f;
  return 18;
}

sub _bfdigits {
  my($wantbf, $xdigits) = (0, 17);
  if (defined $bignum::VERSION || ref($_[0]) =~ /^Math::Big/) {
    do { require Math::BigFloat; Math::BigFloat->import(); }
      if !defined $Math::BigFloat::VERSION;
    if (ref($_[0]) eq 'Math::BigInt') {
      my $xacc = ($_[0])->accuracy();
      $_[0] = Math::BigFloat->new($_[0]);
      ($_[0])->accuracy($xacc) if $xacc;
    }
    $_[0] = Math::BigFloat->new("$_[0]") if ref($_[0]) ne 'Math::BigFloat';
    $wantbf = _find_big_acc($_[0]);
    $xdigits = $wantbf;
  }
  ($wantbf, $xdigits);
}


sub _validate_num {
  my($n, $min, $max) = @_;
  croak "Parameter must be defined" if !defined $n;
  return 0 if ref($n);
  croak "Parameter '$n' must be a positive integer"
          if $n eq '' || ($n =~ tr/0123456789//c && $n !~ /^\+\d+$/);
  croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
  croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
  substr($_[0],0,1,'') if substr($n,0,1) eq '+';
  return 0 unless $n < ~0 || int($n) eq ''.~0;
  1;
}

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 <= BMAX;
  } 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]->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 _validate_integer {
  my($n) = @_;
  croak "Parameter must be defined" if !defined $n;
  if (ref($n) eq 'CODE') {
    $_[0] = $_[0]->();
    $n = $_[0];
  }
  my $poscmp = OLD_PERL_VERSION ?  562949953421312 : ''.~0;
  my $negcmp = OLD_PERL_VERSION ? -562949953421312 : -(~0 >> 1);
  if (ref($n) eq 'Math::BigInt') {
    croak "Parameter '$n' must be an integer" if !$n->is_int();
    $_[0] = _bigint_to_int($_[0]) if $n <= $poscmp && $n >= $negcmp;
  } else {
    my $strn = "$n";
    croak "Parameter '$strn' must be an integer"
      if $strn eq '' || ($strn =~ tr/-0123456789//c && $strn !~ /^[-+]?\d+$/);
    if ($n <= $poscmp && $n >= $negcmp) {
      $_[0] = $strn if ref($n);
    } else {
      $_[0] = Math::BigInt->new($strn)
    }
  }
  $_[0]->upgrade(undef) if ref($_[0]) && $_[0]->upgrade();
  1;
}

sub _binary_search {
  my($n, $lo, $hi, $sub, $exitsub) = @_;
  while ($lo < $hi) {
    my $mid = $lo + int(($hi-$lo) >> 1);
    return $mid if defined $exitsub && $exitsub->($n,$lo,$hi);
    if ($sub->($mid) < $n) { $lo = $mid+1; }
    else                   { $hi = $mid;   }
  }
  return $lo-1;
}

my @_primes_small = (0,2);
{
  my($n, $s, $sieveref) = (7-2, 3, _sieve_erat_string(5003));
  push @_primes_small, 2*pos($$sieveref)-1 while $$sieveref =~ m/0/g;
}
my @_prime_next_small = (
   2,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23,
   29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41,41,43,43,47,
   47,47,47,53,53,53,53,53,53,59,59,59,59,59,59,61,61,67,67,67,67,67,67,71);

# For wheel-30
my @_prime_indices = (1, 7, 11, 13, 17, 19, 23, 29);
my @_nextwheel30 = (1,7,7,7,7,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23,29,29,29,29,29,29,1);
my @_prevwheel30 = (29,29,1,1,1,1,1,1,7,7,7,7,11,11,13,13,13,13,17,17,19,19,19,19,23,23,23,23,23,23);
my @_wheeladvance30 = (1,6,5,4,3,2,1,4,3,2,1,2,1,4,3,2,1,2,1,4,3,2,1,6,5,4,3,2,1,2);
my @_wheelretreat30 = (1,2,1,2,3,4,5,6,1,2,3,4,1,2,1,2,3,4,1,2,1,2,3,4,1,2,3,4,5,6);

sub _tiny_prime_count {
  my($n) = @_;
  return if $n >= $_primes_small[-1];
  my $j = $#_primes_small;
  my $i = 1 + ($n >> 4);
  while ($i < $j) {
    my $mid = ($i+$j)>>1;
    if ($_primes_small[$mid] <= $n) { $i = $mid+1; }
    else                            { $j = $mid;   }
  }
  return $i-1;
}

sub _is_prime7 {  # n must not be divisible by 2, 3, or 5
  my($n) = @_;

  $n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt' && $n <= BMAX;
  if (ref($n) eq 'Math::BigInt') {
    return 0 unless Math::BigInt::bgcd($n, B_PRIM767)->is_one;
    return 0 unless _miller_rabin_2($n);
    my $is_esl_prime = is_extra_strong_lucas_pseudoprime($n);
    return ($is_esl_prime)  ?  (($n <= "18446744073709551615") ? 2 : 1)  :  0;
  }

  if ($n < 61*61) {
    foreach my $i (qw/7 11 13 17 19 23 29 31 37 41 43 47 53 59/) {
      return 2 if $i*$i > $n;
      return 0 if !($n % $i);
    }
    return 2;
  }

  return 0 if !($n %  7) || !($n % 11) || !($n % 13) || !($n % 17) ||
              !($n % 19) || !($n % 23) || !($n % 29) || !($n % 31) ||
              !($n % 37) || !($n % 41) || !($n % 43) || !($n % 47) ||
              !($n % 53) || !($n % 59);

  # We could do:
  #   return is_strong_pseudoprime($n, (2,299417)) if $n < 19471033;
  # or:
  #   foreach my $p (@_primes_small[18..168]) {
  #     last if $p > $limit;
  #     return 0 unless $n % $p;
  #   }
  #   return 2;

  if ($n <= 1_500_000) {
    my $limit = int(sqrt($n));
    my $i = 61;
    while (($i+30) <= $limit) {
      return 0 unless ($n% $i    ) && ($n%($i+ 6)) &&
                      ($n%($i+10)) && ($n%($i+12)) &&
                      ($n%($i+16)) && ($n%($i+18)) &&
                      ($n%($i+22)) && ($n%($i+28));
      $i += 30;
    }
    for my $inc (6,4,2,4,2,4,6,2) {
      last if $i > $limit;
      return 0 if !($n % $i);
      $i += $inc;
    }
    return 2;
  }

  if ($n < 47636622961201) {   # BPSW seems to be faster after this
    # Deterministic set of Miller-Rabin tests.  If the MR routines can handle
    # bases greater than n, then this can be simplified.
    my @bases;
    # n > 1_000_000 because of the previous block.
    if    ($n <         19471033) { @bases = ( 2,  299417); }
    elsif ($n <         38010307) { @bases = ( 2,  9332593); }
    elsif ($n <        316349281) { @bases = ( 11000544, 31481107); }
    elsif ($n <       4759123141) { @bases = ( 2, 7, 61); }
    elsif ($n <     154639673381) { @bases = ( 15, 176006322, 4221622697); }
    elsif ($n <   47636622961201) { @bases = ( 2, 2570940, 211991001, 3749873356); }
    elsif ($n < 3770579582154547) { @bases = ( 2, 2570940, 880937, 610386380, 4130785767); }
    else                          { @bases = ( 2, 325, 9375, 28178, 450775, 9780504, 1795265022); }
    return is_strong_pseudoprime($n, @bases)  ?  2  :  0;
  }

  # Inlined BPSW
  return 0 unless _miller_rabin_2($n);
  return is_almost_extra_strong_lucas_pseudoprime($n) ? 2 : 0;
}

sub is_prime {
  my($n) = @_;
  return 0 if defined($n) && int($n) < 0;
  _validate_positive_integer($n);

  if (ref($n) eq 'Math::BigInt') {
    return 0 unless Math::BigInt::bgcd($n, B_PRIM235)->is_one;
  } else {
    if ($n < 7) { return ($n == 2) || ($n == 3) || ($n == 5) ? 2 : 0; }
    return 0 if !($n % 2) || !($n % 3) || !($n % 5);
  }
  return _is_prime7($n);
}

# is_prob_prime is the same thing for us.
*is_prob_prime = \&is_prime;

# BPSW probable prime.  No composites are known to have passed this test
# since it was published in 1980, though we know infinitely many exist.
# It has also been verified that no 64-bit composite will return true.
# Slow since it's all in PP and uses bigints.
sub is_bpsw_prime {
  my($n) = @_;
  return 0 if defined($n) && int($n) < 0;
  _validate_positive_integer($n);
  return 0 unless _miller_rabin_2($n);
  if ($n <= 18446744073709551615) {
    return is_almost_extra_strong_lucas_pseudoprime($n) ? 2 : 0;
  }
  return is_extra_strong_lucas_pseudoprime($n) ? 1 : 0;
}

sub is_provable_prime {
  my($n) = @_;
  return 0 if defined $n && $n < 2;
  _validate_positive_integer($n);
  if ($n <= 18446744073709551615) {
    return 0 unless _miller_rabin_2($n);
    return 0 unless is_almost_extra_strong_lucas_pseudoprime($n);
    return 2;
  }
  my($is_prime, $cert) = Math::Prime::Util::is_provable_prime_with_cert($n);
  $is_prime;
}

# Possible sieve storage:
#   1) vec with mod-30 wheel:   8 bits  / 30
#   2) vec with mod-2 wheel :  15 bits  / 30
#   3) str with mod-30 wheel:   8 bytes / 30
#   4) str with mod-2 wheel :  15 bytes / 30
#
# It looks like using vecs is about 2x slower than strs, and the strings also
# let us do some fast operations on the results.  E.g.
#   Count all primes:
#      $count += $$sieveref =~ tr/0//;
#   Loop over primes:
#      foreach my $s (split("0", $$sieveref, -1)) {
#        $n += 2 + 2 * length($s);
#        .. do something with the prime $n
#      }
#
# We're using method 4, though sadly it is memory intensive relative to the
# other methods.  I will point out that it is 30-60x less memory than sieves
# using an array, and the performance of this function is over 10x that
# of naive sieves.

sub _sieve_erat_string {
  my($end) = @_;
  $end-- if ($end & 1) == 0;
  my $s_end = $end >> 1;

  my $whole = int( $s_end / 15);   # Prefill with 3 and 5 already marked.
  croak "Sieve too large" if $whole > 1_145_324_612;  # ~32 GB string
  my $sieve = '100010010010110' . '011010010010110' x $whole;
  substr($sieve, $s_end+1) = '';   # Ensure we don't make too many entries
  my ($n, $limit) = ( 7, int(sqrt($end)) );
  while ( $n <= $limit ) {
    for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) {
      substr($sieve, $s, 1) = '1';
    }
    do { $n += 2 } while substr($sieve, $n>>1, 1);
  }
  return \$sieve;
}

# TODO: this should be plugged into precalc, memfree, etc. just like the C code
{
  my $primary_size_limit = 15000;
  my $primary_sieve_size = 0;
  my $primary_sieve_ref;
  sub _sieve_erat {
    my($end) = @_;

    return _sieve_erat_string($end) if $end > $primary_size_limit;

    if ($primary_sieve_size == 0) {
      $primary_sieve_size = $primary_size_limit;
      $primary_sieve_ref = _sieve_erat_string($primary_sieve_size);
    }
    my $sieve = substr($$primary_sieve_ref, 0, ($end+1)>>1);
    return \$sieve;
  }
}


sub _sieve_segment {
  my($beg,$end,$limit) = @_;
  ($beg, $end) = map { _bigint_to_int($_) } ($beg, $end)
    if ref($end) && $end <= BMAX;
  croak "Internal error: segment beg is even" if ($beg % 2) == 0;
  croak "Internal error: segment end is even" if ($end % 2) == 0;
  croak "Internal error: segment end < beg" if $end < $beg;
  croak "Internal error: segment beg should be >= 3" if $beg < 3;
  my $range = int( ($end - $beg) / 2 ) + 1;

  # Prefill with 3 and 5 already marked, and offset to the segment start.
  my $whole = int( ($range+14) / 15);
  my $startp = ($beg % 30) >> 1;
  my $sieve = substr('011010010010110', $startp) . '011010010010110' x $whole;
  # Set 3 and 5 to prime if we're sieving them.
  substr($sieve,0,2) = '00' if $beg == 3;
  substr($sieve,0,1) = '0'  if $beg == 5;
  # Get rid of any extra we added.
  substr($sieve, $range) = '';

  # If the end value is below 7^2, then the pre-sieve is all we needed.
  return \$sieve if $end < 49;

  my $sqlimit = ref($end) ? $end->copy->bsqrt() : int(sqrt($end)+0.0000001);
  $limit = $sqlimit if !defined $limit || $sqlimit < $limit;
  # For large value of end, it's a huge win to just walk primes.

  my($p, $s, $primesieveref) = (7-2, 3, _sieve_erat($limit));
  while ( (my $nexts = 1 + index($$primesieveref, '0', $s)) > 0 ) {
    $p += 2 * ($nexts - $s);
    $s = $nexts;
    my $p2 = $p*$p;
    if ($p2 < $beg) {
      my $f = 1+int(($beg-1)/$p);
      $f++ unless $f % 2;
      $p2 = $p * $f;
    }
    # With large bases and small segments, it's common to find we don't hit
    # the segment at all.  Skip all the setup if we find this now.
    if ($p2 <= $end) {
      # Inner loop marking multiples of p
      # (everything is divided by 2 to keep inner loop simpler)
      my $filter_end = ($end - $beg) >> 1;
      my $filter_p2  = ($p2  - $beg) >> 1;
      while ($filter_p2 <= $filter_end) {
        substr($sieve, $filter_p2, 1) = "1";
        $filter_p2 += $p;
      }
    }
  }
  \$sieve;
}

sub trial_primes {
  my($low,$high) = @_;
  if (!defined $high) {
    $high = $low;
    $low = 2;
  }
  _validate_positive_integer($low);
  _validate_positive_integer($high);
  return if $low > $high;
  my @primes;

  # For a tiny range, just use next_prime calls
  if (($high-$low) < 1000) {
    $low-- if $low >= 2;
    my $curprime = next_prime($low);
    while ($curprime <= $high) {
      push @primes, $curprime;
      $curprime = next_prime($curprime);
    }
    return \@primes;
  }

  # Sieve to 10k then BPSW test
  push @primes, 2  if ($low <= 2) && ($high >= 2);
  push @primes, 3  if ($low <= 3) && ($high >= 3);
  push @primes, 5  if ($low <= 5) && ($high >= 5);
  $low = 7 if $low < 7;
  $low++ if ($low % 2) == 0;
  $high-- if ($high % 2) == 0;
  my $sieveref = _sieve_segment($low, $high, 10000);
  my $n = $low-2;
  while ($$sieveref =~ m/0/g) {
    my $p = $n+2*pos($$sieveref);
    push @primes, $p if _miller_rabin_2($p) && is_extra_strong_lucas_pseudoprime($p);
  }
  return \@primes;
}

sub primes {
  my($low,$high) = @_;
  if (scalar @_ > 1) {
    _validate_positive_integer($low);
    _validate_positive_integer($high);
    $low = 2 if $low < 2;
  } else {
    ($low,$high) = (2, $low);
    _validate_positive_integer($high);
  }
  my $sref = [];
  return $sref if ($low > $high) || ($high < 2);
  return [grep { $_ >= $low && $_ <= $high } @_primes_small]
    if $high <= $_primes_small[-1];

  return [ Math::Prime::Util::GMP::sieve_primes($low, $high, 0) ]
    if $Math::Prime::Util::_GMPfunc{"sieve_primes"} && $Math::Prime::Util::GMP::VERSION >= 0.34;

  # At some point even the pretty-fast pure perl sieve is going to be a
  # dog, and we should move to trials.  This is typical with a small range
  # on a large base.  More thought on the switchover should be done.
  return trial_primes($low, $high) if ref($low)  eq 'Math::BigInt'
                                   || ref($high) eq 'Math::BigInt'
                                   || ($low > 1_000_000_000_000 && ($high-$low) < int($low/1_000_000));

  push @$sref, 2  if ($low <= 2) && ($high >= 2);
  push @$sref, 3  if ($low <= 3) && ($high >= 3);
  push @$sref, 5  if ($low <= 5) && ($high >= 5);
  $low = 7 if $low < 7;
  $low++ if ($low % 2) == 0;
  $high-- if ($high % 2) == 0;
  return $sref if $low > $high;

  my($n,$sieveref);
  if ($low == 7) {
    $n = 0;
    $sieveref = _sieve_erat($high);
    substr($$sieveref,0,3,'111');
  } else {
    $n = $low-1;
    $sieveref = _sieve_segment($low,$high);
  }
  push @$sref, $n+2*pos($$sieveref)-1 while $$sieveref =~ m/0/g;
  $sref;
}

sub sieve_range {
  my($n, $width, $depth) = @_;
  _validate_positive_integer($n);
  _validate_positive_integer($width);
  _validate_positive_integer($depth);

  my @candidates;
  my $start = $n;

  if ($n < 5) {
    push @candidates, (2-$n) if $n <= 2 && $n+$width-1 >= 2;
    push @candidates, (3-$n) if $n <= 3 && $n+$width-1 >= 3;
    push @candidates, (4-$n) if $n <= 4 && $n+$width-1 >= 4 && $depth < 2;
    $start = 5;
    $width -= ($start - $n);
  }

  return @candidates, map {$start+$_-$n } 0 .. $width-1 if $depth < 2;
  return @candidates, map { $_ - $n }
                      grep { ($_ & 1) && ($depth < 3 || ($_ % 3)) }
                      map { $start+$_ }
                      0 .. $width-1                     if $depth < 5;

  if (!($start & 1)) { $start++; $width--; }
  $width-- if !($width&1);
  return @candidates if $width < 1;

  my $sieveref = _sieve_segment($start, $start+$width-1, $depth);
  my $offset = $start - $n - 2;
  while ($$sieveref =~ m/0/g) {
    push @candidates, $offset + (pos($$sieveref) << 1);
  }
  return @candidates;
}

sub sieve_prime_cluster {
  my($lo,$hi,@cl) = @_;
  my $_verbose = Math::Prime::Util::prime_get_config()->{'verbose'};
  _validate_positive_integer($lo);
  _validate_positive_integer($hi);

  if ($Math::Prime::Util::_GMPfunc{"sieve_prime_cluster"}) {
    return map { ($_ > ''.~0) ? Math::BigInt->new(''.$_) : $_ }
           Math::Prime::Util::GMP::sieve_prime_cluster($lo,$hi,@cl);
  }

  return @{primes($lo,$hi)} if scalar(@cl) == 0;

  unshift @cl, 0;
  for my $i (1 .. $#cl) {
    _validate_positive_integer($cl[$i]);
    croak "sieve_prime_cluster: values must be even" if $cl[$i] & 1;
    croak "sieve_prime_cluster: values must be increasing" if $cl[$i] <= $cl[$i-1];
  }
  my($p,$sievelim,@p) = (17, 2000);
  $p = 13 if ($hi-$lo) < 50_000_000;
  $p = 11 if ($hi-$lo) <  1_000_000;
  $p =  7 if ($hi-$lo) <     20_000 && $lo < INTMAX;

  # Add any cases under our sieving point.
  if ($lo <= $sievelim) {
    $sievelim = $hi if $sievelim > $hi;
    for my $n (@{primes($lo,$sievelim)}) {
      my $ac = 1;
      for my $ci (1 .. $#cl) {
        if (!is_prime($n+$cl[$ci])) { $ac = 0; last; }
      }
      push @p, $n if $ac;
    }
    $lo = next_prime($sievelim);
  }
  return @p if $lo > $hi;

  # Compute acceptable residues.
  my $pr = primorial($p);
  my $startpr = _bigint_to_int($lo % $pr);

  my @acc = grep { ($_ & 1) && $_%3 }  ($startpr .. $startpr + $pr - 1);
  for my $c (@cl) {
    if ($p >= 7) {
      @acc = grep { (($_+$c)%3) && (($_+$c)%5) && (($_+$c)%7) } @acc;
    } else {
      @acc = grep { (($_+$c)%3)  && (($_+$c)%5) } @acc;
    }
  }
  for my $c (@cl) {
    @acc = grep { Math::Prime::Util::gcd($_+$c,$pr) == 1 } @acc;
  }
  @acc = map { $_-$startpr } @acc;

  print "cluster sieve using ",scalar(@acc)," residues mod $pr\n" if $_verbose;
  return @p if scalar(@acc) == 0;

  # Prepare table for more sieving.
  my @mprimes = @{primes( $p+1, $sievelim)};
  my @vprem;
  for my $p (@mprimes) {
    for my $c (@cl) {
      $vprem[$p]->[ ($p-($c%$p)) % $p ] = 1;
    }
  }

  # Walk the range in primorial chunks, doing primality tests.
  my $nprim = 0;
  while ($lo <= $hi) {

    my @racc = @acc;

    # Make sure we don't do anything past the limit
    if (($lo+$acc[-1]) > $hi) {
      my $max = _bigint_to_int($hi-$lo);
      @racc = grep { $_ <= $max } @racc;
    }

    # Sieve more values using native math
    foreach my $p (@mprimes) {
      my $rem = _bigint_to_int( $lo % $p );
      @racc = grep { !$vprem[$p]->[ ($rem+$_) % $p ] } @racc;
      last unless scalar(@racc);
    }

    # Do final primality tests.
    for my $c (@cl) {
      last unless scalar(@racc);
      my $loc = $lo + $c;
      $nprim += scalar(@racc);
      @racc = grep { Math::Prime::Util::is_prime($loc+$_) } @racc;
    }

    push @p, map { $lo + $_ } @racc;
    $lo += $pr;
  }
  print "cluster sieve ran $nprim primality tests\n" if $_verbose;
  @p;
}


sub _n_ramanujan_primes {
  my($n) = @_;
  return [] if $n <= 0;
  my $max = nth_prime_upper(int(48/19*$n)+1);
  my @L = (2, (0) x $n-1);
  my $s = 1;
  for (my $k = 7; $k <= $max; $k += 2) {
    $s++ if is_prime($k);
    $L[$s] = $k+1 if $s < $n;
    $s-- if ($k&3) == 1 && is_prime(($k+1)>>1);
    $L[$s] = $k+2 if $s < $n;
  }
  \@L;
}

sub _ramanujan_primes {
  my($low,$high) = @_;
  ($low,$high) = (2, $low) unless defined $high;
  return [] if ($low > $high) || ($high < 2);
  my $nn = prime_count_upper($high) >> 1;
  my $L = _n_ramanujan_primes($nn);
  shift @$L while @$L && $L->[0] < $low;
  pop @$L while @$L && $L->[-1] > $high;
  $L;
}

sub is_ramanujan_prime {
  my($n) = @_;
  return 1 if $n == 2;
  return 0 if $n < 11;
  my $L = _ramanujan_primes($n,$n);
  return (scalar(@$L) > 0) ? 1 : 0;
}

sub nth_ramanujan_prime {
  my($n) = @_;
  return undef if $n <= 0;  ## no critic qw(ProhibitExplicitReturnUndef)
  my $L = _n_ramanujan_primes($n);
  return $L->[$n-1];
}

sub next_prime {
  my($n) = @_;
  _validate_positive_integer($n);
  return $_prime_next_small[$n] if $n <= $#_prime_next_small;
  # This turns out not to be faster.
  # return $_primes_small[1+_tiny_prime_count($n)] if $n < $_primes_small[-1];

  return Math::BigInt->new(MPU_32BIT ? "4294967311" : "18446744073709551629")
    if ref($n) ne 'Math::BigInt' && $n >= MPU_MAXPRIME;
  # n is now either 1) not bigint and < maxprime, or (2) bigint and >= uvmax

  if ($n > 4294967295 && Math::Prime::Util::prime_get_config()->{'gmp'}) {
    return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::next_prime($n));
  }

  if (ref($n) eq 'Math::BigInt') {
    do {
      $n += $_wheeladvance30[$n%30];
    } while !Math::BigInt::bgcd($n, B_PRIM767)->is_one ||
            !_miller_rabin_2($n) || !is_extra_strong_lucas_pseudoprime($n);
  } else {
    do {
      $n += $_wheeladvance30[$n%30];
    } while !($n%7) || !_is_prime7($n);
  }
  $n;
}

sub prev_prime {
  my($n) = @_;
  _validate_positive_integer($n);
  return (undef,undef,undef,2,3,3,5,5,7,7,7,7)[$n] if $n <= 11;
  if ($n > 4294967295 && Math::Prime::Util::prime_get_config()->{'gmp'}) {
    return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::prev_prime($n));
  }

  if (ref($n) eq 'Math::BigInt') {
    do {
      $n -= $_wheelretreat30[$n%30];
    } while !Math::BigInt::bgcd($n, B_PRIM767)->is_one ||
            !_miller_rabin_2($n) || !is_extra_strong_lucas_pseudoprime($n);
  } else {
    do {
      $n -= $_wheelretreat30[$n%30];
    } while !($n%7) || !_is_prime7($n);
  }
  $n;
}

sub partitions {
  my $n = shift;

  my $d = int(sqrt($n+1));
  my @pent = (1, map { (($_*(3*$_+1))>>1, (($_+1)*(3*$_+2))>>1) } 1 .. $d);
  my $ZERO = ($n >= ((~0 > 4294967295) ? 400 : 270)) ? BZERO : 0;
  my @part = ($ZERO+1);
  foreach my $j (scalar @part .. $n) {
    my ($psum1, $psum2, $k) = ($ZERO, $ZERO, 1);
    foreach my $p (@pent) {
      last if $p > $j;
      if ((++$k) & 2) { $psum1 += $part[ $j - $p ] }
      else            { $psum2 += $part[ $j - $p ] }
    }
    $part[$j] = $psum1 - $psum2;
  }
  return $part[$n];
}

sub primorial {
  my $n = shift;

  my @plist = @{primes($n)};
  my $max = (MPU_32BIT) ? 29 : (OLD_PERL_VERSION) ? 43 : 53;

  # If small enough, multiply the small primes.
  if ($n < $max) {
    my $pn = 1;
    $pn *= $_ for @plist;
    return $pn;
  }

  # Otherwise, combine them as UVs, then combine using product tree.
  my $i = 0;
  while ($i < $#plist) {
    my $m = $plist[$i] * $plist[$i+1];
    if ($m <= INTMAX) { splice(@plist, $i, 2, $m); }
    else              { $i++;                      }
  }
  vecprod(@plist);
}

sub consecutive_integer_lcm {
  my $n = shift;

  my $max = (MPU_32BIT) ? 22 : (OLD_PERL_VERSION) ? 37 : 46;
  my $pn = ref($n) ? ref($n)->new(1) : ($n >= $max) ? Math::BigInt->bone() : 1;
  for (my $p = 2; $p <= $n; $p = next_prime($p)) {
    my($p_power, $pmin) = ($p, int($n/$p));
    $p_power *= $p while $p_power <= $pmin;
    $pn *= $p_power;
  }
  $pn = _bigint_to_int($pn) if $pn <= BMAX;
  return $pn;
}

sub jordan_totient {
  my($k, $n) = @_;
  return ($n == 1) ? 1 : 0  if $k == 0;
  return euler_phi($n)      if $k == 1;
  return ($n == 1) ? 1 : 0  if $n <= 1;

  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::jordan_totient($k, $n))
    if $Math::Prime::Util::_GMPfunc{"jordan_totient"};


  my @pe = Math::Prime::Util::factor_exp($n);
  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
  my $totient = BONE->copy;
  foreach my $f (@pe) {
    my ($p, $e) = @$f;
    $p = Math::BigInt->new("$p")->bpow($k);
    $totient->bmul($p->copy->bdec());
    $totient->bmul($p) for 2 .. $e;
  }
  $totient = _bigint_to_int($totient) if $totient->bacmp(BMAX) <= 0;
  return $totient;
}

sub euler_phi {
  return euler_phi_range(@_) if scalar @_ > 1;
  my($n) = @_;

  return Math::Prime::Util::_reftyped($_[0],Math::Prime::Util::GMP::totient($n))
    if $Math::Prime::Util::_GMPfunc{"totient"};

  _validate_positive_integer($n);
  return 0 if $n < 0;
  return $n if $n <= 1;

  my $totient = $n - $n + 1;

  # Fast reduction of multiples of 2, may also reduce n for factoring
  if (ref($n) eq 'Math::BigInt') {
    my $s = 0;
    if ($n->is_even) {
      do { $n->brsft(BONE); $s++; } while $n->is_even;
      $totient->blsft($s-1) if $s > 1;
    }
  } else {
    while (($n % 4) == 0) { $n >>= 1;  $totient <<= 1; }
    if (($n % 2) == 0) { $n >>= 1; }
  }

  my @pe = Math::Prime::Util::factor_exp($n);

  if (ref($n) ne 'Math::BigInt') {
    foreach my $f (@pe) {
      my ($p, $e) = @$f;
      $totient *= $p - 1;
      $totient *= $p for 2 .. $e;
    }
  } else {
    my $zero = $n->copy->bzero;
    foreach my $f (@pe) {
      my ($p, $e) = @$f;
      $p = $zero->copy->badd("$p");
      $totient->bmul($p->copy->bdec());
      $totient->bmul($p) for 2 .. $e;
    }
  }
  $totient = _bigint_to_int($totient) if ref($totient) eq 'Math::BigInt'
                                      && $totient->bacmp(BMAX) <= 0;
  return $totient;
}

sub euler_phi_range {
  my($n, $nend) = @_;
  return () if $nend < $n;
  return euler_phi($n) if $n == $nend;
  my @totients;
  if ($nend > 2**30) {
    while ($n < $nend) {
      push @totients, euler_phi($n++);
    }
  } else {
    @totients = (0 .. $nend);
    foreach my $i (2 .. $nend) {
      next unless $totients[$i] == $i;
      $totients[$i] = $i-1;
      foreach my $j (2 .. int($nend / $i)) {
        $totients[$i*$j] -= $totients[$i*$j]/$i;
      }
    }
    splice(@totients, 0, $n) if $n > 0;
  }
  return @totients;
}

sub moebius {
  return moebius_range(@_) if scalar @_ > 1;
  my($n) = @_;
  return ($n == 1) ? 1 : 0  if $n <= 1;
  return 0 if ($n >= 49) && (!($n % 4) || !($n % 9) || !($n % 25) || !($n%49) );
  my @factors = Math::Prime::Util::factor($n);
  foreach my $i (1 .. $#factors) {
    return 0 if $factors[$i] == $factors[$i-1];
  }
  return ((scalar @factors) % 2) ? -1 : 1;
}
sub is_square_free {
  return (Math::Prime::Util::moebius($_[0]) != 0) ? 1 : 0;
}
sub is_semiprime {
  my($n) = @_;
  _validate_positive_integer($n);
  return ($n == 4) if $n < 6;
  return (Math::Prime::Util::is_prob_prime($n>>1) ? 1 : 0) if ($n % 2) == 0;
  return (Math::Prime::Util::is_prob_prime($n/3)  ? 1 : 0) if ($n % 3) == 0;
  return (Math::Prime::Util::is_prob_prime($n/5)  ? 1 : 0) if ($n % 5) == 0;
  {
    my @f = trial_factor($n, 4999);
    return 0 if @f > 2;
    return (_is_prime7($f[1]) ? 1 : 0) if @f == 2;
  }
  return 0 if _is_prime7($n);
  {
    my @f = pminus1_factor ($n, 250_000);
    return 0 if @f > 2;
    return (_is_prime7($f[1]) ? 1 : 0) if @f == 2;
  }
  {
    my @f = pbrent_factor ($n, 128*1024, 3, 1);
    return 0 if @f > 2;
    return (_is_prime7($f[1]) ? 1 : 0) if @f == 2;
  }
  return (scalar(Math::Prime::Util::factor($n)) == 2) ? 1 : 0;
}


sub moebius_range {
  my($lo, $hi) = @_;
  return () if $hi < $lo;
  return moebius($lo) if $lo == $hi;
  if ($hi > 2**32) {
    my @mu;
    while ($lo <= $hi) {
      push @mu, moebius($lo++);
    }
    return @mu;
  }
  my @mu = map { 1 } $lo .. $hi;
  $mu[0] = 0 if $lo == 0;
  my($p, $sqrtn) = (2, int(sqrt($hi)+0.5));
  while ($p <= $sqrtn) {
    my $i = $p * $p;
    $i = $i * int($lo/$i) + (($lo % $i)  ? $i : 0)  if $i < $lo;
    while ($i <= $hi) {
      $mu[$i-$lo] = 0;
      $i += $p * $p;
    }
    $i = $p;
    $i = $i * int($lo/$i) + (($lo % $i)  ? $i : 0)  if $i < $lo;
    while ($i <= $hi) {
      $mu[$i-$lo] *= -$p;
      $i += $p;
    }
    $p = next_prime($p);
  }
  foreach my $i ($lo .. $hi) {
    my $m = $mu[$i-$lo];
    $m *= -1 if abs($m) != $i;
    $mu[$i-$lo] = ($m>0) - ($m<0);
  }
  return @mu;
}

sub mertens {
  my($n) = @_;
    # This is the most basic Deléglise and Rivat algorithm.  u = n^1/2
  # and no segmenting is done.  Their algorithm uses u = n^1/3, breaks
  # the summation into two parts, and calculates those in segments.  Their
  # computation time growth is half of this code.
  return $n if $n <= 1;
  my $u = int(sqrt($n));
  my @mu = (0, Math::Prime::Util::moebius(1, $u)); # Hold values of mu for 0-u
  my $musum = 0;
  my @M = map { $musum += $_; } @mu;     # Hold values of M for 0-u
  my $sum = $M[$u];
  foreach my $m (1 .. $u) {
    next if $mu[$m] == 0;
    my $inner_sum = 0;
    my $lower = int($u/$m) + 1;
    my $last_nmk = int($n/($m*$lower));
    my ($denom, $this_k, $next_k) = ($m, 0, int($n/($m*1)));
    for my $nmk (1 .. $last_nmk) {
      $denom += $m;
      $this_k = int($n/$denom);
      next if $this_k == $next_k;
      ($this_k, $next_k) = ($next_k, $this_k);
      $inner_sum += $M[$nmk] * ($this_k - $next_k);
    }
    $sum -= $mu[$m] * $inner_sum;
  }
  return $sum;
}

sub ramanujan_sum {
  my($k,$n) = @_;
  return 0 if $k < 1 || $n <  1;
  my $g = $k / Math::Prime::Util::gcd($k,$n);
  my $m = Math::Prime::Util::moebius($g);
  return $m if $m == 0 || $k == $g;
  $m * (Math::Prime::Util::euler_phi($k) / Math::Prime::Util::euler_phi($g));
}

sub liouville {
  my($n) = @_;
  my $l = (-1) ** scalar Math::Prime::Util::factor($n);
  return $l;
}

# Exponential of Mangoldt function (A014963).
# Return p if n = p^m [p prime, m >= 1], 1 otherwise.
sub exp_mangoldt {
  my($n) = @_;
  my $p;
  return 1 unless Math::Prime::Util::is_prime_power($n,\$p);
  $p;
}

sub carmichael_lambda {
  my($n) = @_;
  return euler_phi($n) if $n < 8;                # = phi(n) for n < 8
  return euler_phi($n)/2 if ($n & ($n-1)) == 0;  # = phi(n)/2 for 2^k, k>2

  my @pe = Math::Prime::Util::factor_exp($n);
  $pe[0]->[1]-- if $pe[0]->[0] == 2 && $pe[0]->[1] > 2;

  my $lcm = Math::BigInt::blcm(
    map { $_->[0]->copy->bpow($_->[1]->copy->bdec)->bmul($_->[0]->copy->bdec) }
    map { [ map { Math::BigInt->new("$_") } @$_ ] }
    @pe
  );
  $lcm = _bigint_to_int($lcm) if $lcm->bacmp(BMAX) <= 0;
  return $lcm;
}

sub is_carmichael {
  my($n) = @_;
  _validate_positive_integer($n);

  # This works fine, but very slow
  # return !is_prime($n) && ($n % carmichael_lambda($n)) == 1;

  return 0 if $n < 561 || ($n % 2) == 0;
  return 0 if (!($n % 4) || !($n % 9) || !($n % 25) || !($n%49) || !($n%121));

  my $fn = $n;
  if ($n > 100_000_000) {   # After 100M, this saves time on average
    # Pre-tests which are faster than factoring.
    return 0 if Math::Prime::Util::powmod(2, $n-1, $n) != 1;
    return 0 if Math::Prime::Util::is_prime($n);
    for my $a (3,5,7,11,13,17,19,23,29,31,37) {
      my $gcd = Math::Prime::Util::gcd($a, $fn);
      if ($gcd == 1) {
        return 0 if Math::Prime::Util::powmod($a, $n-1, $n) != 1;
      } else {
        return 0 if $gcd != $a;              # Not square free
        return 0 if (($n-1) % ($a-1)) != 0;  # factor doesn't divide
        $fn /= $a;
      }
    }
  }
  #return 1;
  # Based on pre-tests, it's reasonably likely $n is a Carmichael number.

  # Use probabilistic test if too large to reasonably factor.
  if (length($fn) > 50) {
    for my $t (13 .. 150) {
      my $a = $_primes_small[$t];
      my $gcd = Math::Prime::Util::gcd($a, $fn);
      if ($gcd == 1) {
        return 0 if Math::Prime::Util::powmod($a, $n-1, $n) != 1;
      } else {
        return 0 if $gcd != $a;              # Not square free
        return 0 if (($n-1) % ($a-1)) != 0;  # factor doesn't divide
        $fn /= $a;
      }
    }
    return 1;
  }

  # Verify with factoring.
  my @pe = Math::Prime::Util::factor_exp($n);
  return 0 if scalar(@pe) < 3;
  for my $pe (@pe) {
    return 0 if $pe->[1] > 1 || (($n-1) % ($pe->[0]-1)) != 0;
  }
  1;
}

sub is_quasi_carmichael {
  my($n) = @_;
  _validate_positive_integer($n);

  return 0 if $n < 35;
  return 0 if (!($n % 4) || !($n % 9) || !($n % 25) || !($n%49) || !($n%121));

  my @pe = Math::Prime::Util::factor_exp($n);
  # Not quasi-Carmichael if prime
  return 0 if scalar(@pe) < 2;
  # Not quasi-Carmichael if not square free
  for my $pe (@pe) {
    return 0 if $pe->[1] > 1;
  }
  my @f = map { $_->[0] } @pe;
  my $nbases = 0;
  if ($n < 2000) {
    # In theory for performance, but mainly keeping to show direct method.
    my $lim = $f[-1];
    $lim = (($n-$lim*$lim) + $lim - 1) / $lim;
    for my $b (1 .. $f[0]-1) {
      my $nb = $n - $b;
      $nbases++ if Math::Prime::Util::vecall(sub { $nb % ($_-$b) == 0 }, @f);
    }
    if (scalar(@f) > 2) {
      for my $b (1 .. $lim-1) {
        my $nb = $n + $b;
        $nbases++ if Math::Prime::Util::vecall(sub { $nb % ($_+$b) == 0 }, @f);
      }
    }
  } else {
    my($spf,$lpf) = ($f[0], $f[-1]);
    if (scalar(@f) == 2) {
      foreach my $d (Math::Prime::Util::divisors($n/$spf - 1)) {
        my $k = $spf - $d;
        my $p = $n - $k;
        last if $d >= $spf;
        $nbases++ if Math::Prime::Util::vecall(sub { my $j = $_-$k;  $j && ($p % $j) == 0 }, @f);
      }
    } else {
      foreach my $d (Math::Prime::Util::divisors($lpf * ($n/$lpf - 1))) {
        my $k = $lpf - $d;
        my $p = $n - $k;
        next if $k == 0 || $k >= $spf;
        $nbases++ if Math::Prime::Util::vecall(sub { my $j = $_-$k;  $j && ($p % $j) == 0 }, @f);
      }
    }
  }
  $nbases;
}

sub is_pillai {
  my($p) = @_;
  return 0 if defined($p) && int($p) < 0;
  _validate_positive_integer($p);
  return 0 if $p <= 2;

  my $pm1 = $p-1;
  my $nfac = 5040 % $p;
  for (my $n = 8; $n < $p; $n++) {
    $nfac = Math::Prime::Util::mulmod($nfac, $n, $p);
    return $n if $nfac == $pm1 && ($p % $n) != 1;
  }
  0;
}

sub is_fundamental {
  my($n) = @_;
  _validate_integer($n);
  my $neg = ($n < 0);
  $n = -$n if $neg;
  my $r = $n & 15;
  if ($r) {
    my $r4 = $r & 3;
    if (!$neg) {
      return (($r ==  4) ? 0 : is_square_free($n >> 2)) if $r4 == 0;
      return is_square_free($n) if $r4 == 1;
    } else {
      return (($r == 12) ? 0 : is_square_free($n >> 2)) if $r4 == 0;
      return is_square_free($n) if $r4 == 3;
    }
  }
  0;
}

my @_ds_overflow =  # We'll use BigInt math if the input is larger than this.
  (~0 > 4294967295)
   ? (124, 3000000000000000000, 3000000000, 2487240, 64260, 7026)
   : ( 50,           845404560,      52560,    1548,   252,   84);
sub divisor_sum {
  my($n, $k) = @_;
  return ((defined $k && $k==0) ? 2 : 1) if $n == 0;
  return 1 if $n == 1;

  if (defined $k && ref($k) eq 'CODE') {
    my $sum = $n-$n;
    my $refn = ref($n);
    foreach my $d (Math::Prime::Util::divisors($n)) {
      $sum += $k->( $refn ? $refn->new("$d") : $d );
    }
    return $sum;
  }

  croak "Second argument must be a code ref or number"
    unless !defined $k || _validate_num($k) || _validate_positive_integer($k);
  $k = 1 if !defined $k;

  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::sigma($n, $k))
    if $Math::Prime::Util::_GMPfunc{"sigma"};

  my $will_overflow = ($k == 0) ? (length($n) >= $_ds_overflow[0])
                    : ($k <= 5) ? ($n >= $_ds_overflow[$k])
                    : 1;

  # The standard way is:
  #    my $pk = $f ** $k;  $product *= ($pk ** ($e+1) - 1) / ($pk - 1);
  # But we get less overflow using:
  #    my $pk = $f ** $k;  $product *= $pk**E for E in 0 .. e
  # Also separate BigInt and do fiddly bits for better performance.

  my @factors = Math::Prime::Util::factor_exp($n);
  my $product = (!$will_overflow) ? 1 : BONE->copy;
  if ($k == 0) {
    foreach my $f (@factors) {
      $product *= ($f->[1] + 1);
    }
  } elsif (!$will_overflow) {
    foreach my $f (@factors) {
      my ($p, $e) = @$f;
      my $pk = $p ** $k;
      my $fmult = $pk + 1;
      foreach my $E (2 .. $e) { $fmult += $pk**$E }
      $product *= $fmult;
    }
  } elsif (ref($n) && ref($n) ne 'Math::BigInt') {
    # This can help a lot for Math::GMP, etc.
    $product = ref($n)->new(1);
    foreach my $f (@factors) {
      my ($p, $e) = @$f;
      my $pk = ref($n)->new($p) ** $k;
      my $fmult = $pk;  $fmult++;
      if ($e >= 2) {
        my $pke = $pk;
        for (2 .. $e) { $pke *= $pk; $fmult += $pke; }
      }
      $product *= $fmult;
    }
  } else {
    my $bik = Math::BigInt->new("$k");
    foreach my $f (@factors) {
      my ($p, $e) = @$f;
      my $pk = Math::BigInt->new("$p")->bpow($bik);
      if    ($e == 1) { $pk->binc(); $product->bmul($pk); }
      elsif ($e == 2) { $pk->badd($pk*$pk)->binc(); $product->bmul($pk); }
      else {
        my $fmult = $pk->copy->binc;
        my $pke = $pk->copy;
        for my $E (2 .. $e) {
          $pke->bmul($pk);
          $fmult->badd($pke);
        }
        $product->bmul($fmult);
      }
    }
  }
  $product;
}

#############################################################################
#                       Lehmer prime count
#
#my @_s0 = (0);
#my @_s1 = (0,1);
#my @_s2 = (0,1,1,1,1,2);
my @_s3 = (0,1,1,1,1,1,1,2,2,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,7,7,7,7,8);
my @_s4 = (0,1,1,1,1,1,1,1,1,1,1,2,2,3,3,3,3,4,4,5,5,5,5,6,6,6,6,6,6,7,7,8,8,8,8,8,8,9,9,9,9,10,10,11,11,11,11,12,12,12,12,12,12,13,13,13,13,13,13,14,14,15,15,15,15,15,15,16,16,16,16,17,17,18,18,18,18,18,18,19,19,19,19,20,20,20,20,20,20,21,21,21,21,21,21,21,21,22,22,22,22,23,23,24,24,24,24,25,25,26,26,26,26,27,27,27,27,27,27,27,27,28,28,28,28,28,28,29,29,29,29,30,30,30,30,30,30,31,31,32,32,32,32,33,33,33,33,33,33,34,34,35,35,35,35,35,35,36,36,36,36,36,36,37,37,37,37,38,38,39,39,39,39,40,40,40,40,40,40,41,41,42,42,42,42,42,42,43,43,43,43,44,44,45,45,45,45,46,46,47,47,47,47,47,47,47,47,47,47,48);
sub _tablephi {
  my($x, $a) = @_;
  if ($a == 0) { return $x; }
  elsif ($a == 1) { return $x-int($x/2); }
  elsif ($a == 2) { return $x-int($x/2) - int($x/3) + int($x/6); }
  elsif ($a == 3) { return  8 * int($x /  30) + $_s3[$x %  30]; }
  elsif ($a == 4) { return 48 * int($x / 210) + $_s4[$x % 210]; }
  elsif ($a == 5) { my $xp = int($x/11);
                    return ( (48 * int($x   / 210) + $_s4[$x   % 210]) -
                             (48 * int($xp  / 210) + $_s4[$xp  % 210]) ); }
  else            { my ($xp,$x2) = (int($x/11),int($x/13));
                    my $x2p = int($x2/11);
                    return ( (48 * int($x   / 210) + $_s4[$x   % 210]) -
                             (48 * int($xp  / 210) + $_s4[$xp  % 210]) -
                             (48 * int($x2  / 210) + $_s4[$x2  % 210]) +
                             (48 * int($x2p / 210) + $_s4[$x2p % 210]) ); }
}

sub legendre_phi {
  my ($x, $a, $primes) = @_;
  return _tablephi($x,$a) if $a <= 6;
  $primes = primes(Math::Prime::Util::nth_prime_upper($a+1)) unless defined $primes;
  return ($x > 0 ? 1 : 0) if $x < $primes->[$a];

  my $sum = 0;
  my %vals = ( $x => 1 );
  while ($a > 6) {
    my $primea = $primes->[$a-1];
    my %newvals;
    while (my($v,$c) = each %vals) {
      my $sval = int($v / $primea);
      if ($sval < $primea) {
        $sum -= $c;
      } else {
        $newvals{$sval} -= $c;
      }
    }
    # merge newvals into vals
    while (my($v,$c) = each %newvals) {
      $vals{$v} += $c;
      delete $vals{$v} if $vals{$v} == 0;
    }
    $a--;
  }
  while (my($v,$c) = each %vals) {
    $sum += $c * _tablephi($v, $a);
  }
  return $sum;
}

sub _sieve_prime_count {
  my $high = shift;
  return (0,0,1,2,2,3,3)[$high] if $high < 7;
  $high-- unless ($high & 1);
  return 1 + ${_sieve_erat($high)} =~ tr/0//;
}

sub _count_with_sieve {
  my ($sref, $low, $high) = @_;
  ($low, $high) = (2, $low) if !defined $high;
  my $count = 0;
  if   ($low < 3) { $low = 3; $count++; }
  else            { $low |= 1; }
  $high-- unless ($high & 1);
  return $count if $low > $high;
  my $sbeg = $low >> 1;
  my $send = $high >> 1;

  if ( !defined $sref || $send >= length($$sref) ) {
    # outside our range, so call the segment siever.
    my $seg_ref = _sieve_segment($low, $high);
    return $count + $$seg_ref =~ tr/0//;
  }
  return $count + substr($$sref, $sbeg, $send-$sbeg+1) =~ tr/0//;
}

sub _lehmer_pi {
  my $x = shift;
  return _sieve_prime_count($x) if $x < 1_000;
  do { require Math::BigFloat; Math::BigFloat->import(); }
    if ref($x) eq 'Math::BigInt';
  my $z = (ref($x) ne 'Math::BigInt')
        ? int(sqrt($x+0.5))
        : int(Math::BigFloat->new($x)->badd(0.5)->bsqrt->bfloor->bstr);
  my $a = _lehmer_pi(int(sqrt($z)+0.5));
  my $b = _lehmer_pi($z);
  my $c = _lehmer_pi(int( (ref($x) ne 'Math::BigInt')
                          ? $x**(1/3)+0.5
                          : Math::BigFloat->new($x)->broot(3)->badd(0.5)->bfloor
                     ));
  ($z, $a, $b, $c) = map { (ref($_) =~ /^Math::Big/) ? _bigint_to_int($_) : $_ }
                     ($z, $a, $b, $c);

  # Generate at least b primes.
  my $bth_prime_upper = ($b <= 10) ? 29 : int($b*(log($b) + log(log($b)))) + 1;
  my $primes = primes( $bth_prime_upper );

  my $sum = int(($b + $a - 2) * ($b - $a + 1) / 2);
  $sum += legendre_phi($x, $a, $primes);

  # Get a big sieve for our primecounts.  The C code compromises with either
  # b*10 or x^3/5, as that cuts out all the inner loop sieves and about half
  # of the big outer loop counts.
  # Our sieve count isn't nearly as optimized here, so error on the side of
  # more primes.  This uses a lot more memory but saves a lot of time.
  my $sref = _sieve_erat( int($x / $primes->[$a] / 5) );

  my ($lastw, $lastwpc) = (0,0);
  foreach my $i (reverse $a+1 .. $b) {
    my $w = int($x / $primes->[$i-1]);
    $lastwpc += _count_with_sieve($sref,$lastw+1, $w);
    $lastw = $w;
    $sum -= $lastwpc;
    #$sum -= _count_with_sieve($sref,$w);
    if ($i <= $c) {
      my $bi = _count_with_sieve($sref,int(sqrt($w)+0.5));
      foreach my $j ($i .. $bi) {
        $sum = $sum - _count_with_sieve($sref,int($w / $primes->[$j-1])) + $j - 1;
      }
    }
  }
  $sum;
}
#############################################################################


sub prime_count {
  my($low,$high) = @_;
  if (!defined $high) {
    $high = $low;
    $low = 2;
  }
  _validate_positive_integer($low);
  _validate_positive_integer($high);

  my $count = 0;

  $count++ if ($low <= 2) && ($high >= 2);   # Count 2
  $low = 3 if $low < 3;

  $low++ if ($low % 2) == 0;   # Make low go to odd number.
  $high-- if ($high % 2) == 0; # Make high go to odd number.
  return $count if $low > $high;

  if (   ref($low) eq 'Math::BigInt' || ref($high) eq 'Math::BigInt'
      || ($high-$low) < 10
      || ($high-$low) < int($low/100_000_000_000) ) {
    # Trial primes seems best.  Needs some tuning.
    my $curprime = next_prime($low-1);
    while ($curprime <= $high) {
      $count++;
      $curprime = next_prime($curprime);
    }
    return $count;
  }

  # TODO: Needs tuning
  if ($high > 50_000) {
    if ( ($high / ($high-$low+1)) < 100 ) {
      $count += _lehmer_pi($high);
      $count -= ($low == 3) ? 1 : _lehmer_pi($low-1);
      return $count;
    }
  }

  return (_sieve_prime_count($high) - 1 + $count) if $low == 3;

  my $sieveref = _sieve_segment($low,$high);
  $count += $$sieveref =~ tr/0//;
  return $count;
}


sub nth_prime {
  my($n) = @_;
  _validate_positive_integer($n);

  return undef if $n <= 0;  ## no critic qw(ProhibitExplicitReturnUndef)
  return $_primes_small[$n] if $n <= $#_primes_small;

  if ($n > MPU_MAXPRIMEIDX && ref($n) ne 'Math::BigFloat') {
    do { require Math::BigFloat; Math::BigFloat->import(); }
      if !defined $Math::BigFloat::VERSION;
    $n = Math::BigFloat->new("$n")
  }

  my $prime = 0;
  my $count = 1;
  my $start = 3;

  my $logn = log($n);
  my $loglogn = log($logn);
  my $nth_prime_upper = ($n <= 10) ? 29 : int($n*($logn + $loglogn)) + 1;
  if ($nth_prime_upper > 100000) {
    # Use fast Lehmer prime count combined with lower bound to get close.
    my $nth_prime_lower = int($n * ($logn + $loglogn - 1.0 + (($loglogn-2.10)/$logn)));
    $nth_prime_lower-- unless $nth_prime_lower % 2;
    $count = _lehmer_pi($nth_prime_lower);
    $start = $nth_prime_lower + 2;
  }

  {
    # Make sure incr is an even number.
    my $incr = ($n < 1000) ? 1000 : ($n < 10000) ? 10000 : 100000;
    my $sieveref;
    while (1) {
      $sieveref = _sieve_segment($start, $start+$incr);
      my $segcount = $$sieveref =~ tr/0//;
      last if ($count + $segcount) >= $n;
      $count += $segcount;
      $start += $incr+2;
    }
    # Our count is somewhere in this segment.  Need to look for it.
    $prime = $start - 2;
    while ($count < $n) {
      $prime += 2;
      $count++ if !substr($$sieveref, ($prime-$start)>>1, 1);
    }
  }
  $prime;
}

# The nth prime will be less or equal to this number
sub nth_prime_upper {
  my($n) = @_;
  _validate_positive_integer($n);

  return undef if $n <= 0;  ## no critic qw(ProhibitExplicitReturnUndef)
  return $_primes_small[$n] if $n <= $#_primes_small;

  $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45;

  my $flogn  = log($n);
  my $flog2n = log($flogn);  # Note distinction between log_2(n) and log^2(n)

  my $upper;
  if      ($n >= 46254381) {  # Axler 2017 Corollary 1.2
    $upper = $n * ( $flogn  +  $flog2n-1.0  +  (($flog2n-2.00)/$flogn)  -  (($flog2n*$flog2n - 6*$flog2n + 10.667)/(2*$flogn*$flogn)) );
  } elsif ($n >=  8009824) {  # Axler 2013 page viii Korollar G
    $upper = $n * ( $flogn  +  $flog2n-1.0  +  (($flog2n-2.00)/$flogn)  -  (($flog2n*$flog2n - 6*$flog2n + 10.273)/(2*$flogn*$flogn)) );
  } elsif ($n >=  688383) {   # Dusart 2010 page 2
    $upper = $n * ( $flogn  +  $flog2n - 1.0 + (($flog2n-2.00)/$flogn) );
  } elsif ($n >=  178974) {   # Dusart 2010 page 7
    $upper = $n * ( $flogn  +  $flog2n - 1.0 + (($flog2n-1.95)/$flogn) );
  } elsif ($n >=   39017) {   # Dusart 1999 page 14
    $upper = $n * ( $flogn  +  $flog2n - 0.9484 );
  } elsif ($n >=       6) {   # Modified Robin 1983, for 6-39016 only
    $upper = $n * ( $flogn  +  0.6000 * $flog2n );
  } else {
    $upper = $n * ( $flogn  +  $flog2n );
  }

  return int($upper + 1.0);
}

# The nth prime will be greater than or equal to this number
sub nth_prime_lower {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);

  return undef if $n <= 0;  ## no critic qw(ProhibitExplicitReturnUndef)
  return $_primes_small[$n] if $n <= $#_primes_small;

  $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45;

  my $flogn  = log($n);
  my $flog2n = log($flogn);  # Note distinction between log_2(n) and log^2(n)

  # Dusart 1999 page 14, for all n >= 2
  #my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.25)/$flogn));
  # Dusart 2010 page 2, for all n >= 3
  #my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.10)/$flogn));
  # Axler 2013 page viii Korollar I, for all n >= 2
  #my $lower = $n * ($flogn + $flog2n-1.0 + (($flog2n-2.00)/$flogn) - (($flog2n*$flog2n - 6*$flog2n + 11.847)/(2*$flogn*$flogn)) );
  # Axler 2017 Corollary 1.4
  my $lower = $n * ($flogn + $flog2n-1.0 + (($flog2n-2.00)/$flogn) - (($flog2n*$flog2n - 6*$flog2n + 11.508)/(2*$flogn*$flogn)) );

  return int($lower + 0.999999999);
}

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

  return (0,2,3,5,6,8)[$n] if $n <= 5;
  $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45;
  my $t = $n * log($n);

  # Iterator Halley's method until error term grows
  my $old_term = MPU_INFINITY;
  for my $iter (1 .. 10000) {
    my $dn = Math::Prime::Util::LogarithmicIntegral($t) - $n;
    my $term = $dn * log($t) / (1.0 + $dn/(2*$t));
    last if abs($term) >= abs($old_term);
    $old_term = $term;
    $t -= $term;
    last if abs($term) < 1e-6;
  }
  if (ref($t)) {
    $t = $t->bceil->as_int();
    $t = _bigint_to_int($t) if $t->bacmp(BMAX) <= 0;
  } else {
    $t = int($t+0.999999);
  }
  $t;
}
sub _inverse_R {
  my($n) = @_;
  _validate_num($n) || _validate_positive_integer($n);

  return (0,2,3,5,6,8)[$n] if $n <= 5;
  $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45;
  my $t = $n * log($n);

  # Iterator Halley's method until error term grows
  my $old_term = MPU_INFINITY;
  for my $iter (1 .. 10000) {
    my $dn = Math::Prime::Util::RiemannR($t) - $n;
    my $term = $dn * log($t) / (1.0 + $dn/(2*$t));
    last if abs($term) >= abs($old_term);
    $old_term = $term;
    $t -= $term;
    last if abs($term) < 1e-6;
  }
  if (ref($t)) {
    $t = $t->bceil->as_int();
    $t = _bigint_to_int($t) if $t->bacmp(BMAX) <= 0;
  } else {
    $t = int($t+0.999999);
  }
  $t;
}

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

  return undef if $n <= 0;  ## no critic qw(ProhibitExplicitReturnUndef)
  return $_primes_small[$n] if $n <= $#_primes_small;

  # Once past 10^12 or so, inverse_li gives better results.
  return Math::Prime::Util::inverse_li($n) if $n > 1e12;

  $n = _upgrade_to_float($n)
    if ref($n) eq 'Math::BigInt' || $n >= MPU_MAXPRIMEIDX;

  my $flogn  = log($n);
  my $flog2n = log($flogn);

  # Cipolla 1902:
  #    m=0   fn * ( flogn + flog2n - 1 );
  #    m=1   + ((flog2n - 2)/flogn) );
  #    m=2   - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
  #    + O((flog2n/flogn)^3)
  #
  # Shown in Dusart 1999 page 12, as well as other sources such as:
  #   http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf
  # where the main issue you run into is that you're doing polynomial
  # interpolation, so it oscillates like crazy with many high-order terms.
  # Hence I'm leaving it at m=2.

  my $approx = $n * ( $flogn + $flog2n - 1
                      + (($flog2n - 2)/$flogn)
                      - ((($flog2n*$flog2n) - 6*$flog2n + 11) / (2*$flogn*$flogn))
                    );

  # Apply a correction to help keep values close.
  my $order = $flog2n/$flogn;
  $order = $order*$order*$order * $n;

  if    ($n <        259) { $approx += 10.4 * $order; }
  elsif ($n <        775) { $approx +=  6.3 * $order; }
  elsif ($n <       1271) { $approx +=  5.3 * $order; }
  elsif ($n <       2000) { $approx +=  4.7 * $order; }
  elsif ($n <       4000) { $approx +=  3.9 * $order; }
  elsif ($n <      12000) { $approx +=  2.8 * $order; }
  elsif ($n <     150000) { $approx +=  1.2 * $order; }
  elsif ($n <   20000000) { $approx +=  0.11 * $order; }
  elsif ($n <  100000000) { $approx +=  0.008 * $order; }
  elsif ($n <  500000000) { $approx += -0.038 * $order; }
  elsif ($n < 2000000000) { $approx += -0.054 * $order; }
  else                    { $approx += -0.058 * $order; }
  # If we want the asymptotic approximation to be >= actual, use -0.010.

  return int($approx + 0.5);
}

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

sub prime_count_approx {
  my($x) = @_;
  _validate_num($x) || _validate_positive_integer($x);

  # Turn on high precision FP if they gave us a big number.
  $x = _upgrade_to_float($x) if ref($_[0]) eq 'Math::BigInt';
  #    Method             10^10 %error  10^19 %error
  #    -----------------  ------------  ------------
  #    n/(log(n)-1)        .22%          .058%
  #    n/(ln(n)-1-1/ln(n)) .032%         .0041%
  #    average bounds      .0005%        .0000002%
  #    asymp               .0006%        .00000004%
  #    li(n)               .0007%        .00000004%
  #    li(n)-li(n^.5)/2    .0004%        .00000001%
  #    R(n)                .0004%        .00000001%
  #
  # Also consider: http://trac.sagemath.org/sage_trac/ticket/8135

  # my $result = int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
  # my $result = int( LogarithmicIntegral($x) );
  # my $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2);
  # my $result = RiemannR($x) + 0.5;

  # Sadly my Perl RiemannR function is really slow for big values.
  # However I have written versions in GMP (mpf) and MPFR.  If those are
  # available then we will use them.
  # Otherwise, switch to LiCorr for very big values.  This is hacky and
  # shouldn't be necessary.
  my $result;
  if ( $x < 1e36 || _MPFR_available() || $Math::Prime::Util::_GMPfunc{"riemannr"} ) {
    if (ref($x) eq 'Math::BigFloat') {
      # Make sure we get enough accuracy, and also not too much more than needed
      $x->accuracy(length($x->copy->as_int->bstr())+2);
    }
    $result = RiemannR($x) + 0.5;
  } else {
    # Math::BigInt's default Calc backend takes *ages* to do a cube root, so
    # limit ourselves to just the first two terms.
    $result = int(
        LogarithmicIntegral($x)
      - LogarithmicIntegral(sqrt($x))/2
    #  - LogarithmicIntegral($x**(1.0/3.0))/3
    #  - LogarithmicIntegral($x**(1.0/5.0))/5
    #  + LogarithmicIntegral($x**(1.0/6.0))/6
    #  - LogarithmicIntegral($x**(1.0/7.0))/7
    #  ...
    );
  }
  # Asymp:
  #   my $l1 = log($x);  my $l2 = $l1*$l1;  my $l4 = $l2*$l2;
  #   $result = int( $x/$l1 + $x/$l2 + 2*$x/($l2*$l1) + 6*$x/($l4) + 24*$x/($l4*$l1) + 120*$x/($l4*$l2) + 720*$x/($l4*$l2*$l1) + 5040*$x/($l4*$l4) + 40320*$x/($l4*$l4*$l1) + 0.5 );

  return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
  return int($result);
}

sub prime_count_lower {
  my($x) = @_;
  _validate_num($x) || _validate_positive_integer($x);

  return _tiny_prime_count($x) if $x < $_primes_small[-1];

  $x = _upgrade_to_float($x)
    if ref($x) eq 'Math::BigInt' || ref($_[0]) eq 'Math::BigInt';

  my($result,$a);
  my $fl1 = log($x);
  my $fl2 = $fl1*$fl1;
  my $one = (ref($x) eq 'Math::BigFloat') ? $x->copy->bone : $x-$x+1.0;

  # Chebyshev            1*x/logx       x >= 17
  # Rosser & Schoenfeld  x/(logx-1/2)   x >= 67
  # Dusart 1999          x/logx*(1+1/logx+1.8/logxlogx)  x >= 32299
  # Dusart 2010          x/logx*(1+1/logx+2.0/logxlogx)  x >= 88783
  # Axler 2014 (1.2)     ""+...                          x >= 1332450001
  # Axler 2014 (1.2)     x/(logx-1-1/logx-...)           x >= 1332479531
  # Büthe 2015 (1.9)     li(x)-(sqrtx/logx)*(...)        x <= 10^19
  # Büthe 2014 Th 2      li(x)-logx*sqrtx/8Pi    x > 2657, x <= 1.4*10^25

  if ($x < 599) {                         # Decent for small numbers
    $result = $x / ($fl1 - 0.7);
  } elsif ($x < 52600000) {               # Dusart 2010 tweaked
    if    ($x <       2700) { $a = 0.30; }
    elsif ($x <       5500) { $a = 0.90; }
    elsif ($x <      19400) { $a = 1.30; }
    elsif ($x <      32299) { $a = 1.60; }
    elsif ($x <      88783) { $a = 1.83; }
    elsif ($x <     176000) { $a = 1.99; }
    elsif ($x <     315000) { $a = 2.11; }
    elsif ($x <    1100000) { $a = 2.19; }
    elsif ($x <    4500000) { $a = 2.31; }
    else                    { $a = 2.35; }
    $result = ($x/$fl1) * ($one + $one/$fl1 + $a/$fl2);
  } elsif ($x < 1.4e25 || Math::Prime::Util::prime_get_config()->{'assume_rh'}){
                                          # Büthe 2014/2015
    if (_MPFR_available()) {
      my $wantbf = (defined $bignum::VERSION || ref($x) =~ /^Math::Big/);
      my $xdigits = length($x);
      $xdigits += length(int(log(0.0+"$x"))) + 1;
      my $rnd = 0;  # MPFR_RNDN;
      my $bit_precision = int($xdigits * 3.322) + 4;
      my $rx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
      Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
      my $lix = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($lix, $bit_precision);
      Math::MPFR::Rmpfr_set_str($lix, LogarithmicIntegral($x,1),10,$rnd);
      my $sqx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($sqx, $bit_precision);
      Math::MPFR::Rmpfr_sqrt($sqx, $rx, $bit_precision);
      Math::MPFR::Rmpfr_log($rx, $rx, $rnd);
      # rx = log(x) lix = li(x) sqx = sqrt(x)
      if ($x < 1e19) {                    # Büthe 2015 1.9
        #Math::MPFR::Rmpfr_div($sqx, $sqx, $rx, $rnd);
        #$rx = 1.94 + 3.88/$rx + 27.57/($rx*$rx);
      my $tmp = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($tmp, $bit_precision);
      my $acc = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($acc, $bit_precision);
      Math::MPFR::Rmpfr_set_d($acc, 1.94, $rnd);
      Math::MPFR::Rmpfr_d_div($tmp, 3.88, $rx, $rnd);
      Math::MPFR::Rmpfr_add($acc, $acc, $tmp, $rnd);
      Math::MPFR::Rmpfr_d_div($tmp, 27.57, $rx*$rx, $rnd);
      Math::MPFR::Rmpfr_add($acc, $acc, $tmp, $rnd);
      Math::MPFR::Rmpfr_mul($rx, $acc, $sqx/$rx, $rnd);
        #Math::MPFR::Rmpfr_mul($rx, $rx, $sqx, $rnd);
      } else {                            # Büthe 2014 7.4
        Math::MPFR::Rmpfr_mul($rx, $rx, $sqx, $rnd);
        Math::MPFR::Rmpfr_const_pi($sqx, $rnd);
        Math::MPFR::Rmpfr_mul_ui($sqx, $sqx, 8, $rnd);
        Math::MPFR::Rmpfr_div($rx, $rx, $sqx, $rnd);
      }
      Math::MPFR::Rmpfr_sub($rx, $lix, $rx, $rnd);
      my $strval = Math::MPFR::Rmpfr_integer_string($rx, 10, $rnd);
      $result = ($wantbf) ? Math::BigInt->new($strval) : int($strval);
    } else {
      my $lix = LogarithmicIntegral($x);
      my $sqx = sqrt($x);
      if ($x < 1e19) {
        $result = $lix - ($sqx/$fl1) * (1.94 + 3.88/$fl1 + 27.57/$fl2);
      } else {
        if (ref($x) eq 'Math::BigFloat') {
          my $xdigits = _find_big_acc($x);
          $result = $lix - ($fl1*$sqx / (Math::BigFloat->bpi($xdigits)*8));
        } else {
          $result = $lix - ($fl1*$sqx / PI_TIMES_8);
        }
      }
    }
  } else {                                # Axler 2014 1.4
    if (_MPFR_available()) {
      my $wantbf = (defined $bignum::VERSION || ref($x) =~ /^Math::Big/);
      my $xdigits = length($x);
      $xdigits += length(int(log(0.0+"$x"))) + 1;
      my $rnd = 0;  # MPFR_RNDN;
      my $bit_precision = int($xdigits * 3.322) + 4;
      my $rx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
      Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
      my $term = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($term, $bit_precision);
      my $logx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($logx, $bit_precision);
      my $loglogx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($loglogx, $bit_precision);
      my $div = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($div, $bit_precision);
      Math::MPFR::Rmpfr_log($logx, $rx, $rnd);
      Math::MPFR::Rmpfr_set_ui($loglogx, 1, $rnd);
      Math::MPFR::Rmpfr_sub_ui($div, $logx, 1, $rnd);

      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 1.0, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 2.65, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 13.35, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 70.3, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 465.6275, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 3404.4225, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);

      Math::MPFR::Rmpfr_div($rx, $rx, $div, $rnd);
      my $strval = Math::MPFR::Rmpfr_integer_string($rx, 10, $rnd);
      $result = ($wantbf) ? Math::BigInt->new($strval) : int($strval);
    } else {
      my($fl3,$fl4) = ($fl2*$fl1,$fl2*$fl2);
      my($fl5,$fl6) = ($fl4*$fl1,$fl4*$fl2);
      $result = $x / ($fl1 - $one - $one/$fl1 - 2.65/$fl2 - 13.35/$fl3 - 70.3/$fl4 - 455.6275/$fl5 - 3404.4225/$fl6);
    }
  }

  return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
  return int($result);
}

sub prime_count_upper {
  my($x) = @_;
  _validate_num($x) || _validate_positive_integer($x);

  # Give an exact answer for what we have in our little table.
  return _tiny_prime_count($x) if $x < $_primes_small[-1];

  $x = _upgrade_to_float($x)
    if ref($x) eq 'Math::BigInt' || ref($_[0]) eq 'Math::BigInt';

  # Chebyshev:            1.25506*x/logx       x >= 17
  # Rosser & Schoenfeld:  x/(logx-3/2)         x >= 67
  # Panaitopol 1999:      x/(logx-1.112)       x >= 4
  # Dusart 1999:          x/logx*(1+1/logx+2.51/logxlogx)   x >= 355991
  # Dusart 2010:          x/logx*(1+1/logx+2.334/logxlogx)  x >= 2_953_652_287
  # Axler 2014:           x/(logx-1-1/logx-3.35/logxlogx...) x >= e^3.804
  # Büthe 2014 7.4        Schoenfeld bounds hold to x <= 1.4e25
  # Axler 2017 Prop 2.2   Schoenfeld bounds hold to x <= 5.5e25
  # Skewes                li(x)                x < 1e14

  my($result,$a);
  my $fl1 = log($x);
  my $fl2 = $fl1 * $fl1;
  my $one = (ref($x) eq 'Math::BigFloat') ? $x->copy->bone : $x-$x+1.0;

  if ($x < 15900) {              # Tweaked Rosser-type
    $a = ($x < 1621) ? 1.048 : ($x < 5000) ? 1.071 : 1.098;
    $result = ($x / ($fl1 - $a)) + 1.0;
  } elsif ($x < 821800000) {     # Tweaked Dusart 2010
    if    ($x <      24000) { $a = 2.30; }
    elsif ($x <      59000) { $a = 2.48; }
    elsif ($x <     350000) { $a = 2.52; }
    elsif ($x <     355991) { $a = 2.54; }
    elsif ($x <     356000) { $a = 2.51; }
    elsif ($x <    3550000) { $a = 2.50; }
    elsif ($x <    3560000) { $a = 2.49; }
    elsif ($x <    5000000) { $a = 2.48; }
    elsif ($x <    8000000) { $a = 2.47; }
    elsif ($x <   13000000) { $a = 2.46; }
    elsif ($x <   18000000) { $a = 2.45; }
    elsif ($x <   31000000) { $a = 2.44; }
    elsif ($x <   41000000) { $a = 2.43; }
    elsif ($x <   48000000) { $a = 2.42; }
    elsif ($x <  119000000) { $a = 2.41; }
    elsif ($x <  182000000) { $a = 2.40; }
    elsif ($x <  192000000) { $a = 2.395; }
    elsif ($x <  213000000) { $a = 2.390; }
    elsif ($x <  271000000) { $a = 2.385; }
    elsif ($x <  322000000) { $a = 2.380; }
    elsif ($x <  400000000) { $a = 2.375; }
    elsif ($x <  510000000) { $a = 2.370; }
    elsif ($x <  682000000) { $a = 2.367; }
    elsif ($x < 2953652287) { $a = 2.362; }
    else                    { $a = 2.334; } # Dusart 2010, page 2
    $result = ($x/$fl1) * ($one + $one/$fl1 + $a/$fl2) + $one;
  } elsif ($x < 1e19) {                     # Skewes number lower limit
    $a = ($x < 110e7) ? 0.032 : ($x < 1001e7) ? 0.027 : ($x < 10126e7) ? 0.021 : 0.0;
    $result = LogarithmicIntegral($x) - $a * $fl1*sqrt($x)/PI_TIMES_8;
  } elsif ($x < 5.5e25 || Math::Prime::Util::prime_get_config()->{'assume_rh'}) {
                                            # Schoenfeld / Büthe 2014 Th 7.4
    if (_MPFR_available()) {
      my $wantbf = (defined $bignum::VERSION || ref($x) =~ /^Math::Big/);
      my $xdigits = length($x);
      $xdigits += length(int(log(0.0+"$x"))) + 1;
      my $rnd = 0;  # MPFR_RNDN;
      my $bit_precision = int($xdigits * 3.322) + 4;
      my $rx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
      Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
      my $lix = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($lix, $bit_precision);
      Math::MPFR::Rmpfr_set_str($lix, LogarithmicIntegral($x,1),10,$rnd);
      my $sqx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($sqx, $bit_precision);
      Math::MPFR::Rmpfr_sqrt($sqx, $rx, $bit_precision);
      Math::MPFR::Rmpfr_log($rx, $rx, $rnd);
      Math::MPFR::Rmpfr_mul($rx, $rx, $sqx, $rnd);
      Math::MPFR::Rmpfr_const_pi($sqx, $rnd);
      Math::MPFR::Rmpfr_mul_ui($sqx, $sqx, 8, $rnd);
      Math::MPFR::Rmpfr_div($rx, $rx, $sqx, $rnd);
      Math::MPFR::Rmpfr_add($rx, $lix, $rx, $rnd);
      my $strval = Math::MPFR::Rmpfr_integer_string($rx, 10, $rnd);
      $result = ($wantbf) ? Math::BigInt->new($strval) : int($strval);
    } else {
      my $lix = LogarithmicIntegral($x);
      my $sqx = sqrt($x);
      if (ref($x) eq 'Math::BigFloat') {
        my $xdigits = _find_big_acc($x);
        $result = $lix + ($fl1*$sqx / (Math::BigFloat->bpi($xdigits)*8));
      } else {
        $result = $lix + ($fl1*$sqx / PI_TIMES_8);
      }
    }
  } else {                                  # Axler 2014 1.3
    if (_MPFR_available()) {
      my $wantbf = (defined $bignum::VERSION || ref($x) =~ /^Math::Big/);
      my $xdigits = length($x);
      $xdigits += length(int(log(0.0+"$x"))) + 1;
      my $rnd = 0;  # MPFR_RNDN;
      my $bit_precision = int($xdigits * 3.322) + 4;
      my $rx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
      Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
      my $term = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($term, $bit_precision);
      my $logx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($logx, $bit_precision);
      my $loglogx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($loglogx, $bit_precision);
      my $div = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($div, $bit_precision);
      Math::MPFR::Rmpfr_log($logx, $rx, $rnd);
      Math::MPFR::Rmpfr_set_ui($loglogx, 1, $rnd);
      Math::MPFR::Rmpfr_sub_ui($div, $logx, 1, $rnd);

      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 1.0, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 3.35, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 12.65, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 71.7, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 466.1275, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);
      Math::MPFR::Rmpfr_mul($loglogx, $loglogx, $logx, $rnd);
      Math::MPFR::Rmpfr_d_div($term, 3489.8225, $loglogx, $rnd);
      Math::MPFR::Rmpfr_sub($div, $div, $term, $rnd);

      Math::MPFR::Rmpfr_div($rx, $rx, $div, $rnd);
      my $strval = Math::MPFR::Rmpfr_integer_string($rx, 10, $rnd);
      $result = ($wantbf) ? Math::BigInt->new($strval) : int($strval);
    } else {
      my($fl3,$fl4) = ($fl2*$fl1,$fl2*$fl2);
      my($fl5,$fl6) = ($fl4*$fl1,$fl4*$fl2);
      $result = $x / ($fl1 - $one - $one/$fl1 - 3.35/$fl2 - 12.65/$fl3 - 71.7/$fl4 - 466.1275/$fl5 - 3489.8225/$fl6);
    }
  }

  return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat';
  return int($result);
}

sub twin_prime_count {
  my($low,$high) = @_;
  if (defined $high) { _validate_positive_integer($low); }
  else               { ($low,$high) = (2, $low);         }
  _validate_positive_integer($high);
  my $sum = 0;
  while ($low <= $high) {
    my $seghigh = ($high-$high) + $low + 1e7 - 1;
    $seghigh = $high if $seghigh > $high;
    $sum += scalar(@{Math::Prime::Util::twin_primes($low,$seghigh)});
    $low = $seghigh + 1;
  }
  $sum;
}
sub ramanujan_prime_count {
  my($low,$high) = @_;
  if (defined $high) { _validate_positive_integer($low); }
  else               { ($low,$high) = (2, $low);         }
  _validate_positive_integer($high);
  my $sum = 0;
  while ($low <= $high) {
    my $seghigh = ($high-$high) + $low + 1e9 - 1;
    $seghigh = $high if $seghigh > $high;
    $sum += scalar(@{Math::Prime::Util::ramanujan_primes($low,$seghigh)});
    $low = $seghigh + 1;
  }
  $sum;
}

sub twin_prime_count_approx {
  my($n) = @_;
  return twin_prime_count(3,$n) if $n < 2000;
  $n = _upgrade_to_float($n) if ref($n);
  my $logn = log($n);
  # The loss of full Ei precision is a few orders of magnitude less than the
  # accuracy of the estimate, so save huge time and don't bother.
  my $li2 = Math::Prime::Util::ExponentialIntegral("$logn") + 2.8853900817779268147198494 - ($n/$logn);

  # Empirical correction factor
  my $fm;
  if    ($n <     4000) { $fm = 0.2952; }
  elsif ($n <     8000) { $fm = 0.3151; }
  elsif ($n <    16000) { $fm = 0.3090; }
  elsif ($n <    32000) { $fm = 0.3096; }
  elsif ($n <    64000) { $fm = 0.3100; }
  elsif ($n <   128000) { $fm = 0.3089; }
  elsif ($n <   256000) { $fm = 0.3099; }
  elsif ($n <   600000) { my($x0, $x1, $y0, $y1) = (1e6, 6e5, .3091, .3059);
                          $fm = $y0 + ($n - $x0) * ($y1-$y0) / ($x1 - $x0); }
  elsif ($n <  1000000) { my($x0, $x1, $y0, $y1) = (6e5, 1e6, .3062, .3042);
                          $fm = $y0 + ($n - $x0) * ($y1-$y0) / ($x1 - $x0); }
  elsif ($n <  4000000) { my($x0, $x1, $y0, $y1) = (1e6, 4e6, .3067, .3041);
                          $fm = $y0 + ($n - $x0) * ($y1-$y0) / ($x1 - $x0); }
  elsif ($n < 16000000) { my($x0, $x1, $y0, $y1) = (4e6, 16e6, .3033, .2983);
                          $fm = $y0 + ($n - $x0) * ($y1-$y0) / ($x1 - $x0); }
  elsif ($n < 32000000) { my($x0, $x1, $y0, $y1) = (16e6, 32e6, .2980, .2965);
                          $fm = $y0 + ($n - $x0) * ($y1-$y0) / ($x1 - $x0); }
  $li2 *= $fm * log(12+$logn)  if defined $fm;

  return int(1.32032363169373914785562422 * $li2 + 0.5);
}

sub nth_twin_prime {
  my($n) = @_;
  return undef if $n < 0;  ## no critic qw(ProhibitExplicitReturnUndef)
  return (undef,3,5,11,17,29,41)[$n] if $n <= 6;

  my $p = Math::Prime::Util::nth_twin_prime_approx($n+200);
  my $tp = Math::Prime::Util::twin_primes($p);
  while ($n > scalar(@$tp)) {
    $n -= scalar(@$tp);
    $tp = Math::Prime::Util::twin_primes($p+1,$p+1e5);
    $p += 1e5;
  }
  return $tp->[$n-1];
}

sub nth_twin_prime_approx {
  my($n) = @_;
  _validate_positive_integer($n);
  return nth_twin_prime($n) if $n < 6;
  $n = _upgrade_to_float($n) if ref($n) || $n > 127e14;   # TODO lower for 32-bit
  my $logn = log($n);
  my $nlogn2 = $n * $logn * $logn;

  return int(5.158 * $nlogn2/log(9+log($n*$n))) if $n > 59 && $n <= 1092;

  my $lo = int(0.7 * $nlogn2);
  my $hi = int( ($n > 1e16) ? 1.1 * $nlogn2
              : ($n >  480) ? 1.7 * $nlogn2
                            : 2.3 * $nlogn2 + 3 );

  _binary_search($n, $lo, $hi,
                 sub{Math::Prime::Util::twin_prime_count_approx(shift)},
                 sub{ ($_[2]-$_[1])/$_[1] < 1e-15 } );
}

sub nth_ramanujan_prime_upper {
  my $n = shift;
  return (0,2,11)[$n] if $n <= 2;
  $n = Math::BigInt->new("$n") if $n > (~0/3);
  my $nth = nth_prime_upper(3*$n);
  return $nth if $n < 10000;
  $nth = Math::BigInt->new("$nth") if $nth > (~0/177);
  if ($n < 1000000) { $nth = (177 * $nth) >> 8; }
  elsif ($n < 1e10) { $nth = (175 * $nth) >> 8; }
  else              { $nth = (133 * $nth) >> 8; }
  $nth = _bigint_to_int($nth) if ref($nth) && $nth->bacmp(BMAX) <= 0;
  $nth;
}
sub nth_ramanujan_prime_lower {
  my $n = shift;
  return (0,2,11)[$n] if $n <= 2;
  $n = Math::BigInt->new("$n") if $n > (~0/2);
  my $nth = nth_prime_lower(2*$n);
  $nth = Math::BigInt->new("$nth") if $nth > (~0/275);
  if ($n < 10000)   { $nth = (275 * $nth) >> 8; }
  elsif ($n < 1e10) { $nth = (262 * $nth) >> 8; }
  $nth = _bigint_to_int($nth) if ref($nth) && $nth->bacmp(BMAX) <= 0;
  $nth;
}
sub nth_ramanujan_prime_approx {
  my $n = shift;
  return (0,2,11)[$n] if $n <= 2;
  my($lo,$hi) = (nth_ramanujan_prime_lower($n),nth_ramanujan_prime_upper($n));
  $lo + (($hi-$lo)>>1);
}
sub ramanujan_prime_count_upper {
  my $n = shift;
  return (($n < 2) ? 0 : 1) if $n < 11;
  my $lo = int(prime_count_lower($n) / 3);
  my $hi = prime_count_upper($n) >> 1;
  1+_binary_search($n, $lo, $hi,
                   sub{Math::Prime::Util::nth_ramanujan_prime_lower(shift)});
}
sub ramanujan_prime_count_lower {
  my $n = shift;
  return (($n < 2) ? 0 : 1) if $n < 11;
  my $lo = int(prime_count_lower($n) / 3);
  my $hi = prime_count_upper($n) >> 1;
  _binary_search($n, $lo, $hi,
                 sub{Math::Prime::Util::nth_ramanujan_prime_upper(shift)});
}
sub ramanujan_prime_count_approx {
  my $n = shift;
  return (($n < 2) ? 0 : 1) if $n < 11;
  #$n = _upgrade_to_float($n) if ref($n) || $n > 2e16;
  my $lo = ramanujan_prime_count_lower($n);
  my $hi = ramanujan_prime_count_upper($n);
  _binary_search($n, $lo, $hi,
                 sub{Math::Prime::Util::nth_ramanujan_prime_approx(shift)},
                 sub{ ($_[2]-$_[1])/$_[1] < 1e-15 } );
}

sub _sum_primes_n {
  my $n = shift;
  return (0,0,2,5,5)[$n] if $n < 5;
  my $r = Math::Prime::Util::sqrtint($n);
  my $r2 = $r + int($n/($r+1));
  my(@V,@S);
  for my $k (0 .. $r2) {
    my $v = ($k <= $r) ? $k : int($n/($r2-$k+1));
    $V[$k] = $v;
    $S[$k] = (($v*($v+1)) >> 1) - 1;
  }
  Math::Prime::Util::forprimes( sub { my $p = $_;
    my $sp = $S[$p-1];
    my $p2 = $p*$p;
    for my $v (reverse @V) {
      last if $v < $p2;
      my($a,$b) = ($v,int($v/$p));
      $a = $r2 - int($n/$a) + 1 if $a > $r;
      $b = $r2 - int($n/$b) + 1 if $b > $r;
      $S[$a] -= $p * ($S[$b] - $sp);
    }
  }, 2, $r);
  $S[$r2];
}

sub sum_primes {
  my($low,$high) = @_;
  if (defined $high) { _validate_positive_integer($low); }
  else               { ($low,$high) = (2, $low);         }
  _validate_positive_integer($high);
  my $sum = 0;
  $sum = BZERO->copy if ( (MPU_32BIT && $high >        323_380) ||
                          (MPU_64BIT && $high > 29_505_444_490) );

  # It's very possible we're here because they've counted too high.  Skip fwd.
  if ($low <= 2 && $high >= 29505444491) {
    $low = 29505444503;
    $sum = Math::BigInt->new("18446744087046669523");
  }

  return $sum if $low > $high;

  # We have to make some decision about whether to use our PP prime sum or loop
  # doing the XS sieve.  TODO: Be smarter here?
  if (!Math::Prime::Util::prime_get_config()->{'xs'} && !ref($sum) && !MPU_32BIT && ($high-$low) > 1000000) {
    # Unfortunately with bigints this is horrifically slow, but we have to do it.
    $high = BZERO->copy + $high  if $high >= (1 << (MPU_MAXBITS/2))-1;
    $sum = _sum_primes_n($high);
    $sum -= _sum_primes_n($low-1) if $low > 2;
    return $sum;
  }

  my $xssum = (MPU_64BIT && $high < 6e14 && Math::Prime::Util::prime_get_config()->{'xs'});
  my $step = ($xssum && $high > 5e13) ? 1_000_000 : 11_000_000;
  Math::Prime::Util::prime_precalc(sqrtint($high));
  while ($low <= $high) {
    my $next = $low + $step - 1;
    $next = $high if $next > $high;
    $sum += ($xssum) ? Math::Prime::Util::sum_primes($low,$next)
                     : Math::Prime::Util::vecsum( @{Math::Prime::Util::primes($low,$next)} );
    last if $next == $high;
    $low = $next+1;
  }
  $sum;
}
sub print_primes {
  my($low,$high,$fd) = @_;
  if (defined $high) { _validate_positive_integer($low); }
  else               { ($low,$high) = (2, $low);         }
  _validate_positive_integer($high);

  $fd = fileno(STDOUT) unless defined $fd;
  open(my $fh, ">>&=", $fd);  # TODO .... or die

  if ($high >= $low) {
    my $p1 = $low;
    while ($p1 <= $high) {
      my $p2 = $p1 + 15_000_000 - 1;
      $p2 = $high if $p2 > $high;
      if ($Math::Prime::Util::_GMPfunc{"sieve_primes"}) {
        print $fh "$_\n" for Math::Prime::Util::GMP::sieve_primes($p1,$p2,0);
      } else {
        print $fh "$_\n" for @{primes($p1,$p2)};
      }
      $p1 = $p2+1;
    }
  }
  close($fh);
}


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

sub _mulmod {
  my($x, $y, $n) = @_;
  return (($x * $y) % $n) if ($x|$y) < MPU_HALFWORD;
  #return (($x * $y) % $n) if ($x|$y) < MPU_HALFWORD || $y == 0 || $x < int(~0/$y);
  my $r = 0;
  $x %= $n if $x >= $n;
  $y %= $n if $y >= $n;
  ($x,$y) = ($y,$x) if $x < $y;
  if ($n <= (~0 >> 1)) {
    while ($y > 1) {
      if ($y & 1) { $r += $x;  $r -= $n if $r >= $n; }
      $y >>= 1;
      $x += $x;  $x -= $n if $x >= $n;
    }
    if ($y & 1) { $r += $x;  $r -= $n if $r >= $n; }
  } else {
    while ($y > 1) {
      if ($y & 1) { $r = $n-$r;  $r = ($x >= $r) ? $x-$r : $n-$r+$x; }
      $y >>= 1;
      $x = ($x > ($n - $x))  ?  ($x - $n) + $x  :  $x + $x;
    }
    if ($y & 1) { $r = $n-$r;  $r = ($x >= $r) ? $x-$r : $n-$r+$x; }
  }
  $r;
}
sub _addmod {
  my($x, $y, $n) = @_;
  $x %= $n if $x >= $n;
  $y %= $n if $y >= $n;
  if (($n-$x) <= $y) {
    ($x,$y) = ($y,$x) if $y > $x;
    $x -= $n;
  }
  $x + $y;
}

# Note that Perl 5.6.2 with largish 64-bit numbers will break.  As usual.
sub _native_powmod {
  my($n, $power, $m) = @_;
  my $t = 1;
  $n = $n % $m;
  while ($power) {
    $t = ($t * $n) % $m if ($power & 1);
    $power >>= 1;
    $n = ($n * $n) % $m if $power;
  }
  $t;
}

sub _powmod {
  my($n, $power, $m) = @_;
  my $t = 1;

  $n %= $m if $n >= $m;
  if ($m < MPU_HALFWORD) {
    while ($power) {
      $t = ($t * $n) % $m if ($power & 1);
      $power >>= 1;
      $n = ($n * $n) % $m if $power;
    }
  } else {
    while ($power) {
      $t = _mulmod($t, $n, $m) if ($power & 1);
      $power >>= 1;
      $n = _mulmod($n, $n, $m) if $power;
    }
  }
  $t;
}

# Make sure to work around RT71548, Math::BigInt::Lite,
# and use correct lcm semantics.
sub gcd {
  # First see if all inputs are non-bigints  5-10x faster if so.
  if (0 == scalar(grep { ref($_) } @_)) {
    my($x,$y) = (shift || 0, 0);
    while (@_) {
      $y = shift;
      while ($y) {  ($x,$y) = ($y, $x % $y);  }
      $x = -$x if $x < 0;
    }
    return $x;
  }
  my $gcd = Math::BigInt::bgcd( map {
    my $v = (($_ < 2147483647 && !ref($_)) || ref($_) eq 'Math::BigInt') ? $_ : "$_";
    $v;
  } @_ );
  $gcd = _bigint_to_int($gcd) if $gcd->bacmp(BMAX) <= 0;
  return $gcd;
}
sub lcm {
  return 0 unless @_;
  my $lcm = Math::BigInt::blcm( map {
    my $v = (($_ < 2147483647 && !ref($_)) || ref($_) eq 'Math::BigInt') ? $_ : "$_";
    return 0 if $v == 0;
    $v = -$v if $v < 0;
    $v;
  } @_ );
  $lcm = _bigint_to_int($lcm) if $lcm->bacmp(BMAX) <= 0;
  return $lcm;
}
sub gcdext {
  my($x,$y) = @_;
  if ($x == 0) { return (0, (-1,0,1)[($y>=0)+($y>0)], abs($y)); }
  if ($y == 0) { return ((-1,0,1)[($x>=0)+($x>0)], 0, abs($x)); }

  if ($Math::Prime::Util::_GMPfunc{"gcdext"}) {
    my($a,$b,$g) = Math::Prime::Util::GMP::gcdext($x,$y);
    $a = Math::Prime::Util::_reftyped($_[0], $a);
    $b = Math::Prime::Util::_reftyped($_[0], $b);
    $g = Math::Prime::Util::_reftyped($_[0], $g);
    return ($a,$b,$g);
  }

  my($a,$b,$g,$u,$v,$w);
  if (abs($x) < (~0>>1) && abs($y) < (~0>>1)) {
    $x = _bigint_to_int($x) if ref($x) eq 'Math::BigInt';
    $y = _bigint_to_int($y) if ref($y) eq 'Math::BigInt';
    ($a,$b,$g,$u,$v,$w) = (1,0,$x,0,1,$y);
    while ($w != 0) {
      my $r = $g % $w;
      my $q = int(($g-$r)/$w);
      ($a,$b,$g,$u,$v,$w) = ($u,$v,$w,$a-$q*$u,$b-$q*$v,$r);
    }
  } else {
    ($a,$b,$g,$u,$v,$w) = (BONE->copy,BZERO->copy,Math::BigInt->new("$x"),
                           BZERO->copy,BONE->copy,Math::BigInt->new("$y"));
    while ($w != 0) {
      # Using the array bdiv is logical, but is the wrong sign.
      my $r = $g->copy->bmod($w);
      my $q = $g->copy->bsub($r)->bdiv($w);
      ($a,$b,$g,$u,$v,$w) = ($u,$v,$w,$a-$q*$u,$b-$q*$v,$r);
    }
    $a = _bigint_to_int($a) if $a->bacmp(BMAX) <= 0;
    $b = _bigint_to_int($b) if $b->bacmp(BMAX) <= 0;
    $g = _bigint_to_int($g) if $g->bacmp(BMAX) <= 0;
  }
  if ($g < 0) { ($a,$b,$g) = (-$a,-$b,-$g); }
  return ($a,$b,$g);
}

sub chinese {
  return 0 unless scalar @_;
  return $_[0]->[0] % $_[0]->[1] if scalar @_ == 1;
  my($lcm, $sum);

  if ($Math::Prime::Util::_GMPfunc{"chinese"} && $Math::Prime::Util::GMP::VERSION >= 0.42) {
    $sum = Math::Prime::Util::GMP::chinese(@_);
    if (defined $sum) {
      $sum = Math::BigInt->new("$sum");
      $sum = _bigint_to_int($sum) if ref($sum) && $sum->bacmp(BMAX) <= 0;
    }
    return $sum;
  }
  foreach my $aref (sort { $b->[1] <=> $a->[1] } @_) {
    my($ai, $ni) = @$aref;
    $ai = Math::BigInt->new("$ai") if !ref($ai) && (abs($ai) > (~0>>1) || OLD_PERL_VERSION);
    $ni = Math::BigInt->new("$ni") if !ref($ni) && (abs($ni) > (~0>>1) || OLD_PERL_VERSION);
    if (!defined $lcm) {
      ($sum,$lcm) = ($ai % $ni, $ni);
      next;
    }
    # gcdext
    my($u,$v,$g,$s,$t,$w) = (1,0,$lcm,0,1,$ni);
    while ($w != 0) {
      my $r = $g % $w;
      my $q = ref($g)  ?  $g->copy->bsub($r)->bdiv($w)  :  int(($g-$r)/$w);
      ($u,$v,$g,$s,$t,$w) = ($s,$t,$w,$u-$q*$s,$v-$q*$t,$r);
    }
    ($u,$v,$g) = (-$u,-$v,-$g)  if $g < 0;
    return if $g != 1 && ($sum % $g) != ($ai % $g);  # Not co-prime
    $s = -$s if $s < 0;
    $t = -$t if $t < 0;
    # Convert to bigint if necessary.  Performance goes to hell.
    if (!ref($lcm) && ($lcm*$s) > ~0) { $lcm = Math::BigInt->new("$lcm"); }
    if (ref($lcm)) {
      $lcm->bmul("$s");
      my $m1 = Math::BigInt->new("$v")->bmul("$s")->bmod($lcm);
      my $m2 = Math::BigInt->new("$u")->bmul("$t")->bmod($lcm);
      $m1->bmul("$sum")->bmod($lcm);
      $m2->bmul("$ai")->bmod($lcm);
      $sum = $m1->badd($m2)->bmod($lcm);
    } else {
      $lcm *= $s;
      $u += $lcm if $u < 0;
      $v += $lcm if $v < 0;
      my $vs = _mulmod($v,$s,$lcm);
      my $ut = _mulmod($u,$t,$lcm);
      my $m1 = _mulmod($sum,$vs,$lcm);
      my $m2 = _mulmod($ut,$ai % $lcm,$lcm);
      $sum = _addmod($m1, $m2, $lcm);
    }
  }
  $sum = _bigint_to_int($sum) if ref($sum) && $sum->bacmp(BMAX) <= 0;
  $sum;
}

sub _from_128 {
  my($hi, $lo) = @_;
  return 0 unless defined $hi && defined $lo;
  #print "hi $hi lo $lo\n";
  (Math::BigInt->new("$hi") << MPU_MAXBITS) + $lo;
}

sub vecsum {
  return Math::Prime::Util::_reftyped($_[0], @_ ? $_[0] : 0)  if @_ <= 1;

  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::vecsum(@_))
    if $Math::Prime::Util::_GMPfunc{"vecsum"};
  my $sum = 0;
  my $neglim = -(INTMAX >> 1) - 1;
  foreach my $v (@_) {
    $sum += $v;
    if ($sum > (INTMAX-250) || $sum < $neglim) {
      $sum = BZERO->copy;
      $sum->badd("$_") for @_;
      return $sum;
    }
  }
  $sum;
}

sub vecprod {
  return 1 unless @_;
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::vecprod(@_))
    if $Math::Prime::Util::_GMPfunc{"vecprod"};
  # Product tree:
  my $prod = _product(0, $#_, [map { Math::BigInt->new("$_") } @_]);
  # Linear:
  # my $prod = BONE->copy;  $prod *= "$_" for @_;
  $prod = _bigint_to_int($prod) if $prod->bacmp(BMAX) <= 0 && $prod->bcmp(-(BMAX>>1)) > 0;
  $prod;
}

sub vecmin {
  return unless @_;
  my $min = shift;
  for (@_) { $min = $_ if $_ < $min; }
  $min;
}
sub vecmax {
  return unless @_;
  my $max = shift;
  for (@_) { $max = $_ if $_ > $max; }
  $max;
}

sub vecextract {
  my($aref, $mask) = @_;

  return @$aref[@$mask] if ref($mask) eq 'ARRAY';

  # This is concise but very slow.
  # map { $aref->[$_] }  grep { $mask & (1 << $_) }  0 .. $#$aref;

  my($i, @v) = (0);
  while ($mask) {
    push @v, $i if $mask & 1;
    $mask >>= 1;
    $i++;
  }
  @$aref[@v];
}

sub sumdigits {
  my($n,$base) = @_;
  my $sum = 0;
  $base =  2 if !defined $base && $n =~ s/^0b//;
  $base = 16 if !defined $base && $n =~ s/^0x//;
  if (!defined $base || $base == 10) {
    $n =~ tr/0123456789//cd;
    $sum += $_ for (split(//,$n));
  } else {
    croak "sumdigits: invalid base $base" if $base < 2;
    my $cmap = substr("0123456789abcdefghijklmnopqrstuvwxyz",0,$base);
    for my $c (split(//,lc($n))) {
      my $p = index($cmap,$c);
      $sum += $p if $p > 0;
    }
  }
  $sum;
}

sub invmod {
  my($a,$n) = @_;
  return if $n == 0 || $a == 0;
  return 0 if $n == 1;
  $n = -$n if $n < 0;  # Pari semantics
  if ($n > ~0) {
    my $invmod = Math::BigInt->new("$a")->bmodinv("$n");
    return if !defined $invmod || $invmod->is_nan;
    $invmod = _bigint_to_int($invmod) if $invmod->bacmp(BMAX) <= 0;
    return $invmod;
  }
  my($t,$nt,$r,$nr) = (0, 1, $n, $a % $n);
  while ($nr != 0) {
    # Use mod before divide to force correct behavior with high bit set
    my $quot = int( ($r-($r % $nr))/$nr );
    ($nt,$t) = ($t-$quot*$nt,$nt);
    ($nr,$r) = ($r-$quot*$nr,$nr);
  }
  return if $r > 1;
  $t += $n if $t < 0;
  $t;
}

sub _verify_sqrtmod {
  my($r,$a,$n) = @_;
  if (ref($r)) {
    return if $r->copy->bmul($r)->bmod($n)->bcmp($a);
    $r = _bigint_to_int($r) if $r->bacmp(BMAX) <= 0;
  } else {
    return unless (($r*$r) % $n) == $a;
  }
  $r = $n-$r if $n-$r < $r;
  $r;
}

sub sqrtmod {
  my($a,$n) = @_;
  return if $n == 0;
  if ($n <= 2 || $a <= 1) {
    $a %= $n;
    return ((($a*$a) % $n) == $a) ? $a : undef;
  }

  if ($n < 10000000) {
    # Horrible trial search
    $a = _bigint_to_int($a);
    $n = _bigint_to_int($n);
    $a %= $n;
    return 1 if $a == 1;
    my $lim = ($n+1) >> 1;
    for my $r (2 .. $lim) {
      return $r if (($r*$r) % $n) == $a;
    }
    undef;
  }

  $a = Math::BigInt->new("$a") unless ref($a) eq 'Math::BigInt';
  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
  $a->bmod($n);
  my $r;

  if (($n % 4) == 3) {
    $r = $a->copy->bmodpow(($n+1)>>2, $n);
    return _verify_sqrtmod($r, $a, $n);
  }
  if (($n % 8) == 5) {
    my $q = $a->copy->bmodpow(($n-1)>>2, $n);
    if ($q->is_one) {
      $r = $a->copy->bmodpow(($n+3)>>3, $n);
    } else {
      my $v = $a->copy->bmul(4)->bmodpow(($n-5)>>3, $n);
      $r = $a->copy->bmul(2)->bmul($v)->bmod($n);
    }
    return _verify_sqrtmod($r, $a, $n);
  }

  return if $n->is_odd && !$a->copy->bmodpow(($n-1)>>1,$n)->is_one();

  # Horrible trial search.  Need to use Tonelli-Shanks here.
  $r = Math::BigInt->new(2);
  my $lim = int( ($n+1) / 2 );
  while ($r < $lim) {
    return $r if $r->copy->bmul($r)->bmod($n) == $a;
    $r++;
  }
  undef;
}

sub addmod {
  my($a, $b, $n) = @_;
  return 0 if $n <= 1;
  return _addmod($a,$b,$n) if $n < INTMAX && $a>=0 && $a=0 && $b
  my $ret = Math::BigInt->new("$a")->badd("$b")->bmod("$n");
  $ret = _bigint_to_int($ret) if $ret->bacmp(BMAX) <= 0;
  $ret;
}

sub mulmod {
  my($a, $b, $n) = @_;
  return 0 if $n <= 1;
  return _mulmod($a,$b,$n) if $n < INTMAX && $a>0 && $a0 && $b
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::mulmod($a,$b,$n))
    if $Math::Prime::Util::_GMPfunc{"mulmod"};
  my $ret = Math::BigInt->new("$a")->bmod("$n")->bmul("$b")->bmod("$n");
  $ret = _bigint_to_int($ret) if $ret->bacmp(BMAX) <= 0;
  $ret;
}
sub divmod {
  my($a, $b, $n) = @_;
  return 0 if $n <= 1;
  my $ret = Math::BigInt->new("$b")->bmodinv("$n")->bmul("$a")->bmod("$n");
  if ($ret->is_nan) {
    $ret = undef;
  } else {
    $ret = _bigint_to_int($ret) if $ret->bacmp(BMAX) <= 0;
  }
  $ret;
}
sub powmod {
  my($a, $b, $n) = @_;
  return 0 if $n <= 1;
  if ($Math::Prime::Util::_GMPfunc{"powmod"}) {
    my $r = Math::Prime::Util::GMP::powmod($a,$b,$n);
    return (defined $r) ? Math::Prime::Util::_reftyped($_[0], $r) : undef;
  }
  my $ret = Math::BigInt->new("$a")->bmod("$n")->bmodpow("$b","$n");
  if ($ret->is_nan) {
    $ret = undef;
  } else {
    $ret = _bigint_to_int($ret) if $ret->bacmp(BMAX) <= 0;
  }
  $ret;
}

# no validation, x is allowed to be negative, y must be >= 0
sub _gcd_ui {
  my($x, $y) = @_;
  if ($y < $x) { ($x, $y) = ($y, $x); }
  elsif ($x < 0) { $x = -$x; }
  while ($y > 0) {
    ($x, $y) = ($y, $x % $y);
  }
  $x;
}

sub is_power {
  my ($n, $a, $refp) = @_;
  croak("is_power third argument not a scalar reference") if defined($refp) && !ref($refp);
  _validate_integer($n);
  return 0 if abs($n) <= 3 && !$a;

  if ($Math::Prime::Util::_GMPfunc{"is_power"} &&
      ($Math::Prime::Util::GMP::VERSION >= 0.42 ||
       ($Math::Prime::Util::GMP::VERSION >= 0.28 && $n > 0))) {
    $a = 0 unless defined $a;
    my $k = Math::Prime::Util::GMP::is_power($n,$a);
    return 0 unless $k > 0;
    if (defined $refp) {
      $a = $k unless $a;
      my $isneg = ($n < 0);
      $n =~ s/^-// if $isneg;
      $$refp = Math::Prime::Util::rootint($n, $a);
      $$refp = Math::Prime::Util::_reftyped($_[0], $$refp) if $$refp > INTMAX;
      $$refp = -$$refp if $isneg;
    }
    return $k;
  }

  if (defined $a && $a != 0) {
    return 1 if $a == 1;                  # Everything is a 1st power
    return 0 if $n < 0 && $a % 2 == 0;    # Negative n never an even power
    if ($a == 2) {
      if (_is_perfect_square($n)) {
        $$refp = int(sqrt($n)) if defined $refp;
        return 1;
      }
    } else {
      $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
      my $root = $n->copy->babs->broot($a)->bfloor;
      $root->bneg if $n->is_neg;
      if ($root->copy->bpow($a) == $n) {
        $$refp = $root if defined $refp;
        return 1;
      }
    }
  } else {
    $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
    if ($n < 0) {
      my $absn = $n->copy->babs;
      my $root = is_power($absn, 0, $refp);
      return 0 unless $root;
      if ($root % 2 == 0) {
        my $power = valuation($root, 2);
        $root >>= $power;
        return 0 if $root == 1;
        $power = BTWO->copy->bpow($power);
        $$refp = $$refp ** $power if defined $refp;
      }
      $$refp = -$$refp if defined $refp;
      return $root;
    }
    my $e = 2;
    while (1) {
      my $root = $n->copy()->broot($e)->bfloor;
      last if $root->is_one();
      if ($root->copy->bpow($e) == $n) {
        my $next = is_power($root, 0, $refp);
        $$refp = $root if !$next && defined $refp;
        $e *= $next if $next != 0;
        return $e;
      }
      $e = next_prime($e);
    }
  }
  0;
}

sub is_square {
  my($n) = @_;
  return 0 if $n < 0;
  #is_power($n,2);
  _validate_integer($n);
  _is_perfect_square($n);
}

sub is_prime_power {
  my ($n, $refp) = @_;
  croak("is_prime_power second argument not a scalar reference") if defined($refp) && !ref($refp);
  return 0 if $n <= 1;

  if (Math::Prime::Util::is_prime($n)) { $$refp = $n if defined $refp; return 1; }
  my $r;
  my $k = Math::Prime::Util::is_power($n,0,\$r);
  if ($k) {
    $r = _bigint_to_int($r) if ref($r) && $r->bacmp(BMAX) <= 0;
    return 0 unless Math::Prime::Util::is_prime($r);
    $$refp = $r if defined $refp;
  }
  $k;
}

sub is_polygonal {
  my ($n, $k, $refp) = @_;
  croak("is_polygonal third argument not a scalar reference") if defined($refp) && !ref($refp);
  croak("is_polygonal: k must be >= 3") if $k < 3;
  return 0 if $n <= 0;
  if ($n == 1) { $$refp = 1 if defined $refp; return 1; }
  my($D,$R);
  if ($k == 4) {
    return 0 unless _is_perfect_square($n);
    $$refp = sqrtint($n) if defined $refp;
    return 1;
  }
  if ($n <= MPU_HALFWORD && $k <= MPU_HALFWORD) {
    $D = ($k==3) ? 1+($n<<3) : (8*$k-16)*$n + ($k-4)*($k-4);
    return 0 unless _is_perfect_square($D);
    $D = $k-4 + Math::Prime::Util::sqrtint($D);
    $R = 2*$k-4;
  } else {
    if ($k == 3) {
      $D = vecsum(1, vecprod($n, 8));
    } else {
      $D = vecsum(vecprod($n, vecprod(8, $k) - 16),  vecprod($k-4,$k-4));;
    }
    return 0 unless _is_perfect_square($D);
    $D = vecsum( sqrtint($D), $k-4 );
    $R = vecprod(2, $k) - 4;
  }
  return 0 if ($D % $R) != 0;
  $$refp = $D / $R if defined $refp;
  1;
}

sub valuation {
  my($n, $k) = @_;
  return 0 if $n < 2 || $k < 2;
  my $v = 0;
  if ($k == 2) { # Accelerate power of 2
    if (ref($n) eq 'Math::BigInt') {   # This can pay off for big inputs
      return 0 unless $n->is_even;
      my $s = $n->as_bin;              # We could do same for k=10
      return length($s) - rindex($s,'1') - 1;
    }
    while (!($n & 0xFFFF) ) {  $n >>=16;  $v +=16;  }
    while (!($n & 0x000F) ) {  $n >>= 4;  $v += 4;  }
  }
  while ( !($n % $k) ) {
    $n /= $k;
    $v++;
  }
  $v;
}

sub hammingweight {
  my $n = shift;
  return 0 + (Math::BigInt->new("$n")->as_bin() =~ tr/1//);
}

my @_digitmap = (0..9, 'a'..'z');
my %_mapdigit = map { $_digitmap[$_] => $_ } 0 .. $#_digitmap;
sub _splitdigits {
  my($n, $base, $len) = @_;    # n is num or bigint, base is in range
  my @d;
  if ($base == 10) {
    @d = split(//,"$n");
  } elsif ($base == 2) {
    @d = split(//,substr(Math::BigInt->new("$n")->as_bin,2));
  } elsif ($base == 16) {
    @d = map { $_mapdigit{$_} } split(//,substr(Math::BigInt->new("$n")->as_hex,2));
  } else {
    while ($n >= 1) {
      my $rem = $n % $base;
      unshift @d, $rem;
      $n = ($n-$rem)/$base;    # Always an exact division
    }
  }
  if ($len >= 0 && $len != scalar(@d)) {
    while (@d < $len) { unshift @d, 0; }
    while (@d > $len) { shift @d; }
  }
  @d;
}

sub todigits {
  my($n,$base,$len) = @_;
  $base = 10 unless defined $base;
  $len = -1 unless defined $len;
  die "Invalid base: $base" if $base < 2;
  return if $n == 0;
  $n = -$n if $n < 0;
  _validate_num($n) || _validate_positive_integer($n);
  _splitdigits($n, $base, $len);
}

sub todigitstring {
  my($n,$base,$len) = @_;
  $base = 10 unless defined $base;
  $len = -1 unless defined $len;
  $n =~ s/^-//;
  return substr(Math::BigInt->new("$n")->as_bin,2) if $base ==  2 && $len < 0;
  return substr(Math::BigInt->new("$n")->as_oct,1) if $base ==  8 && $len < 0;
  return substr(Math::BigInt->new("$n")->as_hex,2) if $base == 16 && $len < 0;
  my @d = ($n == 0) ? () : _splitdigits($n, $base, $len);
  return join("", @d) if $base <= 10;
  die "Invalid base for string: $base" if $base > 36;
  join("", map { $_digitmap[$_] } @d);
}

sub fromdigits {
  my($r, $base) = @_;
  $base = 10 unless defined $base;
  return $r if $base == 10 && ref($r) =~ /^Math::/;
  my $n;
  if (ref($r) && ref($r) !~ /^Math::/) {
    croak "fromdigits first argument must be a string or array reference"
      unless ref($r) eq 'ARRAY';
    ($n,$base) = (BZERO->copy, BZERO + $base);
    for my $d (@$r) {
      $n = $n * $base + $d;
    }
  } elsif ($base == 2) {
    $n = Math::BigInt->from_bin("0b$r");
  } elsif ($base == 8) {
    $n = Math::BigInt->from_oct("0$r");
  } elsif ($base == 16) {
    $n = Math::BigInt->from_hex("0x$r");
  } else {
    $r =~ s/^0*//;
    ($n,$base) = (BZERO->copy, BZERO + $base);
    #for my $d (map { $_mapdigit{$_} } split(//,$r)) {
    #  croak "Invalid digit for base $base" unless defined $d && $d < $base;
    #  $n = $n * $base + $d;
    #}
    for my $c (split(//, lc($r))) {
      $n->bmul($base);
      if ($c ne '0') {
        my $d = index("0123456789abcdefghijklmnopqrstuvwxyz", $c);
        croak "Invalid digit for base $base" unless $d >= 0;
        $n->badd($d);
      }
    }
  }
  $n = _bigint_to_int($n) if $n->bacmp(BMAX) <= 0;
  $n;
}

sub sqrtint {
  my($n) = @_;
  my $sqrt = Math::BigInt->new("$n")->bsqrt;
  return Math::Prime::Util::_reftyped($_[0], "$sqrt");
}

sub rootint {
  my ($n, $k, $refp) = @_;
  croak "rootint: k must be > 0" unless $k > 0;
  # Math::BigInt returns NaN for any root of a negative n.
  my $root = Math::BigInt->new("$n")->babs->broot("$k");
  if (defined $refp) {
    croak("logint third argument not a scalar reference") unless ref($refp);
    $$refp = $root->copy->bpow($k);
  }
  return Math::Prime::Util::_reftyped($_[0], "$root");
}

sub logint {
  my ($n, $b, $refp) = @_;
  croak("logint third argument not a scalar reference") if defined($refp) && !ref($refp);
  croak "logint: n must be > 0" unless $n > 0;
  croak "logint: missing base" unless defined $b;
  if ($b == 10) {
    my $e = length($n)-1;
    $$refp = Math::BigInt->new("1" . "0"x$e) if defined $refp;
    return $e;
  }
  if ($b == 2) {
    my $e = length(Math::BigInt->new("$n")->as_bin)-2-1;
    $$refp = Math::BigInt->from_bin("1" . "0"x$e) if defined $refp;
    return $e;
  }
  croak "logint: base must be > 1" unless $b > 1;

  my $e = Math::BigInt->new("$n")->blog("$b");
  $$refp = Math::BigInt->new("$b")->bpow($e) if defined $refp;
  return Math::Prime::Util::_reftyped($_[0], "$e");
}

# Seidel (Luschny), core using Trizen's simplications from Math::BigNum.
# http://oeis.org/wiki/User:Peter_Luschny/ComputationAndAsymptoticsOfBernoulliNumbers#Bernoulli_numbers__after_Seidel
sub _bernoulli_seidel {
  my($n) = @_;
  return (1,1) if $n == 0;
  return (0,1) if $n > 1 && $n % 2;

  my $oacc = Math::BigInt->accuracy();  Math::BigInt->accuracy(undef);
  my @D = (BZERO->copy, BONE->copy, map { BZERO->copy } 1 .. ($n>>1)-1);
  my ($h, $w) = (1, 1);

  foreach my $i (0 .. $n-1) {
    if ($w ^= 1) {
      $D[$_]->badd($D[$_-1]) for 1 .. $h-1;
    } else {
      $w = $h++;
      $D[$w]->badd($D[$w+1]) while --$w;
    }
  }
  my $num = $D[$h-1];
  my $den = BONE->copy->blsft($n+1)->bsub(BTWO);
  my $gcd = Math::BigInt::bgcd($num, $den);
  $num /= $gcd;
  $den /= $gcd;
  $num->bneg() if ($n % 4) == 0;
  Math::BigInt->accuracy($oacc);
  ($num,$den);
}

sub bernfrac {
  my $n = shift;
  return (BONE,BONE) if $n == 0;
  return (BONE,BTWO) if $n == 1;    # We're choosing 1/2 instead of -1/2
  return (BZERO,BONE) if $n < 0 || $n & 1;

  # We should have used one of the GMP functions.  At this point we could
  # replicate that with Math::MPFR, but the chance that they have the latter
  # but not the former is very small.

  _bernoulli_seidel($n);
}

sub stirling {
  my($n, $m, $type) = @_;
  return 1 if $m == $n;
  return 0 if $n == 0 || $m == 0 || $m > $n;
  $type = 1 unless defined $type;
  croak "stirling type must be 1, 2, or 3" unless $type == 1 || $type == 2 || $type == 3;
  if ($m == 1) {
    return 1 if $type == 2;
    return factorial($n) if $type == 3;
    return factorial($n-1) if $n&1;
    return vecprod(-1, factorial($n-1));
  }
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::stirling($n,$m,$type))
    if $Math::Prime::Util::_GMPfunc{"stirling"};
  my $s = BZERO->copy;
  if ($type == 3) {
    $s = Math::Prime::Util::vecprod( Math::Prime::Util::binomial($n,$m), Math::Prime::Util::binomial($n-1,$m-1), Math::Prime::Util::factorial($n-$m) );
  } elsif ($type == 2) {
    for my $j (1 .. $m) {
      my $t = Math::Prime::Util::vecprod(
                Math::BigInt->new($j) ** $n,
                Math::Prime::Util::binomial($m,$j)
              );
      $s = (($m-$j) & 1)  ?  $s - $t  :  $s + $t;
    }
    $s /= factorial($m);
  } else {
    for my $k (1 .. $n-$m) {
      my $t = Math::Prime::Util::vecprod(
        Math::Prime::Util::binomial($k + $n - 1, $k + $n - $m),
        Math::Prime::Util::binomial(2 * $n - $m, $n - $k - $m),
        Math::Prime::Util::stirling($k - $m + $n, $k, 2),
      );
      $s = ($k & 1)  ?  $s - $t  :  $s + $t;
    }
  }
  $s;
}

sub _harmonic_split { # From Fredrik Johansson
  my($a,$b) = @_;
  return (BONE, $a) if $b - $a == BONE;
  return ($a+$a+BONE, $a*$a+$a) if $b - $a == BTWO;   # Cut down recursion
  my $m = $a->copy->badd($b)->brsft(BONE);
  my ($p,$q) = _harmonic_split($a, $m);
  my ($r,$s) = _harmonic_split($m, $b);
  ($p*$s+$q*$r, $q*$s);
}

sub harmfrac {
  my($n) = @_;
  return (BZERO,BONE) if $n <= 0;
  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
  my($p,$q) = _harmonic_split($n-$n+1, $n+1);
  my $gcd = Math::BigInt::bgcd($p,$q);
  ( scalar $p->bdiv($gcd), scalar $q->bdiv($gcd) );
}

sub harmreal {
  my($n, $precision) = @_;

  do { require Math::BigFloat; Math::BigFloat->import(); } unless defined $Math::BigFloat::VERSION;
  return Math::BigFloat->bzero if $n <= 0;

  if (_MPFR_available(3,0)) {
    $precision = _find_big_acc($n) unless defined $precision;
    my $rnd = 0;  # MPFR_RNDN;
    my $bit_precision = int("$precision" * 3.322) + 7;
    my($n_mpfr, $euler, $psi) = map { Math::MPFR->new() } 1..3;
    Math::MPFR::Rmpfr_set_str($n_mpfr, "$n", 10, $rnd);
    Math::MPFR::Rmpfr_set_prec($euler, $bit_precision);
    Math::MPFR::Rmpfr_set_prec($psi, $bit_precision);
    Math::MPFR::Rmpfr_const_euler($euler, $rnd);
    Math::MPFR::Rmpfr_digamma($psi, $n_mpfr+1, $rnd);
    Math::MPFR::Rmpfr_add($psi, $psi, $euler, $rnd);
    my $strval = Math::MPFR::Rmpfr_get_str($psi, 10, 0, $rnd);
    return Math::BigFloat->new($strval,$precision);
  }

  # Use asymptotic formula for larger $n if possible.  Saves lots of time if
  # the default Calc backend is being used.
  {
    my $sprec = $precision;
    $sprec = Math::BigFloat->precision unless defined $sprec;
    $sprec = 40 unless defined $sprec;
    if ( ($sprec <= 23 && $n >    54) ||
         ($sprec <= 30 && $n >   348) ||
         ($sprec <= 40 && $n >  2002) ||
         ($sprec <= 50 && $n > 12644) ) {
      $n = Math::BigFloat->new($n, $sprec+5);
      my($n2, $one, $h) = ($n*$n, Math::BigFloat->bone, Math::BigFloat->bzero);
      my $nt = $n2;
      my $eps = Math::BigFloat->new(10)->bpow(-$sprec-4);
      foreach my $d (-12, 120, -252, 240, -132, 32760, -12, 8160, -14364, 6600, -276, 65520, -12) { # OEIS A006593
        my $term = $one/($d * $nt);
        last if $term->bacmp($eps) < 0;
        $h += $term;
        $nt *= $n2;
      }
      $h->badd(scalar $one->copy->bdiv(2*$n));
      $h->badd('0.57721566490153286060651209008240243104215933593992359880576723488486772677766467');
      $h->badd($n->copy->blog);
      $h->round($sprec);
      return $h;
    }
  }

  my($num,$den) = Math::Prime::Util::harmfrac($n);
  # Note, with Calc backend this can be very, very slow
  scalar Math::BigFloat->new($num)->bdiv($den, $precision);
}

sub is_pseudoprime {
  my($n, @bases) = @_;
  return 0 if int($n) < 0;
  _validate_positive_integer($n);
  croak("No bases given to is_pseudoprime") unless scalar(@bases) > 0;
  return 0+($n >= 2) if $n < 4;

  foreach my $base (@bases) {
    croak "Base $base is invalid" if $base < 2;
    $base = $base % $n if $base >= $n;
    if ($base > 1 && $base != $n-1) {
      my $x = (ref($n) eq 'Math::BigInt')
        ? $n->copy->bzero->badd($base)->bmodpow($n-1,$n)->is_one
        : _powmod($base, $n-1, $n);
      return 0 unless $x == 1;
    }
  }
  1;
}

sub is_euler_pseudoprime {
  my($n, @bases) = @_;
  return 0 if int($n) < 0;
  _validate_positive_integer($n);
  croak("No bases given to is_euler_pseudoprime") unless scalar(@bases) > 0;
  return 0+($n >= 2) if $n < 4;

  foreach my $base (@bases) {
    croak "Base $base is invalid" if $base < 2;
    $base = $base % $n if $base >= $n;
    if ($base > 1 && $base != $n-1) {
      my $j = kronecker($base, $n);
      return 0 if $j == 0;
      $j = ($j > 0) ? 1 : $n-1;
      my $x = (ref($n) eq 'Math::BigInt')
        ? $n->copy->bzero->badd($base)->bmodpow(($n-1)/2,$n)
        : _powmod($base, ($n-1)>>1, $n);
      return 0 unless $x == $j;
    }
  }
  1;
}

sub is_euler_plumb_pseudoprime {
  my($n) = @_;
  return 0 if int($n) < 0;
  _validate_positive_integer($n);
  return 0+($n >= 2) if $n < 4;
  return 0 if ($n % 2) == 0;
  my $nmod8 = $n % 8;
  my $exp = 1 + ($nmod8 == 1);
  my $ap = Math::Prime::Util::powmod(2, ($n-1) >> $exp, $n);
  if ($ap ==    1) { return ($nmod8 == 1 || $nmod8 == 7); }
  if ($ap == $n-1) { return ($nmod8 == 1 || $nmod8 == 3 || $nmod8 == 5); }
  0;
}

sub _miller_rabin_2 {
  my($n, $nm1, $s, $d) = @_;

  if ( ref($n) eq 'Math::BigInt' ) {

    if (!defined $nm1) {
      $nm1 = $n->copy->bdec();
      $s = 0;
      $d = $nm1->copy;
      do {
        $s++;
        $d->brsft(BONE);
      } while $d->is_even;
    }
    my $x = BTWO->copy->bmodpow($d,$n);
    return 1 if $x->is_one || $x->bcmp($nm1) == 0;
    foreach my $r (1 .. $s-1) {
      $x->bmul($x)->bmod($n);
      last if $x->is_one;
      return 1 if $x->bcmp($nm1) == 0;
    }

  } else {

    if (!defined $nm1) {
      $nm1 = $n-1;
      $s = 0;
      $d = $nm1;
      while ( ($d & 1) == 0 ) {
        $s++;
        $d >>= 1;
      }
    }

    if ($n < MPU_HALFWORD) {
      my $x = _native_powmod(2, $d, $n);
      return 1 if $x == 1 || $x == $nm1;
      foreach my $r (1 .. $s-1) {
        $x = ($x*$x) % $n;
        last if $x == 1;
        return 1 if $x == $n-1;
      }
    } else {
      my $x = _powmod(2, $d, $n);
      return 1 if $x == 1 || $x == $nm1;
      foreach my $r (1 .. $s-1) {
        $x = ($x < MPU_HALFWORD) ? ($x*$x) % $n : _mulmod($x, $x, $n);
        last if $x == 1;
        return 1 if $x == $n-1;
      }
    }
  }
  0;
}

sub is_strong_pseudoprime {
  my($n, @bases) = @_;
  return 0 if int($n) < 0;
  _validate_positive_integer($n);
  croak("No bases given to is_strong_pseudoprime") unless scalar(@bases) > 0;

  return 0+($n >= 2) if $n < 4;
  return 0 if ($n % 2) == 0;

  if ($bases[0] == 2) {
    return 0 unless _miller_rabin_2($n);
    shift @bases;
    return 1 unless @bases;
  }

  my @newbases;
  for my $base (@bases) {
    croak "Base $base is invalid" if $base < 2;
    $base %= $n if $base >= $n;
    return 0 if $base == 0 || ($base == $n-1 && ($base % 2) == 1);
    push @newbases, $base;
  }
  @bases = @newbases;

  if ( ref($n) eq 'Math::BigInt' ) {

    my $nminus1 = $n->copy->bdec();
    my $s = 0;
    my $d = $nminus1->copy;
    do {  # n is > 3 and odd, so n-1 must be even
      $s++;
      $d->brsft(BONE);
    } while $d->is_even;
    # Different way of doing the above.  Fewer function calls, slower on ave.
    #my $dbin = $nminus1->as_bin;
    #my $last1 = rindex($dbin, '1');
    #my $s = length($dbin)-2-$last1+1;
    #my $d = $nminus1->copy->brsft($s);

    foreach my $ma (@bases) {
      my $x = $n->copy->bzero->badd($ma)->bmodpow($d,$n);
      next if $x->is_one || $x->bcmp($nminus1) == 0;
      foreach my $r (1 .. $s-1) {
        $x->bmul($x); $x->bmod($n);
        return 0 if $x->is_one;
        do { $ma = 0; last; } if $x->bcmp($nminus1) == 0;
      }
      return 0 if $ma != 0;
    }

  } else {

   my $s = 0;
   my $d = $n - 1;
   while ( ($d & 1) == 0 ) {
     $s++;
     $d >>= 1;
   }

   if ($n < MPU_HALFWORD) {
    foreach my $ma (@bases) {
      my $x = _native_powmod($ma, $d, $n);
      next if ($x == 1) || ($x == ($n-1));
      foreach my $r (1 .. $s-1) {
        $x = ($x*$x) % $n;
        return 0 if $x == 1;
        last if $x == $n-1;
      }
      return 0 if $x != $n-1;
    }
   } else {
    foreach my $ma (@bases) {
      my $x = _powmod($ma, $d, $n);
      next if ($x == 1) || ($x == ($n-1));

      foreach my $r (1 .. $s-1) {
        $x = ($x < MPU_HALFWORD) ? ($x*$x) % $n : _mulmod($x, $x, $n);
        return 0 if $x == 1;
        last if $x == $n-1;
      }
      return 0 if $x != $n-1;
    }
   }

  }
  1;
}


# Calculate Kronecker symbol (a|b).  Cohen Algorithm 1.4.10.
# Extension of the Jacobi symbol, itself an extension of the Legendre symbol.
sub kronecker {
  my($a, $b) = @_;
  return (abs($a) == 1) ? 1 : 0  if $b == 0;
  my $k = 1;
  if ($b % 2 == 0) {
    return 0 if $a % 2 == 0;
    my $v = 0;
    do { $v++; $b /= 2; } while $b % 2 == 0;
    $k = -$k if $v % 2 == 1 && ($a % 8 == 3 || $a % 8 == 5);
  }
  if ($b < 0) {
    $b = -$b;
    $k = -$k if $a < 0;
  }
  if ($a < 0) { $a = -$a; $k = -$k if $b % 4 == 3; }
  $b = _bigint_to_int($b) if ref($b) eq 'Math::BigInt' && $b <= BMAX;
  $a = _bigint_to_int($a) if ref($a) eq 'Math::BigInt' && $a <= BMAX;
  # Now:  b > 0, b odd, a >= 0
  while ($a != 0) {
    if ($a % 2 == 0) {
      my $v = 0;
      do { $v++; $a /= 2; } while $a % 2 == 0;
      $k = -$k if $v % 2 == 1 && ($b % 8 == 3 || $b % 8 == 5);
    }
    $k = -$k if $a % 4 == 3 && $b % 4 == 3;
    ($a, $b) = ($b % $a, $a);
    # If a,b are bigints and now small enough, finish as native.
    if (   ref($a) eq 'Math::BigInt' && $a <= BMAX
        && ref($b) eq 'Math::BigInt' && $b <= BMAX) {
      return $k * kronecker(_bigint_to_int($a),_bigint_to_int($b));
    }
  }
  return ($b == 1) ? $k : 0;
}

sub _binomialu {
  my($r, $n, $k) = (1, @_);
  return ($k == $n) ? 1 : 0 if $k >= $n;
  $k = $n - $k if $k > ($n >> 1);
  foreach my $d (1 .. $k) {
    if ($r >= int(~0/$n)) {
      my($g, $nr, $dr);
      $g = _gcd_ui($n, $d);   $nr = int($n/$g);   $dr = int($d/$g);
      $g = _gcd_ui($r, $dr);  $r  = int($r/$g);   $dr = int($dr/$g);
      return 0 if $r >= int(~0/$nr);
      $r *= $nr;
      $r = int($r/$dr);
    } else {
      $r *= $n;
      $r = int($r/$d);
    }
    $n--;
  }
  $r;
}

sub binomial {
  my($n, $k) = @_;

  # 1. Try GMP
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::binomial($n,$k))
    if $Math::Prime::Util::_GMPfunc{"binomial"};

  # 2. Exit early for known 0 cases, and adjust k to be positive.
  if ($n >= 0) {  return 0 if $k < 0 || $k > $n;  }
  else         {  return 0 if $k < 0 && $k > $n;  }
  $k = $n - $k if $k < 0;

  # 3. Try to do in integer Perl
  my $r;
  if ($n >= 0) {
    $r = _binomialu($n, $k);
    return $r  if $r > 0;
  } else {
    $r = _binomialu(-$n+$k-1, $k);
    return $r   if $r > 0 && !($k & 1);
    return -$r  if $r > 0 && $r <= (~0>>1);
  }

  # 4. Overflow.  Solve using Math::BigInt
  return 1 if $k == 0;        # Work around bug in old
  return $n if $k == $n-1;    # Math::BigInt (fixed in 1.90)
  if ($n >= 0) {
    $r = Math::BigInt->new(''.$n)->bnok($k);
    $r = _bigint_to_int($r) if $r->bacmp(BMAX) <= 0;
  } else { # Math::BigInt is incorrect for negative n
    $r = Math::BigInt->new(''.(-$n+$k-1))->bnok($k);
    if ($k & 1) {
      $r->bneg;
      $r = _bigint_to_int($r) if $r->bacmp(''.(~0>>1)) <= 0;
    } else {
      $r = _bigint_to_int($r) if $r->bacmp(BMAX) <= 0;
    }
  }
  $r;
}

sub _product {
  my($a, $b, $r) = @_;
  if ($b <= $a) {
    $r->[$a];
  } elsif ($b == $a+1) {
    $r->[$a] -> bmul( $r->[$b] );
  } elsif ($b == $a+2) {
    $r->[$a] -> bmul( $r->[$a+1] ) -> bmul( $r->[$a+2] );
  } else {
    my $c = $a + (($b-$a+1)>>1);
    _product($a, $c-1, $r);
    _product($c, $b, $r);
    $r->[$a] -> bmul( $r->[$c] );
  }
}

sub factorial {
  my($n) = @_;
  return (1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800,479001600)[$n] if $n <= 12;
  return Math::GMP::bfac($n) if ref($n) eq 'Math::GMP';
  do { my $r = Math::GMPz->new(); Math::GMPz::Rmpz_fac_ui($r,$n); return $r; }
    if ref($n) eq 'Math::GMPz';
  if (Math::BigInt->config()->{lib} !~ /GMP|Pari/) {
    # It's not a GMP or GMPz object, and we have a slow bigint library.
    my $r;
    if (defined $Math::GMPz::VERSION) {
      $r = Math::GMPz->new(); Math::GMPz::Rmpz_fac_ui($r,$n);
    } elsif (defined $Math::GMP::VERSION) {
      $r = Math::GMP::bfac($n);
    } elsif (defined &Math::Prime::Util::GMP::factorial && Math::Prime::Util::prime_get_config()->{'gmp'}) {
      $r = Math::Prime::Util::GMP::factorial($n);
    }
    return Math::Prime::Util::_reftyped($_[0], $r)    if defined $r;
  }
  my $r = Math::BigInt->new($n)->bfac();
  $r = _bigint_to_int($r) if $r->bacmp(BMAX) <= 0;
  $r;
}

sub factorialmod {
  my($n,$m) = @_;

  return Math::Prime::Util::GMP::factorialmod($n,$m)
    if $Math::Prime::Util::_GMPfunc{"factorialmod"};

  if ($n > 10) {
    my($s,$t,$e) = (1);
    Math::Prime::Util::forprimes( sub {
      ($t,$e) = ($n,0);
      while ($t > 0) {
        $t = int($t/$_);
        $e += $t;
      }
      $s = Math::Prime::Util::mulmod($s, Math::Prime::Util::powmod($_,$e,$m), $m);
    }, 2, $n >> 1);
    Math::Prime::Util::forprimes( sub {
      $s = Math::Prime::Util::mulmod($s, $_, $m);
    }, ($n >> 1)+1, $n);
    return $s;
  }

  return factorial($n) % $m;
}

sub _is_perfect_square {
  my($n) = @_;

  if (ref($n) eq 'Math::BigInt') {
    my $mc = _bigint_to_int($n & 31);
    if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
      my $sq = $n->copy->bsqrt->bfloor;
      $sq->bmul($sq);
      return 1 if $sq == $n;
    }
  } else {
    my $mc = $n & 31;
    if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
      my $sq = int(sqrt($n));
      return 1 if ($sq*$sq) == $n;
    }
  }
  0;
}

sub is_primitive_root {
  my($a, $n) = @_;
  $n = -$n if $n < 0;  # Ignore sign of n
  return ($n==1) ? 1 : 0 if $n <= 1;
  $a %= $n if $a < 0 || $a >= $n;

  return Math::Prime::Util::GMP::is_primitive_root($a,$n)
    if $Math::Prime::Util::_GMPfunc{"is_primitive_root"};

  if ($Math::Prime::Util::_GMPfunc{"znorder"} && $Math::Prime::Util::_GMPfunc{"totient"}) {
    my $order = Math::Prime::Util::GMP::znorder($a,$n);
    return 0 unless defined $order;
    my $totient = Math::Prime::Util::GMP::totient($n);
    return ($order eq $totient) ? 1 : 0;
  }

  return 0 if Math::Prime::Util::gcd($a, $n) != 1;
  my $s = Math::Prime::Util::euler_phi($n);
  return 0 if ($s % 2) == 0 && Math::Prime::Util::powmod($a, $s/2, $n) == 1;
  return 0 if ($s % 3) == 0 && Math::Prime::Util::powmod($a, $s/3, $n) == 1;
  return 0 if ($s % 5) == 0 && Math::Prime::Util::powmod($a, $s/5, $n) == 1;
  foreach my $f (Math::Prime::Util::factor_exp($s)) {
    my $fp = $f->[0];
    return 0 if $fp > 5 && Math::Prime::Util::powmod($a, $s/$fp, $n) == 1;
  }
  1;
}

sub znorder {
  my($a, $n) = @_;
  return if $n <= 0;
  return 1 if $n == 1;
  return if $a <= 0;
  return 1 if $a == 1;

  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::znorder($a,$n))
    if $Math::Prime::Util::_GMPfunc{"znorder"};

  # Sadly, Calc/FastCalc are horrendously slow for this function.
  return if Math::Prime::Util::gcd($a, $n) > 1;

  # The answer is one of the divisors of phi(n) and lambda(n).
  my $lambda = Math::Prime::Util::carmichael_lambda($n);
  $a = Math::BigInt->new("$a") unless ref($a) eq 'Math::BigInt';

  # This is easy and usually fast, but can bog down with too many divisors.
  if ($lambda <= 2**64) {
    foreach my $k (Math::Prime::Util::divisors($lambda)) {
      return $k if Math::Prime::Util::powmod($a,$k,$n) == 1;
    }
    return;
  }

  # Algorithm 1.7 from A. Das applied to Carmichael Lambda.
  $lambda = Math::BigInt->new("$lambda") unless ref($lambda) eq 'Math::BigInt';
  my $k = Math::BigInt->bone;
  foreach my $f (Math::Prime::Util::factor_exp($lambda)) {
    my($pi, $ei, $enum) = (Math::BigInt->new("$f->[0]"), $f->[1], 0);
    my $phidiv = $lambda / ($pi**$ei);
    my $b = Math::Prime::Util::powmod($a,$phidiv,$n);
    while ($b != 1) {
      return if $enum++ >= $ei;
      $b = Math::Prime::Util::powmod($b,$pi,$n);
      $k *= $pi;
    }
  }
  $k = _bigint_to_int($k) if $k->bacmp(BMAX) <= 0;
  return $k;
}

sub _dlp_trial {
  my ($a,$g,$p,$limit) = @_;
  $limit = $p if !defined $limit || $limit > $p;
  my $t = $g->copy;

  if ($limit < 1_000_000_000) {
    for my $k (1 .. $limit) {
      return $k if $t == $a;
      $t = Math::Prime::Util::mulmod($t, $g, $p);
    }
    return 0;
  }

  for (my $k = BONE->copy; $k < $limit; $k->binc) {
    if ($t == $a) {
      $k = _bigint_to_int($k) if $k->bacmp(BMAX) <= 0;
      return $k;
    }
    $t->bmul($g)->bmod($p);
  }
  0;
}
sub _dlp_bsgs {
  my ($a,$g,$p,$n,$_verbose) = @_;
  my $invg = invmod($g, $p);
  return unless defined $invg;
  my $maxm = Math::Prime::Util::sqrtint($n)+1;
  my $b = ($p + $maxm - 1) / $maxm;
  # Limit for time and space.
  $b = ($b > 4_000_000) ? 4_000_000 : int("$b");
  $maxm = ($maxm > $b) ? $b : int("$maxm");

  my %hash;
  my $am = BONE->copy;
  my $gm = Math::Prime::Util::powmod($invg, $maxm, $p);
  my $key = $a->copy;
  my $r;

  foreach my $m (0 .. $b) {
    # Baby Step
    if ($m <= $maxm) {
      $r = $hash{"$am"};
      if (defined $r) {
        print "  bsgs found in stage 1 after $m tries\n" if $_verbose;
        $r = Math::Prime::Util::addmod($m, Math::Prime::Util::mulmod($r,$maxm,$p), $p);
        return $r;
      }
      $hash{"$am"} = $m;
      $am = Math::Prime::Util::mulmod($am,$g,$p);
      if ($am == $a) {
        print "  bsgs found during bs\n" if $_verbose;
        return $m+1;
      }
    }

    # Giant Step
    $r = $hash{"$key"};
    if (defined $r) {
      print "  bsgs found in stage 2 after $m tries\n" if $_verbose;
      $r = Math::Prime::Util::addmod($r, Math::Prime::Util::mulmod($m,$maxm,$p), $p);
      return $r;
    }
    $hash{"$key"} = $m if $m <= $maxm;
    $key = Math::Prime::Util::mulmod($key,$gm,$p);
  }
  0;
}

sub znlog {
  my ($a,$g,$p) =
    map { ref($_) eq 'Math::BigInt' ? $_ : Math::BigInt->new("$_") } @_;
  $a->bmod($p);
  $g->bmod($p);
  return 0 if $a == 1 || $g == 0 || $p < 2;
  my $_verbose = Math::Prime::Util::prime_get_config()->{'verbose'};

  # For large p, znorder can be very slow.  Do trial test first.
  my $x = _dlp_trial($a, $g, $p, 200);
  if ($x == 0) {
    my $n = znorder($g, $p);
    if (defined $n && $n > 1000) {
      $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
      $x = _dlp_bsgs($a, $g, $p, $n, $_verbose);
      $x = _bigint_to_int($x) if ref($x) && $x->bacmp(BMAX) <= 0;
      return $x if $x > 0 && $g->copy->bmodpow($x, $p) == $a;
      print "  BSGS giving up\n" if $x == 0 && $_verbose;
      print "  BSGS incorrect answer $x\n" if $x > 0 && $_verbose > 1;
    }
    $x = _dlp_trial($a,$g,$p);
  }
  $x = _bigint_to_int($x) if ref($x) && $x->bacmp(BMAX) <= 0;
  return ($x == 0) ? undef : $x;
}

sub znprimroot {
  my($n) = @_;
  $n = -$n if $n < 0;
  if ($n <= 4) {
    return if $n == 0;
    return $n-1;
  }
  return if $n % 4 == 0;
  my $a = 1;
  my $phi = $n-1;
  if (!is_prob_prime($n)) {
    $phi = euler_phi($n);
    # Check that a primitive root exists.
    return if $phi != Math::Prime::Util::carmichael_lambda($n);
  }
  my @exp = map { Math::BigInt->new("$_") }
            map { int($phi/$_->[0]) }
            Math::Prime::Util::factor_exp($phi);
  #print "phi: $phi  factors: ", join(",",factor($phi)), "\n";
  #print "  exponents: ", join(",", @exp), "\n";
  while (1) {
    my $fail = 0;
    do { $a++ } while Math::Prime::Util::kronecker($a,$n) == 0;
    return if $a >= $n;
    foreach my $f (@exp) {
      if (Math::Prime::Util::powmod($a,$f,$n) == 1) {
        $fail = 1;
        last;
      }
    }
    return $a if !$fail;
  }
}


# Find first D in sequence (5,-7,9,-11,13,-15,...) where (D|N) == -1
sub _lucas_selfridge_params {
  my($n) = @_;

  # D is typically quite small: 67 max for N < 10^19.  However, it is
  # theoretically possible D could grow unreasonably.  I'm giving up at 4000M.
  my $d = 5;
  my $sign = 1;
  while (1) {
    my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($d, $n)
                                          : _gcd_ui($d, $n);
    return (0,0,0) if $gcd > 1 && $gcd != $n;  # Found divisor $d
    my $j = kronecker($d * $sign, $n);
    last if $j == -1;
    $d += 2;
    croak "Could not find Jacobi sequence for $n" if $d > 4_000_000_000;
    $sign = -$sign;
  }
  my $D = $sign * $d;
  my $P = 1;
  my $Q = int( (1 - $D) / 4 );
  ($P, $Q, $D)
}

sub _lucas_extrastrong_params {
  my($n, $increment) = @_;
  $increment = 1 unless defined $increment;

  my ($P, $Q, $D) = (3, 1, 5);
  while (1) {
    my $gcd = (ref($n) eq 'Math::BigInt') ? Math::BigInt::bgcd($D, $n)
                                          : _gcd_ui($D, $n);
    return (0,0,0) if $gcd > 1 && $gcd != $n;  # Found divisor $d
    last if kronecker($D, $n) == -1;
    $P += $increment;
    croak "Could not find Jacobi sequence for $n" if $P > 65535;
    $D = $P*$P - 4;
  }
  ($P, $Q, $D);
}

# returns U_k, V_k, Q_k all mod n
sub lucas_sequence {
  my($n, $P, $Q, $k) = @_;

  croak "lucas_sequence: n must be >= 2" if $n < 2;
  croak "lucas_sequence: k must be >= 0" if $k < 0;
  croak "lucas_sequence: P out of range" if abs($P) >= $n;
  croak "lucas_sequence: Q out of range" if abs($Q) >= $n;

  if ($Math::Prime::Util::_GMPfunc{"lucas_sequence"} && $Math::Prime::Util::GMP::VERSION >= 0.30) {
    return map { ($_ > ''.~0) ? Math::BigInt->new(''.$_) : $_ }
           Math::Prime::Util::GMP::lucas_sequence($n, $P, $Q, $k);
  }

  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';

  my $ZERO = $n->copy->bzero;
  $P = $ZERO+$P unless ref($P) eq 'Math::BigInt';
  $Q = $ZERO+$Q unless ref($Q) eq 'Math::BigInt';
  my $D = $P*$P - BTWO*BTWO*$Q;
  if ($D->is_zero) {
    my $S = ($ZERO+$P) >> 1;
    my $U = $S->copy->bmodpow($k-1,$n)->bmul($k)->bmod($n);
    my $V = $S->copy->bmodpow($k,$n)->bmul(BTWO)->bmod($n);
    my $Qk = ($ZERO+$Q)->bmodpow($k, $n);
    return ($U, $V, $Qk);
  }
  my $U = BONE->copy;
  my $V = $P->copy;
  my $Qk = $Q->copy;

  return (BZERO->copy, BTWO->copy, $Qk) if $k == 0;
  $k = Math::BigInt->new("$k") unless ref($k) eq 'Math::BigInt';
  my $kstr = substr($k->as_bin, 2);
  my $bpos = 0;

  if (($n % 2)==0) {
    $P->bmod($n);
    $Q->bmod($n);
    my($Uh,$Vl, $Vh, $Ql, $Qh) = (BONE->copy, BTWO->copy, $P->copy, BONE->copy, BONE->copy);
    my ($b,$s) = (length($kstr)-1, 0);
    if ($kstr =~ /(0+)$/) { $s = length($1); }
    for my $bpos (0 .. $b-$s-1) {
      $Ql->bmul($Qh)->bmod($n);
      if (substr($kstr,$bpos,1)) {
        $Qh = $Ql * $Q;
        $Uh->bmul($Vh)->bmod($n);
        $Vl->bmul($Vh)->bsub($P * $Ql)->bmod($n);
        $Vh->bmul($Vh)->bsub(BTWO * $Qh)->bmod($n);
      } else {
        $Qh = $Ql->copy;
        $Uh->bmul($Vl)->bsub($Ql)->bmod($n);
        $Vh->bmul($Vl)->bsub($P * $Ql)->bmod($n);
        $Vl->bmul($Vl)->bsub(BTWO * $Ql)->bmod($n);
      }
    }
    $Ql->bmul($Qh);
    $Qh = $Ql * $Q;
    $Uh->bmul($Vl)->bsub($Ql)->bmod($n);
    $Vl->bmul($Vh)->bsub($P * $Ql)->bmod($n);
    $Ql->bmul($Qh)->bmod($n);
    for (1 .. $s) {
      $Uh->bmul($Vl)->bmod($n);
      $Vl->bmul($Vl)->bsub(BTWO * $Ql)->bmod($n);
      $Ql->bmul($Ql)->bmod($n);
    }
    ($U, $V, $Qk) = ($Uh, $Vl, $Ql);
  } elsif ($Q->is_one) {
    my $Dinverse = $D->copy->bmodinv($n);
    if ($P > BTWO && !$Dinverse->is_nan) {
      # Calculate V_k with U=V_{k+1}
      $U = $P->copy->bmul($P)->bsub(BTWO)->bmod($n);
      while (++$bpos < length($kstr)) {
        if (substr($kstr,$bpos,1)) {
          $V->bmul($U)->bsub($P  )->bmod($n);
          $U->bmul($U)->bsub(BTWO)->bmod($n);
        } else {
          $U->bmul($V)->bsub($P  )->bmod($n);
          $V->bmul($V)->bsub(BTWO)->bmod($n);
        }
      }
      # Crandall and Pomerance eq 3.13: U_n = D^-1 (2V_{n+1} - PV_n)
      $U = $Dinverse * (BTWO*$U - $P*$V);
    } else {
      while (++$bpos < length($kstr)) {
        $U->bmul($V)->bmod($n);
        $V->bmul($V)->bsub(BTWO)->bmod($n);
        if (substr($kstr,$bpos,1)) {
          my $T1 = $U->copy->bmul($D);
          $U->bmul($P)->badd( $V);
          $U->badd($n) if $U->is_odd;
          $U->brsft(BONE);
          $V->bmul($P)->badd($T1);
          $V->badd($n) if $V->is_odd;
          $V->brsft(BONE);
        }
      }
    }
  } else {
    my $qsign = ($Q == -1) ? -1 : 0;
    while (++$bpos < length($kstr)) {
      $U->bmul($V)->bmod($n);
      if    ($qsign ==  1) { $V->bmul($V)->bsub(BTWO)->bmod($n); }
      elsif ($qsign == -1) { $V->bmul($V)->badd(BTWO)->bmod($n); }
      else { $V->bmul($V)->bsub($Qk->copy->blsft(BONE))->bmod($n); }
      if (substr($kstr,$bpos,1)) {
        my $T1 = $U->copy->bmul($D);
        $U->bmul($P)->badd( $V);
        $U->badd($n) if $U->is_odd;
        $U->brsft(BONE);

        $V->bmul($P)->badd($T1);
        $V->badd($n) if $V->is_odd;
        $V->brsft(BONE);

        if ($qsign != 0) { $qsign = -1; }
        else             { $Qk->bmul($Qk)->bmul($Q)->bmod($n); }
      } else {
        if ($qsign != 0) { $qsign = 1; }
        else             { $Qk->bmul($Qk)->bmod($n); }
      }
    }
    if    ($qsign ==  1) { $Qk->bneg; }
    elsif ($qsign == -1) { $Qk = $n->copy->bdec; }
  }
  $U->bmod($n);
  $V->bmod($n);
  return ($U, $V, $Qk);
}
sub _lucasuv {
  my($P, $Q, $k) = @_;

  croak "lucas_sequence: k must be >= 0" if $k < 0;
  return (0,2) if $k == 0;

  $P = Math::BigInt->new("$P") unless ref($P) eq 'Math::BigInt';
  $Q = Math::BigInt->new("$Q") unless ref($Q) eq 'Math::BigInt';

  # Simple way, very slow as k increases:
  #my($U0, $U1) = (BZERO->copy, BONE->copy);
  #my($V0, $V1) = (BTWO->copy, Math::BigInt->new("$P"));
  #for (2 .. $k) {
  #  ($U0,$U1) = ($U1, $P*$U1 - $Q*$U0);
  #  ($V0,$V1) = ($V1, $P*$V1 - $Q*$V0);
  #}
  #return ($U1, $V1);

  my($Uh,$Vl, $Vh, $Ql, $Qh) = (BONE->copy, BTWO->copy, $P->copy, BONE->copy, BONE->copy);
  $k = Math::BigInt->new("$k") unless ref($k) eq 'Math::BigInt';
  my $kstr = substr($k->as_bin, 2);
  my ($n,$s) = (length($kstr)-1, 0);
  if ($kstr =~ /(0+)$/) { $s = length($1); }

  if ($Q == -1) {
    # This could be simplified, and it's running 10x slower than it should.
    my ($ql,$qh) = (1,1);
    for my $bpos (0 .. $n-$s-1) {
      $ql *= $qh;
      if (substr($kstr,$bpos,1)) {
        $qh = -$ql;
        $Uh->bmul($Vh);
        if ($ql == 1) {
          $Vl->bmul($Vh)->bsub( $P );
          $Vh->bmul($Vh)->badd( BTWO );
        } else {
          $Vl->bmul($Vh)->badd( $P );
          $Vh->bmul($Vh)->bsub( BTWO );
        }
      } else {
        $qh = $ql;
        if ($ql == 1) {
          $Uh->bmul($Vl)->bdec;
          $Vh->bmul($Vl)->bsub($P);
          $Vl->bmul($Vl)->bsub(BTWO);
        } else {
          $Uh->bmul($Vl)->binc;
          $Vh->bmul($Vl)->badd($P);
          $Vl->bmul($Vl)->badd(BTWO);
        }
      }
    }
    $ql *= $qh;
    $qh = -$ql;
    if ($ql == 1) {
      $Uh->bmul($Vl)->bdec;
      $Vl->bmul($Vh)->bsub($P);
    } else {
      $Uh->bmul($Vl)->binc;
      $Vl->bmul($Vh)->badd($P);
    }
    $ql *= $qh;
    for (1 .. $s) {
      $Uh->bmul($Vl);
      if ($ql == 1) { $Vl->bmul($Vl)->bsub(BTWO); $ql *= $ql; }
      else          { $Vl->bmul($Vl)->badd(BTWO); $ql *= $ql; }
    }
    return map { ($_ > ''.~0) ? Math::BigInt->new(''.$_) : $_ } ($Uh, $Vl);
  }

  for my $bpos (0 .. $n-$s-1) {
    $Ql->bmul($Qh);
    if (substr($kstr,$bpos,1)) {
      $Qh = $Ql * $Q;
      #$Uh = $Uh * $Vh;
      #$Vl = $Vh * $Vl - $P * $Ql;
      #$Vh = $Vh * $Vh - BTWO * $Qh;
      $Uh->bmul($Vh);
      $Vl->bmul($Vh)->bsub($P * $Ql);
      $Vh->bmul($Vh)->bsub(BTWO * $Qh);
    } else {
      $Qh = $Ql->copy;
      #$Uh = $Uh * $Vl - $Ql;
      #$Vh = $Vh * $Vl - $P * $Ql;
      #$Vl = $Vl * $Vl - BTWO * $Ql;
      $Uh->bmul($Vl)->bsub($Ql);
      $Vh->bmul($Vl)->bsub($P * $Ql);
      $Vl->bmul($Vl)->bsub(BTWO * $Ql);
    }
  }
  $Ql->bmul($Qh);
  $Qh = $Ql * $Q;
  $Uh->bmul($Vl)->bsub($Ql);
  $Vl->bmul($Vh)->bsub($P * $Ql);
  $Ql->bmul($Qh);
  for (1 .. $s) {
    $Uh->bmul($Vl);
    $Vl->bmul($Vl)->bsub(BTWO * $Ql);
    $Ql->bmul($Ql);
  }
  return map { ($_ > ''.~0) ? Math::BigInt->new(''.$_) : $_ } ($Uh, $Vl, $Ql);
}
sub lucasu { (_lucasuv(@_))[0] }
sub lucasv { (_lucasuv(@_))[1] }

sub is_lucas_pseudoprime {
  my($n) = @_;

  return 0+($n >= 2) if $n < 4;
  return 0 if ($n % 2) == 0 || _is_perfect_square($n);

  my ($P, $Q, $D) = _lucas_selfridge_params($n);
  return 0 if $D == 0;  # We found a divisor in the sequence
  die "Lucas parameter error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);

  my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $n+1);
  return ($U == 0) ? 1 : 0;
}

sub is_strong_lucas_pseudoprime {
  my($n) = @_;

  return 0+($n >= 2) if $n < 4;
  return 0 if ($n % 2) == 0 || _is_perfect_square($n);

  my ($P, $Q, $D) = _lucas_selfridge_params($n);
  return 0 if $D == 0;  # We found a divisor in the sequence
  die "Lucas parameter error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);

  my $m = $n+1;
  my($s, $k) = (0, $m);
  while ( $k > 0 && !($k % 2) ) {
    $s++;
    $k >>= 1;
  }
  my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k);

  return 1 if $U == 0;
  $V = Math::BigInt->new("$V") unless ref($V) eq 'Math::BigInt';
  $Qk = Math::BigInt->new("$Qk") unless ref($Qk) eq 'Math::BigInt';
  foreach my $r (0 .. $s-1) {
    return 1 if $V->is_zero;
    if ($r < ($s-1)) {
      $V->bmul($V)->bsub(BTWO*$Qk)->bmod($n);
      $Qk->bmul($Qk)->bmod($n);
    }
  }
  return 0;
}

sub is_extra_strong_lucas_pseudoprime {
  my($n) = @_;

  return 0+($n >= 2) if $n < 4;
  return 0 if ($n % 2) == 0 || _is_perfect_square($n);

  my ($P, $Q, $D) = _lucas_extrastrong_params($n);
  return 0 if $D == 0;  # We found a divisor in the sequence
  die "Lucas parameter error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);

  # We have to convert n to a bigint or Math::BigInt::GMP's stupid set_si bug
  # (RT 71548) will hit us and make the test $V == $n-2 always return false.
  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';

  my($s, $k) = (0, $n->copy->binc);
  while ($k->is_even && !$k->is_zero) {
    $s++;
    $k->brsft(BONE);
  }

  my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $k);

  return 1 if $U == 0 && ($V == BTWO || $V == ($n - BTWO));
  $V = Math::BigInt->new("$V") unless ref($V) eq 'Math::BigInt';
  foreach my $r (0 .. $s-2) {
    return 1 if $V->is_zero;
    $V->bmul($V)->bsub(BTWO)->bmod($n);
  }
  return 0;
}

sub is_almost_extra_strong_lucas_pseudoprime {
  my($n, $increment) = @_;
  $increment = 1 unless defined $increment;

  return 0+($n >= 2) if $n < 4;
  return 0 if ($n % 2) == 0 || _is_perfect_square($n);

  my ($P, $Q, $D) = _lucas_extrastrong_params($n, $increment);
  return 0 if $D == 0;  # We found a divisor in the sequence
  die "Lucas parameter error: $D, $P, $Q\n" if ($D != $P*$P - 4*$Q);

  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';

  my $ZERO = $n->copy->bzero;
  my $TWO = $ZERO->copy->binc->binc;
  my $V = $ZERO + $P;           # V_{k}
  my $W = $ZERO + $P*$P-$TWO;   # V_{k+1}
  my $kstr = substr($n->copy->binc()->as_bin, 2);
  $kstr =~ s/(0*)$//;
  my $s = length($1);
  my $bpos = 0;
  while (++$bpos < length($kstr)) {
    if (substr($kstr,$bpos,1)) {
      $V->bmul($W)->bsub($P  )->bmod($n);
      $W->bmul($W)->bsub($TWO)->bmod($n);
    } else {
      $W->bmul($V)->bsub($P  )->bmod($n);
      $V->bmul($V)->bsub($TWO)->bmod($n);
    }
  }

  return 1 if $V == 2 || $V == ($n-$TWO);
  foreach my $r (0 .. $s-2) {
    return 1 if $V->is_zero;
    $V->bmul($V)->bsub($TWO)->bmod($n);
  }
  return 0;
}

sub is_frobenius_khashin_pseudoprime {
  my($n) = @_;
  return 0+($n >= 2) if $n < 4;
  return 0 unless $n % 2;
  return 0 if _is_perfect_square($n);

  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';

  my $k;
  my $c = 1;
  do {
    $c += 2;
    $k = kronecker($c, $n);
  } while $k == 1;
  return 0 if $k == 0;

  my($ra,$rb,$a,$b,$d) = (1,1,1,1,$n-1);
  while (!$d->is_zero) {
    if ($d->is_odd()) {
      ($ra, $rb) = ( (($ra*$a)%$n + ((($rb*$b)%$n)*$c)%$n) % $n,
                     (($rb*$a)%$n + ($ra*$b)%$n) % $n );
    }
    $d >>= 1;
    if (!$d->is_zero) {
      ($a, $b) = ( (($a*$a)%$n + ((($b*$b)%$n)*$c)%$n) % $n,
                   (($b*$a)%$n + ($a*$b)%$n) % $n );
    }
  }
  return ($ra == 1 && $rb == $n-1) ? 1 : 0;
}

sub is_frobenius_underwood_pseudoprime {
  my($n) = @_;
  return 0+($n >= 2) if $n < 4;
  return 0 unless $n % 2;

  my($a, $temp1, $temp2);
  if ($n % 4 == 3) {
    $a = 0;
  } else {
    for ($a = 1; $a < 1000000; $a++) {
      next if $a==2 || $a==4 || $a==7 || $a==8 || $a==10 || $a==14 || $a==16 || $a==18;
      my $j = kronecker($a*$a - 4, $n);
      last if $j == -1;
      return 0 if $j == 0 || ($a == 20 && _is_perfect_square($n));
    }
  }
  $temp1 = Math::Prime::Util::gcd(($a+4)*(2*$a+5), $n);
  return 0 if $temp1 != 1 && $temp1 != $n;

  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
  my $ZERO = $n->copy->bzero;
  my $ONE = $ZERO->copy->binc;
  my $TWO = $ONE->copy->binc;
  my($s, $t) = ($ONE->copy, $TWO->copy);

  my $ap2 = $TWO + $a;
  my $np1string = substr( $n->copy->binc->as_bin, 2);
  my $np1len = length($np1string);

  foreach my $bit (1 .. $np1len-1) {
    $temp2 = $t+$t;
    $temp2 += ($s * $a)  if $a != 0;
    $temp1 = $temp2 * $s;
    $temp2 = $t - $s;
    $s += $t;
    $t = ($s * $temp2) % $n;
    $s = $temp1 % $n;
    if ( substr( $np1string, $bit, 1 ) ) {
      if ($a == 0)  { $temp1 = $s + $s; }
      else          { $temp1 = $s * $ap2; }
      $temp1 += $t;
      $t->badd($t)->bsub($s);   # $t = ($t+$t) - $s;
      $s = $temp1;
    }
  }
  $temp1 = (2*$a+5) % $n;
  return ($s == 0 && $t == $temp1) ? 1 : 0;
}

sub _perrin_signature {
  my($n) = @_;
  my @S = (1,$n-1,3, 3,0,2);
  return @S if $n <= 1;

  my @nbin = todigits($n,2);
  shift @nbin;

  while (@nbin) {
    my @T = map { addmod(addmod(Math::Prime::Util::mulmod($S[$_],$S[$_],$n), $n-$S[5-$_],$n), $n-$S[5-$_],$n); } 0..5;
    my $T01 = addmod($T[2], $n-$T[1], $n);
    my $T34 = addmod($T[5], $n-$T[4], $n);
    my $T45 = addmod($T34, $T[3], $n);
    if (shift @nbin) {
      @S = ($T[0], $T01, $T[1], $T[4], $T45, $T[5]);
    } else {
      @S = ($T01, $T[1], addmod($T01,$T[0],$n), $T34, $T[4], $T45);
    }
  }
  @S;
}

sub is_perrin_pseudoprime {
  my($n, $restrict) = @_;
  $restrict = 0 unless defined $restrict;
  return 0+($n >= 2) if $n < 4;
  return 0 if $restrict > 2 && ($n % 2) == 0;

  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';

  my @S = _perrin_signature($n);
  return 0 unless $S[4] == 0;
  return 1 if $restrict == 0;
  return 0 unless $S[1] == $n-1;
  return 1 if $restrict == 1;
  my $j = kronecker(-23,$n);
  if ($j == -1) {
    my $B = $S[2];
    my $B2 = mulmod($B,$B,$n);
    my $A = addmod(addmod(1,mulmod(3,$B,$n),$n),$n-$B2,$n);
    my $C = addmod(mulmod(3,$B2,$n),$n-2,$n);
    return 1 if $S[0] == $A && $S[2] == $B && $S[3] == $B && $S[5] == $C && $B != 3 && addmod(mulmod($B2,$B,$n),$n-$B,$n) == 1;
  } else {
    return 0 if $j == 0 && $n != 23 && $restrict > 2;
    return 1 if $S[0] == 1 && $S[2] == 3 && $S[3] == 3 && $S[5] == 2;
    return 1 if $S[0] == 0 && $S[5] == $n-1 && $S[2] != $S[3] && addmod($S[2],$S[3],$n) == $n-3 && mulmod(addmod($S[2],$n-$S[3],$n),addmod($S[2],$n-$S[3],$n),$n) == $n-(23%$n);
  }
  0;
}

sub is_catalan_pseudoprime {
  my($n) = @_;
  return 0+($n >= 2) if $n < 4;
  my $m = ($n-1)>>1;
  return (binomial($m<<1,$m) % $n) == (($m&1) ? $n-1 : 1) ? 1 : 0;
}

sub is_frobenius_pseudoprime {
  my($n, $P, $Q) = @_;
  ($P,$Q) = (0,0) unless defined $P && defined $Q;
  return 0+($n >= 2) if $n < 4;

  $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
  return 0 if $n->is_even;

  my($k, $Vcomp, $D, $Du) = (0, 4);
  if ($P == 0 && $Q == 0) {
    ($P,$Q) = (-1,2);
    while ($k != -1) {
      $P += 2;
      $P = 5 if $P == 3;  # Skip 3
      $D = $P*$P-4*$Q;
      $Du = ($D >= 0) ? $D : -$D;
      last if $P >= $n || $Du >= $n;   # TODO: remove?
      $k = kronecker($D, $n);
      return 0 if $k == 0;
      return 0 if $P == 10001 && _is_perfect_square($n);
    }
  } else {
    $D = $P*$P-4*$Q;
    $Du = ($D >= 0) ? $D : -$D;
    croak "Frobenius invalid P,Q: ($P,$Q)" if _is_perfect_square($Du);
  }
  return (is_prime($n) ? 1 : 0) if $n <= $Du || $n <= abs($Q) || $n <= abs($P);
  return 0 if Math::Prime::Util::gcd(abs($P*$Q*$D), $n) > 1;

  if ($k == 0) {
    $k = kronecker($D, $n);
    return 0 if $k == 0;
    my $Q2 = (2*abs($Q)) % $n;
    $Vcomp = ($k == 1) ? 2 : ($Q >= 0) ? $Q2 : $n-$Q2;
  }

  my($U, $V, $Qk) = lucas_sequence($n, $P, $Q, $n-$k);
  return 1 if $U == 0 && $V == $Vcomp;
  0;
}

# Since people have graciously donated millions of CPU years to doing these
# tests, it would be rude of us not to use the results.  This means we don't
# actually use the pretest and Lucas-Lehmer test coded below for any reasonable
# size number.
# See: http://www.mersenne.org/report_milestones/
my %_mersenne_primes;
undef @_mersenne_primes{2,3,5,7,13,17,19,31,61,89,107,127,521,607,1279,2203,2281,3217,4253,4423,9689,9941,11213,19937,21701,23209,44497,86243,110503,132049,216091,756839,859433,1257787,1398269,2976221,3021377,6972593,13466917,20996011,24036583,25964951,30402457,32582657,37156667,42643801,43112609,57885161,74207281};

sub is_mersenne_prime {
  my $p = shift;

  # Use the known Mersenne primes
  return 1 if exists $_mersenne_primes{$p};
  return 0 if $p < 34007399; # GIMPS has checked all below
  # Past this we do a generic Mersenne prime test

  return 1 if $p == 2;
  return 0 unless is_prob_prime($p);
  return 0 if $p > 3 && $p % 4 == 3 && $p < ((~0)>>1) && is_prob_prime($p*2+1);
  my $mp = BONE->copy->blsft($p)->bdec;

  # Definitely faster than using Math::BigInt
  return (0 == (Math::Prime::Util::GMP::lucas_sequence($mp, 4, 1, $mp+1))[0])
    if $Math::Prime::Util::_GMPfunc{"lucas_sequence"};

  my $V = Math::BigInt->new(4);
  for my $k (3 .. $p) {
    $V->bmul($V)->bsub(BTWO)->bmod($mp);
  }
  return $V->is_zero;
}


my $_poly_bignum;
sub _poly_new {
  my @poly = @_;
  push @poly, 0 unless scalar @poly;
  if ($_poly_bignum) {
    @poly = map { (ref $_ eq 'Math::BigInt')
                  ?  $_->copy
                  :  Math::BigInt->new("$_"); } @poly;
  }
  return \@poly;
}

#sub _poly_print {
#  my($poly) = @_;
#  carp "poly has null top degree" if $#$poly > 0 && !$poly->[-1];
#  foreach my $d (reverse 1 .. $#$poly) {
#    my $coef = $poly->[$d];
#    print "", ($coef != 1) ? $coef : "", ($d > 1) ? "x^$d" : "x", " + "
#      if $coef;
#  }
#  my $p0 = $poly->[0] || 0;
#  print "$p0\n";
#}

sub _poly_mod_mul {
  my($px, $py, $r, $n) = @_;

  my $px_degree = $#$px;
  my $py_degree = $#$py;
  my @res = map { $_poly_bignum ? Math::BigInt->bzero : 0 } 0 .. $r-1;

  # convolve(px, py) mod (X^r-1,n)
  my @indices_y = grep { $py->[$_] } (0 .. $py_degree);
  foreach my $ix (0 .. $px_degree) {
    my $px_at_ix = $px->[$ix];
    next unless $px_at_ix;
    if ($_poly_bignum) {
      foreach my $iy (@indices_y) {
        my $rindex = ($ix + $iy) % $r;  # reduce mod X^r-1
        $res[$rindex]->badd($px_at_ix->copy->bmul($py->[$iy]))->bmod($n);
      }
    } else {
      foreach my $iy (@indices_y) {
        my $rindex = ($ix + $iy) % $r;  # reduce mod X^r-1
        $res[$rindex] = ($res[$rindex] + $px_at_ix * $py->[$iy]) % $n;
      }
    }
  }
  # In case we had upper terms go to zero after modulo, reduce the degree.
  pop @res while !$res[-1];
  return \@res;
}

sub _poly_mod_pow {
  my($pn, $power, $r, $mod) = @_;
  my $res = _poly_new(1);
  my $p = $power;

  while ($p) {
    $res = _poly_mod_mul($res, $pn, $r, $mod) if ($p & 1);
    $p >>= 1;
    $pn  = _poly_mod_mul($pn,  $pn, $r, $mod) if $p;
  }
  return $res;
}

sub _test_anr {
  my($a, $n, $r) = @_;
  my $pp = _poly_mod_pow(_poly_new($a, 1), $n, $r, $n);
  $pp->[$n % $r] = (($pp->[$n % $r] || 0) -  1) % $n;  # subtract X^(n%r)
  $pp->[      0] = (($pp->[      0] || 0) - $a) % $n;  # subtract a
  return 0 if scalar grep { $_ } @$pp;
  1;
}

sub is_aks_prime {
  my $n = shift;
  return 0 if $n < 2 || is_power($n);

  my($log2n, $limit);
  if ($n > 2**48) {
    do { require Math::BigFloat; Math::BigFloat->import(); }
      if !defined $Math::BigFloat::VERSION;
    # limit = floor( log2(n) * log2(n) ).  o_r(n) must be larger than this
    my $floatn = Math::BigFloat->new("$n");
    #my $sqrtn = _bigint_to_int($floatn->copy->bsqrt->bfloor);
    # The following line seems to trigger a memory leak in Math::BigFloat::blog
    # (the part where $MBI is copied to $int) if $n is a Math::BigInt::GMP.
    $log2n = $floatn->copy->blog(2);
    $limit = _bigint_to_int( ($log2n * $log2n)->bfloor );
  } else {
    $log2n = log($n)/log(2) + 0.0001;      # Error on large side.
    $limit = int( $log2n*$log2n + 0.0001 );
  }

  my $r = next_prime($limit);
  foreach my $f (@{primes(0,$r-1)}) {
    return 1 if $f == $n;
    return 0 if !($n % $f);
  }

  while ($r < $n) {
    return 0 if !($n % $r);
    #return 1 if $r >= $sqrtn;
    last if znorder($n, $r) > $limit;  # Note the arguments!
    $r = next_prime($r);
  }

  return 1 if $r >= $n;

  # Since r is a prime, phi(r) = r-1
  my $rlimit = (ref($log2n) eq 'Math::BigFloat')
             ? _bigint_to_int( Math::BigFloat->new("$r")->bdec()
                                           ->bsqrt->bmul($log2n)->bfloor)
             : int( (sqrt(($r-1)) * $log2n) + 0.001 );

  $_poly_bignum = 1;
  if ( $n < (MPU_HALFWORD-1) ) {
    $_poly_bignum = 0;
    #$n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt';
  } else {
    $n = Math::BigInt->new("$n") unless ref($n) eq 'Math::BigInt';
  }

  my $_verbose = Math::Prime::Util::prime_get_config()->{'verbose'};
  print "# aks r = $r  s = $rlimit\n" if $_verbose;
  local $| = 1 if $_verbose > 1;
  for (my $a = 1; $a <= $rlimit; $a++) {
    return 0 unless _test_anr($a, $n, $r);
    print "." if $_verbose > 1;
  }
  print "\n" if $_verbose > 1;

  return 1;
}


sub _basic_factor {
  # MODIFIES INPUT SCALAR
  return ($_[0] == 1) ? () : ($_[0])   if $_[0] < 4;

  my @factors;
  if (ref($_[0]) ne 'Math::BigInt') {
    while ( !($_[0] % 2) ) { push @factors, 2;  $_[0] = int($_[0] / 2); }
    while ( !($_[0] % 3) ) { push @factors, 3;  $_[0] = int($_[0] / 3); }
    while ( !($_[0] % 5) ) { push @factors, 5;  $_[0] = int($_[0] / 5); }
  } else {
    # Without this, the bdivs will try to convert the results to BigFloat
    # and lose precision.
    $_[0]->upgrade(undef) if ref($_[0]) && $_[0]->upgrade();
    if (!Math::BigInt::bgcd($_[0], B_PRIM235)->is_one) {
      while ( $_[0]->is_even)   { push @factors, 2;  $_[0]->brsft(BONE); }
      foreach my $div (3, 5) {
        my ($q, $r) = $_[0]->copy->bdiv($div);
        while ($r->is_zero) {
          push @factors, $div;
          $_[0] = $q;
          ($q, $r) = $_[0]->copy->bdiv($div);
        }
      }
    }
    $_[0] = _bigint_to_int($_[0]) if $] >= 5.008 && $_[0] <= BMAX;
  }

  if ( ($_[0] > 1) && _is_prime7($_[0]) ) {
    push @factors, $_[0];
    $_[0] = 1;
  }
  @factors;
}

sub trial_factor {
  my($n, $limit) = @_;

  # Don't use _basic_factor here -- they want a trial forced.
  my @factors;
  if ($n < 4) {
    @factors = ($n == 1) ? () : ($n);
    return @factors;
  }

  my $start_idx = 1;
  # Expand small primes if it would help.
  push @_primes_small, @{primes($_primes_small[-1]+1, 100_003)}
    if $n > 400_000_000
    && $_primes_small[-1] < 99_000
    && (!defined $limit || $limit > $_primes_small[-1]);

  # Do initial bigint reduction.  Hopefully reducing it to native int.
  if (ref($n) eq 'Math::BigInt') {
    $n = $n->copy;  # Don't modify their original input!
    my $newlim = $n->copy->bsqrt;
    $limit = $newlim if !defined $limit || $limit > $newlim;
    while ($start_idx <= $#_primes_small) {
      my $f = $_primes_small[$start_idx++];
      last if $f > $limit;
      if ($n->copy->bmod($f)->is_zero) {
        do {
          push @factors, $f;
          $n->bdiv($f)->bfloor();
        } while $n->copy->bmod($f)->is_zero;
        last if $n < BMAX;
        my $newlim = $n->copy->bsqrt;
        $limit = $newlim if $limit > $newlim;
      }
    }
    return @factors if $n->is_one;
    $n = _bigint_to_int($n) if $n <= BMAX;
    return (@factors,$n) if $start_idx <= $#_primes_small && $_primes_small[$start_idx] > $limit;
  }

  {
    my $newlim = (ref($n) eq 'Math::BigInt') ? $n->copy->bsqrt : int(sqrt($n) + 0.001);
    $limit = $newlim if !defined $limit || $limit > $newlim;
  }

  if (ref($n) ne 'Math::BigInt') {
    for my $i ($start_idx .. $#_primes_small) {
      my $p = $_primes_small[$i];
      last if $p > $limit;
      if (($n % $p) == 0) {
        do { push @factors, $p;  $n = int($n/$p); } while ($n % $p) == 0;
        last if $n == 1;
        my $newlim = int( sqrt($n) + 0.001);
        $limit = $newlim if $newlim < $limit;
      }
    }
    if ($_primes_small[-1] < $limit) {
      my $inc = (($_primes_small[-1] % 6) == 1) ? 4 : 2;
      my $p = $_primes_small[-1] + $inc;
      while ($p <= $limit) {
        if (($n % $p) == 0) {
          do { push @factors, $p;  $n = int($n/$p); } while ($n % $p) == 0;
          last if $n == 1;
          my $newlim = int( sqrt($n) + 0.001);
          $limit = $newlim if $newlim < $limit;
        }
        $p += ($inc ^= 6);
      }
    }
  } else {   # n is a bigint.  Use mod-210 wheel trial division.
    # Generating a wheel mod $w starting at $s:
    # mpu 'my($s,$w,$t)=(11,2*3*5); say join ",",map { ($t,$s)=($_-$s,$_); $t; } grep { gcd($_,$w)==1 } $s+1..$s+$w;'
    # Should start at $_primes_small[$start_idx], do 11 + next multiple of 210.
    my @incs = map { Math::BigInt->new($_) } (2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10);
    my $f = 11; while ($f <= $_primes_small[$start_idx-1]-210) { $f += 210; }
    ($f, $limit) = map { Math::BigInt->new("$_") } ($f, $limit);
    SEARCH: while ($f <= $limit) {
      foreach my $finc (@incs) {
        if ($n->copy->bmod($f)->is_zero && $f->bacmp($limit) <= 0) {
          my $sf = ($f <= BMAX) ? _bigint_to_int($f) : $f->copy;
          do {
            push @factors, $sf;
            $n->bdiv($f)->bfloor();
          } while $n->copy->bmod($f)->is_zero;
          last SEARCH if $n->is_one;
          my $newlim = $n->copy->bsqrt;
          $limit = $newlim if $limit > $newlim;
        }
        $f->badd($finc);
      }
    }
  }
  push @factors, $n  if $n > 1;
  @factors;
}

my $_holf_r;
my @_fsublist = (
  [ "pbrent 32k", sub { pbrent_factor (shift,   32*1024, 1, 1) } ],
  [ "p-1 1M",     sub { pminus1_factor(shift, 1_000_000, undef, 1); } ],
  [ "ECM 1k",     sub { ecm_factor    (shift,     1_000,   5_000, 15) } ],
  [ "pbrent 512k",sub { pbrent_factor (shift,  512*1024, 7, 1) } ],
  [ "p-1 4M",     sub { pminus1_factor(shift, 4_000_000, undef, 1); } ],
  [ "ECM 10k",    sub { ecm_factor    (shift,    10_000,  50_000, 10) } ],
  [ "pbrent 512k",sub { pbrent_factor (shift,  512*1024, 11, 1) } ],
  [ "HOLF 256k",  sub { holf_factor   (shift, 256*1024, $_holf_r); $_holf_r += 256*1024; } ],
  [ "p-1 20M",    sub { pminus1_factor(shift,20_000_000); } ],
  [ "ECM 100k",   sub { ecm_factor    (shift,   100_000, 800_000, 10) } ],
  [ "HOLF 512k",  sub { holf_factor   (shift, 512*1024, $_holf_r); $_holf_r += 512*1024; } ],
  [ "pbrent 2M",  sub { pbrent_factor (shift, 2048*1024, 13, 1) } ],
  [ "HOLF 2M",    sub { holf_factor   (shift, 2048*1024, $_holf_r); $_holf_r += 2048*1024; } ],
  [ "ECM 1M",     sub { ecm_factor    (shift, 1_000_000, 1_000_000, 10) } ],
  [ "p-1 100M",   sub { pminus1_factor(shift, 100_000_000, 500_000_000); } ],
);

sub factor {
  my($n) = @_;
  _validate_positive_integer($n);
  my @factors;

  if ($n < 4) {
    @factors = ($n == 1) ? () : ($n);
    return @factors;
  }
  $n = $n->copy if ref($n) eq 'Math::BigInt';
  my $lim = 4999;  # How much trial factoring to do

  # For native integers, we could save a little time by doing hardcoded trials
  # by 2-29 here.  Skipping it.

  push @factors, trial_factor($n, $lim);
  return @factors if $factors[-1] < $lim*$lim;
  $n = pop(@factors);

  my @nstack = ($n);
  while (@nstack) {
    $n = pop @nstack;
    # Don't use bignum on $n if it has gotten small enough.
    $n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt' && $n <= BMAX;
    #print "Looking at $n with stack ", join(",",@nstack), "\n";
    while ( ($n >= ($lim*$lim)) && !_is_prime7($n) ) {
      my @ftry;
      $_holf_r = 1;
      foreach my $sub (@_fsublist) {
        last if scalar @ftry >= 2;
        print "  starting $sub->[0]\n" if Math::Prime::Util::prime_get_config()->{'verbose'} > 1;
        @ftry = $sub->[1]->($n);
      }
      if (scalar @ftry > 1) {
        #print "  split into ", join(",",@ftry), "\n";
        $n = shift @ftry;
        $n = _bigint_to_int($n) if ref($n) eq 'Math::BigInt' && $n <= BMAX;
        push @nstack, @ftry;
      } else {
        #warn "trial factor $n\n";
        push @factors, trial_factor($n);
        #print "  trial into ", join(",",@factors), "\n";
        $n = 1;
        last;
      }
    }
    push @factors, $n  if $n != 1;
  }
  @factors = sort {$a<=>$b} @factors;
  return @factors;
}

sub _found_factor {
  my($f, $n, $what, @factors) = @_;
  if ($f == 1 || $f == $n) {
    push @factors, $n;
  } else {
    # Perl 5.6.2 needs things spelled out for it.
    my $f2 = (ref($n) eq 'Math::BigInt') ? $n->copy->bdiv($f)->as_int
                                         : int($n/$f);
    push @factors, $f;
    push @factors, $f2;
    croak "internal error in $what" unless $f * $f2 == $n;
    # MPU::GMP prints this type of message if verbose, so do the same.
    print "$what found factor $f\n" if Math::Prime::Util::prime_get_config()->{'verbose'} > 0;
  }
  @factors;
}

# TODO:
sub squfof_factor { trial_factor(@_) }

sub prho_factor {
  my($n, $rounds, $pa, $skipbasic) = @_;
  $rounds = 4*1024*1024 unless defined $rounds;
  $pa = 3 unless defined $pa;

  my @factors;
  if (!$skipbasic) {
    @factors = _basic_factor($n);
    return @factors if $n < 4;
  }

  my $inloop = 0;
  my $U = 7;
  my $V = 7;

  if ( ref($n) eq 'Math::BigInt' ) {

    my $zero = $n->copy->bzero;
    $pa = $zero->badd("$pa");
    $U = $zero->copy->badd($U);
    $V = $zero->copy->badd($V);
    for my $i (1 .. $rounds) {
      # Would use bmuladd here, but old Math::BigInt's barf with scalar $pa.
      $U->bmul($U)->badd($pa)->bmod($n);
      $V->bmul($V)->badd($pa);
      $V->bmul($V)->badd($pa)->bmod($n);
      my $f = Math::BigInt::bgcd($U-$V, $n);
      if ($f->bacmp($n) == 0) {
        last if $inloop++;  # We've been here before
      } elsif (!$f->is_one) {
        return _found_factor($f, $n, "prho", @factors);
      }
    }

  } elsif ($n < MPU_HALFWORD) {

    my $inner = 32;
    $rounds = int( ($rounds + $inner-1) / $inner );
    while ($rounds-- > 0) {
      my($m, $oldU, $oldV, $f) = (1, $U, $V);
      for my $i (1 .. $inner) {
        $U = ($U * $U + $pa) % $n;
        $V = ($V * $V + $pa) % $n;
        $V = ($V * $V + $pa) % $n;
        $f = ($U > $V) ? $U-$V : $V-$U;
        $m = ($m * $f) % $n;
      }
      $f = _gcd_ui( $m, $n );
      next if $f == 1;
      if ($f == $n) {
        ($U, $V) = ($oldU, $oldV);
        for my $i (1 .. $inner) {
          $U = ($U * $U + $pa) % $n;
          $V = ($V * $V + $pa) % $n;
          $V = ($V * $V + $pa) % $n;
          $f = ($U > $V) ? $U-$V : $V-$U;
          $f = _gcd_ui( $f, $n);
          last if $f != 1;
        }
        last if $f == 1 || $f == $n;
      }
      return _found_factor($f, $n, "prho", @factors);
    }

  } else {

    for my $i (1 .. $rounds) {
      if ($n <= (~0 >> 1)) {
       $U = _mulmod($U, $U, $n);  $U += $pa;  $U -= $n if $U >= $n;
       $V = _mulmod($V, $V, $n);  $V += $pa;  # Let the mulmod handle it
       $V = _mulmod($V, $V, $n);  $V += $pa;  $V -= $n if $V >= $n;
      } else {
       #$U = _mulmod($U, $U, $n); $U=$n-$U; $U = ($pa>=$U) ? $pa-$U : $n-$U+$pa;
       #$V = _mulmod($V, $V, $n); $V=$n-$V; $V = ($pa>=$V) ? $pa-$V : $n-$V+$pa;
       #$V = _mulmod($V, $V, $n); $V=$n-$V; $V = ($pa>=$V) ? $pa-$V : $n-$V+$pa;
       $U = _mulmod($U, $U, $n);  $U = _addmod($U, $pa, $n);
       $V = _mulmod($V, $V, $n);  $V = _addmod($V, $pa, $n);
       $V = _mulmod($V, $V, $n);  $V = _addmod($V, $pa, $n);
      }
      my $f = _gcd_ui( $U-$V,  $n );
      if ($f == $n) {
        last if $inloop++;  # We've been here before
      } elsif ($f != 1) {
        return _found_factor($f, $n, "prho", @factors);
      }
    }

  }
  push @factors, $n;
  @factors;
}

sub pbrent_factor {
  my($n, $rounds, $pa, $skipbasic) = @_;
  $rounds = 4*1024*1024 unless defined $rounds;
  $pa = 3 unless defined $pa;

  my @factors;
  if (!$skipbasic) {
    @factors = _basic_factor($n);
    return @factors if $n < 4;
  }

  my $Xi = 2;
  my $Xm = 2;

  if ( ref($n) eq 'Math::BigInt' ) {

    # Same code as the GMP version, but runs *much* slower.  Even with
    # Math::BigInt::GMP it's >200x slower.  With the default Calc backend
    # it's thousands of times slower.
    my $inner = 32;
    my $zero = $n->copy->bzero;
    my $saveXi;
    my $f;
    $Xi = $zero->copy->badd($Xi);
    $Xm = $zero->copy->badd($Xm);
    $pa = $zero->copy->badd($pa);
    my $r = 1;
    while ($rounds > 0) {
      my $rleft = ($r > $rounds) ? $rounds : $r;
      while ($rleft > 0) {
        my $dorounds = ($rleft > $inner) ? $inner : $rleft;
        my $m = $zero->copy->bone;
        $saveXi = $Xi->copy;
        foreach my $i (1 .. $dorounds) {
          $Xi->bmul($Xi)->badd($pa)->bmod($n);
          $m->bmul($Xi->copy->bsub($Xm));
        }
        $rleft -= $dorounds;
        $rounds -= $dorounds;
        $m->bmod($n);
        $f = Math::BigInt::bgcd($m,  $n);
        last unless $f->is_one;
      }
      if ($f->is_one) {
        $r *= 2;
        $Xm = $Xi->copy;
        next;
      }
      if ($f == $n) {  # back up to determine the factor
        $Xi = $saveXi->copy;
        do {
          $Xi->bmul($Xi)->badd($pa)->bmod($n);
          $f = Math::BigInt::bgcd($Xm-$Xi, $n);
        } while ($f != 1 && $r-- != 0);
        last if $f == 1 || $f == $n;
      }
      return _found_factor($f, $n, "pbrent", @factors);
    }

  } elsif ($n < MPU_HALFWORD) {

    # Doing the gcd batching as above works pretty well here, but it's a lot
    # of code for not much gain for general users.
    for my $i (1 .. $rounds) {
      $Xi = ($Xi * $Xi + $pa) % $n;
      my $f = _gcd_ui( ($Xi>$Xm) ? $Xi-$Xm : $Xm-$Xi, $n);
      return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
      $Xm = $Xi if ($i & ($i-1)) == 0;  # i is a power of 2
    }

  } else {

    for my $i (1 .. $rounds) {
      $Xi = _addmod( _mulmod($Xi, $Xi, $n), $pa, $n);
      my $f = _gcd_ui( ($Xi>$Xm) ? $Xi-$Xm : $Xm-$Xi, $n);
      return _found_factor($f, $n, "pbrent", @factors) if $f != 1 && $f != $n;
      $Xm = $Xi if ($i & ($i-1)) == 0;  # i is a power of 2
    }

  }
  push @factors, $n;
  @factors;
}

sub pminus1_factor {
  my($n, $B1, $B2, $skipbasic) = @_;

  my @factors;
  if (!$skipbasic) {
    @factors = _basic_factor($n);
    return @factors if $n < 4;
  }

  if ( ref($n) ne 'Math::BigInt' ) {
    # Stage 1 only
    $B1 = 10_000_000 unless defined $B1;
    my $pa = 2;
    my $f = 1;
    my($pc_beg, $pc_end, @bprimes);
    $pc_beg = 2;
    $pc_end = $pc_beg + 100_000;
    my $sqrtb1 = int(sqrt($B1));
    while (1) {
      $pc_end = $B1 if $pc_end > $B1;
      @bprimes = @{ primes($pc_beg, $pc_end) };
      foreach my $q (@bprimes) {
        my $k = $q;
        if ($q <= $sqrtb1) {
          my $kmin = int($B1 / $q);
          while ($k <= $kmin) { $k *= $q; }
        }
        $pa = _powmod($pa, $k, $n);
        if ($pa == 0) { push @factors, $n; return @factors; }
        my $f = _gcd_ui( $pa-1, $n );
        return _found_factor($f, $n, "pminus1", @factors) if $f != 1;
      }
      last if $pc_end >= $B1;
      $pc_beg = $pc_end+1;
      $pc_end += 500_000;
    }
    push @factors, $n;
    return @factors;
  }

  # Stage 2 isn't really any faster than stage 1 for the examples I've tried.
  # Perl's overhead is greater than the savings of multiply vs. powmod

  if (!defined $B1) {
    for my $mul (1, 100, 1000, 10_000, 100_000, 1_000_000) {
      $B1 = 1000 * $mul;
      $B2 = 1*$B1;
      #warn "Trying p-1 with $B1 / $B2\n";
      my @nf = pminus1_factor($n, $B1, $B2);
      if (scalar @nf > 1) {
        push @factors, @nf;
        return @factors;
      }
    }
    push @factors, $n;
    return @factors;
  }
  $B2 = 1*$B1 unless defined $B2;

  my $one = $n->copy->bone;
  my ($j, $q, $saveq) = (32, 2, 2);
  my $t = $one->copy;
  my $pa = $one->copy->binc();
  my $savea = $pa->copy;
  my $f = $one->copy;
  my($pc_beg, $pc_end, @bprimes);

  $pc_beg = 2;
  $pc_end = $pc_beg + 100_000;
  while (1) {
    $pc_end = $B1 if $pc_end > $B1;
    @bprimes = @{ primes($pc_beg, $pc_end) };
    foreach my $q (@bprimes) {
      my($k, $kmin) = ($q, int($B1 / $q));
      while ($k <= $kmin) { $k *= $q; }
      $t *= $k;                         # accumulate powers for a
      if ( ($j++ % 64) == 0) {
        next if $pc_beg > 2 && ($j-1) % 256;
        $pa->bmodpow($t, $n);
        $t = $one->copy;
        if ($pa == 0) { push @factors, $n; return @factors; }
        $f = Math::BigInt::bgcd( $pa->copy->bdec, $n );
        last if $f == $n;
        return _found_factor($f, $n, "pminus1", @factors) unless $f->is_one;
        $saveq = $q;
        $savea = $pa->copy;
      }
    }
    $q = $bprimes[-1];
    last if !$f->is_one || $pc_end >= $B1;
    $pc_beg = $pc_end+1;
    $pc_end += 500_000;
  }
  undef @bprimes;
  $pa->bmodpow($t, $n);
  if ($pa == 0) { push @factors, $n; return @factors; }
  $f = Math::BigInt::bgcd( $pa-1, $n );
  if ($f == $n) {
    $q = $saveq;
    $pa = $savea->copy;
    while ($q <= $B1) {
      my ($k, $kmin) = ($q, int($B1 / $q));
      while ($k <= $kmin) { $k *= $q; }
      $pa->bmodpow($k, $n);
      my $f = Math::BigInt::bgcd( $pa-1, $n );
      if ($f == $n) { push @factors, $n; return @factors; }
      last if !$f->is_one;
      $q = next_prime($q);
    }
  }
  # STAGE 2
  if ($f->is_one && $B2 > $B1) {
    my $bm = $pa->copy;
    my $b = $one->copy;
    my @precomp_bm;
    $precomp_bm[0] = ($bm * $bm) % $n;
    foreach my $j (1..19) {
      $precomp_bm[$j] = ($precomp_bm[$j-1] * $bm * $bm) % $n;
    }
    $pa->bmodpow($q, $n);
    my $j = 1;
    $pc_beg = $q+1;
    $pc_end = $pc_beg + 100_000;
    while (1) {
      $pc_end = $B2 if $pc_end > $B2;
      @bprimes = @{ primes($pc_beg, $pc_end) };
      foreach my $i (0 .. $#bprimes) {
        my $diff = $bprimes[$i] - $q;
        $q = $bprimes[$i];
        my $qdiff = ($diff >> 1) - 1;
        if (!defined $precomp_bm[$qdiff]) {
          $precomp_bm[$qdiff] = $bm->copy->bmodpow($diff, $n);
        }
        $pa->bmul($precomp_bm[$qdiff])->bmod($n);
        if ($pa == 0) { push @factors, $n; return @factors; }
        $b->bmul($pa-1);
        if (($j++ % 128) == 0) {
          $b->bmod($n);
          $f = Math::BigInt::bgcd( $b, $n );
          last if !$f->is_one;
        }
      }
      last if !$f->is_one || $pc_end >= $B2;
      $pc_beg = $pc_end+1;
      $pc_end += 500_000;
    }
    $f = Math::BigInt::bgcd( $b, $n );
  }
  return _found_factor($f, $n, "pminus1", @factors);
}

sub holf_factor {
  my($n, $rounds, $startrounds) = @_;
  $rounds = 64*1024*1024 unless defined $rounds;
  $startrounds = 1 unless defined $startrounds;
  $startrounds = 1 if $startrounds < 1;

  my @factors = _basic_factor($n);
  return @factors if $n < 4;

  if ( ref($n) eq 'Math::BigInt' ) {
    for my $i ($startrounds .. $rounds) {
      my $ni = $n->copy->bmul($i);
      my $s = $ni->copy->bsqrt->bfloor->as_int;
      if ($s * $s == $ni) {
        # s^2 = n*i, so m = s^2 mod n = 0.  Hence f = GCD(n, s) = GCD(n, n*i)
        my $f = Math::BigInt::bgcd($ni, $n);
        return _found_factor($f, $n, "HOLF", @factors);
      }
      $s->binc;
      my $m = ($s * $s) - $ni;
      # Check for perfect square
      my $mc = _bigint_to_int($m & 31);
      next unless $mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25;
      my $f = $m->copy->bsqrt->bfloor->as_int;
      next unless ($f*$f) == $m;
      $f = Math::BigInt::bgcd( ($s > $f) ? $s-$f : $f-$s,  $n);
      return _found_factor($f, $n, "HOLF ($i rounds)", @factors);
    }
  } else {
    for my $i ($startrounds .. $rounds) {
      my $s = int(sqrt($n * $i));
      $s++ if ($s * $s) != ($n * $i);
      my $m = ($s < MPU_HALFWORD) ? ($s*$s) % $n : _mulmod($s, $s, $n);
      # Check for perfect square
      my $mc = $m & 31;
      next unless $mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25;
      my $f = int(sqrt($m));
      next unless $f*$f == $m;
      $f = _gcd_ui($s - $f,  $n);
      return _found_factor($f, $n, "HOLF ($i rounds)", @factors);
    }
  }
  push @factors, $n;
  @factors;
}

sub fermat_factor {
  my($n, $rounds) = @_;
  $rounds = 64*1024*1024 unless defined $rounds;

  my @factors = _basic_factor($n);
  return @factors if $n < 4;

  if ( ref($n) eq 'Math::BigInt' ) {
    my $pa = $n->copy->bsqrt->bfloor->as_int;
    return _found_factor($pa, $n, "Fermat", @factors) if $pa*$pa == $n;
    $pa++;
    my $b2 = $pa*$pa - $n;
    my $lasta = $pa + $rounds;
    while ($pa <= $lasta) {
      my $mc = _bigint_to_int($b2 & 31);
      if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
        my $s = $b2->copy->bsqrt->bfloor->as_int;
        if ($s*$s == $b2) {
          my $i = $pa-($lasta-$rounds)+1;
          return _found_factor($pa - $s, $n, "Fermat ($i rounds)", @factors);
        }
      }
      $pa++;
      $b2 = $pa*$pa-$n;
    }
  } else {
    my $pa = int(sqrt($n));
    return _found_factor($pa, $n, "Fermat", @factors) if $pa*$pa == $n;
    $pa++;
    my $b2 = $pa*$pa - $n;
    my $lasta = $pa + $rounds;
    while ($pa <= $lasta) {
      my $mc = $b2 & 31;
      if ($mc==0||$mc==1||$mc==4||$mc==9||$mc==16||$mc==17||$mc==25) {
        my $s = int(sqrt($b2));
        if ($s*$s == $b2) {
          my $i = $pa-($lasta-$rounds)+1;
          return _found_factor($pa - $s, $n, "Fermat ($i rounds)", @factors);
        }
      }
      $pa++;
      $b2 = $pa*$pa-$n;
    }
  }
  push @factors, $n;
  @factors;
}


sub ecm_factor {
  my($n, $B1, $B2, $ncurves) = @_;
  _validate_positive_integer($n);

  my @factors = _basic_factor($n);
  return @factors if $n < 4;

  if ($Math::Prime::Util::_GMPfunc{"ecm_factor"}) {
    $B1 = 0 if !defined $B1;
    $ncurves = 0 if !defined $ncurves;
    my @ef = Math::Prime::Util::GMP::ecm_factor($n, $B1, $ncurves);
    if (@ef > 1) {
      my $ecmfac = Math::Prime::Util::_reftyped($n, $ef[-1]);
      return _found_factor($ecmfac, $n, "ECM (GMP) B1=$B1 curves $ncurves", @factors);
    }
    push @factors, $n;
    return @factors;
  }

  $ncurves = 10 unless defined $ncurves;

  if (!defined $B1) {
    for my $mul (1, 10, 100, 1000, 10_000, 100_000, 1_000_000) {
      $B1 = 100 * $mul;
      $B2 = 10*$B1;
      #warn "Trying ecm with $B1 / $B2\n";
      my @nf = ecm_factor($n, $B1, $B2, $ncurves);
      if (scalar @nf > 1) {
        push @factors, @nf;
        return @factors;
      }
    }
    push @factors, $n;
    return @factors;
  }

  $B2 = 10*$B1 unless defined $B2;
  my $sqrt_b1 = int(sqrt($B1)+1);

  # Affine code.  About 3x slower than the projective, and no stage 2.
  #
  #if (!defined $Math::Prime::Util::ECAffinePoint::VERSION) {
  #  eval { require Math::Prime::Util::ECAffinePoint; 1; }
  #  or do { croak "Cannot load Math::Prime::Util::ECAffinePoint"; };
  #}
  #my @bprimes = @{ primes(2, $B1) };
  #my $irandf = Math::Prime::Util::_get_rand_func();
  #foreach my $curve (1 .. $ncurves) {
  #  my $a = $irandf->($n-1);
  #  my $b = 1;
  #  my $ECP = Math::Prime::Util::ECAffinePoint->new($a, $b, $n, 0, 1);
  #  foreach my $q (@bprimes) {
  #    my $k = $q;
  #    if ($k < $sqrt_b1) {
  #      my $kmin = int($B1 / $q);
  #      while ($k <= $kmin) { $k *= $q; }
  #    }
  #    $ECP->mul($k);
  #    my $f = $ECP->f;
  #    if ($f != 1) {
  #      last if $f == $n;
  #      warn "ECM found factors with B1 = $B1 in curve $curve\n";
  #      return _found_factor($f, $n, "ECM B1=$B1 curve $curve", @factors);
  #    }
  #    last if $ECP->is_infinity;
  #  }
  #}

  require Math::Prime::Util::ECProjectivePoint;
  require Math::Prime::Util::RandomPrimes;

  # With multiple curves, it's better to get all the primes at once.
  # The downside is this can kill memory with a very large B1.
  my @bprimes = @{ primes(3, $B1) };
  foreach my $q (@bprimes) {
    last if $q > $sqrt_b1;
    my($k,$kmin) = ($q, int($B1/$q));
    while ($k <= $kmin) { $k *= $q; }
    $q = $k;
  }
  my @b2primes = ($B2 > $B1) ? @{primes($B1+1, $B2)} : ();

  foreach my $curve (1 .. $ncurves) {
    my $sigma = Math::Prime::Util::urandomm($n-6) + 6;
    my ($u, $v) = ( ($sigma*$sigma - 5) % $n, (4 * $sigma) % $n );
    my ($x, $z) = ( ($u*$u*$u) % $n,  ($v*$v*$v) % $n );
    my $cb = (4 * $x * $v) % $n;
    my $ca = ( (($v-$u)**3) * (3*$u + $v) ) % $n;
    my $f = Math::BigInt::bgcd( $cb, $n );
    $f = Math::BigInt::bgcd( $z, $n ) if $f == 1;
    next if $f == $n;
    return _found_factor($f,$n, "ECM B1=$B1 curve $curve", @factors) if $f != 1;
    $cb = Math::BigInt->new("$cb") unless ref($cb) eq 'Math::BigInt';
    $u = $cb->copy->bmodinv($n);
    $ca = (($ca*$u) - 2) % $n;

    my $ECP = Math::Prime::Util::ECProjectivePoint->new($ca, $n, $x, $z);
    my $fm = $n-$n+1;
    my $i = 15;

    for (my $q = 2; $q < $B1; $q *= 2) { $ECP->double(); }
    foreach my $k (@bprimes) {
      $ECP->mul($k);
      $fm = ($fm * $ECP->x() ) % $n;
      if ($i++ % 32 == 0) {
        $f = Math::BigInt::bgcd($fm, $n);
        last if $f != 1;
      }
    }
    $f = Math::BigInt::bgcd($fm, $n);
    next if $f == $n;

    if ($f == 1 && $B2 > $B1) { # BEGIN STAGE 2
      my $D = int(sqrt($B2/2));  $D++ if $D % 2;
      my $one = $n - $n + 1;
      my $g = $one;

      my $S2P = $ECP->copy->normalize;
      $f = $S2P->f;
      if ($f != 1) {
        next if $f == $n;
        #warn "ECM S2 normalize f=$f\n" if $f != 1;
        return _found_factor($f, $n, "ECM S2 B1=$B1 curve $curve");
      }
      my $S2x = $S2P->x;
      my $S2d = $S2P->d;
      my @nqx = ($n-$n, $S2x);

      foreach my $i (2 .. 2*$D) {
        my($x2, $z2);
        if ($i % 2) {
          ($x2, $z2) = Math::Prime::Util::ECProjectivePoint::_addx($nqx[($i-1)/2], $nqx[($i+1)/2], $S2x, $n);
        } else {
          ($x2, $z2) = Math::Prime::Util::ECProjectivePoint::_double($nqx[$i/2], $one, $n, $S2d);
        }
        $nqx[$i] = $x2;
        #($f, $u, undef) = _extended_gcd($z2, $n);
        $f = Math::BigInt::bgcd( $z2, $n );
        last if $f != 1;
        $u = $z2->copy->bmodinv($n);
        $nqx[$i] = ($x2 * $u) % $n;
      }
      if ($f != 1) {
        next if $f == $n;
        #warn "ECM S2 1: B1 $B1 B2 $B2 curve $curve f=$f\n";
        return _found_factor($f, $n, "ECM S2 B1=$B1 curve $curve", @factors);
      }

      $x = $nqx[2*$D-1];
      my $m = 1;
      while ($m < ($B2+$D)) {
        if ($m != 1) {
          my $oldx = $S2x;
          my ($x1, $z1) = Math::Prime::Util::ECProjectivePoint::_addx($nqx[2*$D], $S2x, $x, $n);
          $f = Math::BigInt::bgcd( $z1, $n );
          last if $f != 1;
          $u = $z1->copy->bmodinv($n);
          $S2x = ($x1 * $u) % $n;
          $x = $oldx;
          last if $f != 1;
        }
        if ($m+$D > $B1) {
          my @p = grep { $_ >= $m-$D && $_ <= $m+$D } @b2primes;
          foreach my $i (@p) {
            last if $i >= $m;
            $g = ($g * ($S2x - $nqx[$m+$D-$i])) % $n;
          }
          foreach my $i (@p) {
            next unless $i > $m;
            next if $i > ($m+$m) || is_prime($m+$m-$i);
            $g = ($g * ($S2x - $nqx[$i-$m])) % $n;
          }
          $f = Math::BigInt::bgcd($g, $n);
          #warn "ECM S2 3: found $f in stage 2\n" if $f != 1;
          last if $f != 1;
        }
        $m += 2*$D;
      }
    } # END STAGE 2

    next if $f == $n;
    if ($f != 1) {
      #warn "ECM found factors with B1 = $B1 in curve $curve\n";
      return _found_factor($f, $n, "ECM B1=$B1 curve $curve", @factors);
    }
    # end of curve loop
  }
  push @factors, $n;
  @factors;
}

sub divisors {
  my($n) = @_;
  _validate_positive_integer($n);
  my(@factors, @d, @t);

  # In scalar context, returns sigma_0(n).  Very fast.
  return Math::Prime::Util::divisor_sum($n,0) unless wantarray;
  return ($n == 0) ? (0,1) : (1)  if $n <= 1;

  if ($Math::Prime::Util::_GMPfunc{"divisors"}) {
    # This trips an erroneous compile time error without the eval.
    eval ' @d = Math::Prime::Util::GMP::divisors($n); ';  ## no critic qw(ProhibitStringyEval)
    @d = map { $_ <= ~0 ? $_ : ref($n)->new($_) } @d   if ref($n);
    return @d;
  }

  @factors = Math::Prime::Util::factor($n);
  return (1,$n) if scalar @factors == 1;

  my $bigint = ref($n);
  @factors = map { $bigint->new("$_") } @factors  if $bigint;
  @d = $bigint ? ($bigint->new(1)) : (1);

  while (my $p = shift @factors) {
    my $e = 1;
    while (@factors && $p == $factors[0]) { $e++; shift(@factors); }
    push @d,  @t = map { $_ * $p } @d;               # multiply through once
    push @d,  @t = map { $_ * $p } @t   for 2 .. $e; # repeat
  }

  @d = map { $_ <= INTMAX ? _bigint_to_int($_) : $_ } @d   if $bigint;
  @d = sort { $a <=> $b } @d;
  @d;
}


sub chebyshev_theta {
  my($n,$low) = @_;
  $low = 2 unless defined $low;
  my($sum,$high) = (0.0, 0);
  while ($low <= $n) {
    $high = $low + 1e6;
    $high = $n if $high > $n;
    $sum += log($_) for @{primes($low,$high)};
    $low = $high+1;
  }
  $sum;
}

sub chebyshev_psi {
  my($n) = @_;
  return 0 if $n <= 1;
  my ($sum, $logn, $sqrtn) = (0.0, log($n), int(sqrt($n)));

  # Sum the log of prime powers first
  for my $p (@{primes($sqrtn)}) {
    my $logp = log($p);
    $sum += $logp * int($logn/$logp+1e-15);
  }
  # The rest all have exponent 1: add them in using the segmenting theta code
  $sum += chebyshev_theta($n, $sqrtn+1);

  $sum;
}

sub hclassno {
  my $n = shift;

  return -1 if $n == 0;
  return 0 if $n < 0 || ($n % 4) == 1 || ($n % 4) == 2;
  return 2 * (2,3,6,6,6,8,12,9,6,12,18,12,8,12,18,18,12,15,24,12,6,24,30,20,12,12,24,24,18,24)[($n>>1)-1] if $n <= 60;

  my ($h, $square, $b, $b2) = (0, 0, $n & 1, ($n+1) >> 2);

  if ($b == 0) {
    my $lim = int(sqrt($b2));
    if (_is_perfect_square($b2)) {
      $square = 1;
      $lim--;
    }
    #$h += scalar(grep { $_ <= $lim } divisors($b2));
    for my $i (1 .. $lim) { $h++ unless $b2 % $i; }
    ($b,$b2) = (2, ($n+4) >> 2);
  }
  while ($b2 * 3 < $n) {
    $h++ unless $b2 % $b;
    my $lim = int(sqrt($b2));
    if (_is_perfect_square($b2)) {
      $h++;
      $lim--;
    }
    #$h += 2 * scalar(grep { $_ > $b && $_ <= $lim } divisors($b2));
    for my $i ($b+1 .. $lim) { $h += 2 unless $b2 % $i; }
    $b += 2;
    $b2 = ($n+$b*$b) >> 2;
  }
  return (($b2*3 == $n) ? 2*(3*$h+1) : $square ? 3*(2*$h+1) : 6*$h) << 1;
}

# Sigma method for prime powers
sub _taup {
  my($p, $e, $n) = @_;
  my($bp) = Math::BigInt->new("".$p);
  if ($e == 1) {
    return (0,1,-24,252,-1472,4830,-6048,-16744,84480)[$p] if $p <= 8;
    my $ds5  = $bp->copy->bpow( 5)->binc();  # divisor_sum(p,5)
    my $ds11 = $bp->copy->bpow(11)->binc();  # divisor_sum(p,11)
    my $s    = Math::BigInt->new("".vecsum(map { vecprod(BTWO,Math::Prime::Util::divisor_sum($_,5), Math::Prime::Util::divisor_sum($p-$_,5)) } 1..($p-1)>>1));
    $n = ( 65*$ds11 + 691*$ds5 - (691*252)*$s ) / 756;
  } else {
    my $t = Math::BigInt->new(""._taup($p,1));
    $n = $t->copy->bpow($e);
    if ($e == 2) {
      $n -= $bp->copy->bpow(11);
    } elsif ($e == 3) {
      $n -= BTWO * $t * $bp->copy->bpow(11);
    } else {
      $n += vecsum( map { vecprod( ($_&1) ? - BONE : BONE,
                                   $bp->copy->bpow(11*$_),
                                   binomial($e-$_, $e-2*$_),
                                   $t ** ($e-2*$_) ) } 1 .. ($e>>1) );
    }
  }
  $n = _bigint_to_int($n) if ref($n) && $n->bacmp(BMAX) <= 0;
  $n;
}

# Cohen's method using Hurwitz class numbers
# The two hclassno calls could be collapsed with some work
sub _tauprime {
  my $p = shift;
  return -24 if $p == 2;
  my $sum = Math::BigInt->new(0);
  if ($p < (MPU_32BIT ?  300  :  1600)) {
    my($p9,$pp7) = (9*$p, 7*$p*$p);
    for my $t (1 .. Math::Prime::Util::sqrtint($p)) {
      my $t2 = $t * $t;
      my $v = $p - $t2;
      $sum += $t2**3 * (4*$t2*$t2 - $p9*$t2 + $pp7) * (Math::Prime::Util::hclassno(4*$v) + 2 * Math::Prime::Util::hclassno($v));
    }
    $p = Math::BigInt->new("$p");
  } else {
    $p = Math::BigInt->new("$p");
    my($p9,$pp7) = (9*$p, 7*$p*$p);
    for my $t (1 .. Math::Prime::Util::sqrtint($p)) {
      my $t2 = Math::BigInt->new("$t") ** 2;
      my $v = $p - $t2;
      $sum += $t2**3 * (4*$t2*$t2 - $p9*$t2 + $pp7) * (Math::Prime::Util::hclassno(4*$v) + 2 * Math::Prime::Util::hclassno($v));
    }
  }
  28*$p**6 - 28*$p**5 - 90*$p**4 - 35*$p**3 - 1 - 32 * ($sum/3);
}

# Recursive method for handling prime powers
sub _taupower {
  my($p, $e) = @_;
  return 1 if $e <= 0;
  return _tauprime($p) if $e == 1;
  $p = Math::BigInt->new("$p");
  my($tp, $p11) = ( _tauprime($p), $p**11 );
  return $tp ** 2 - $p11 if $e == 2;
  return $tp ** 3 - 2 * $tp * $p11 if $e == 3;
  return $tp ** 4 - 3 * $tp**2 * $p11 + $p11**2 if $e == 4;
  # Recurse -3
  ($tp**3 - 2*$tp*$p11) * _taupower($p,$e-3) + ($p11*$p11 - $tp*$tp*$p11) * _taupower($p,$e-4);
}

sub ramanujan_tau {
  my $n = shift;
  return 0 if $n <= 0;

  # Use GMP if we have no XS or if size is small
  if ($n < 100000 || !Math::Prime::Util::prime_get_config()->{'xs'}) {
    if ($Math::Prime::Util::_GMPfunc{"ramanujan_tau"}) {
      return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::ramanujan_tau($n));
    }
  }

  # _taup is faster for small numbers, but gets very slow.  It's not a huge
  # deal, and the GMP code will probably get run for small inputs anyway.
  vecprod(map { _taupower($_->[0],$_->[1]) } Math::Prime::Util::factor_exp($n));
}


sub ExponentialIntegral {
  my($x) = @_;
  return - MPU_INFINITY if $x == 0;
  return 0              if $x == - MPU_INFINITY;
  return MPU_INFINITY   if $x == MPU_INFINITY;

  # Gotcha -- MPFR decided to make negative inputs return NaN.  Grrr.
  if ($x > 0 && _MPFR_available()) {
    my($wantbf,$xdigits) = _bfdigits($x);
    my $rnd = 0;  # MPFR_RNDN;
    my $bit_precision = int($xdigits * 3.322) + 4;
    my $rx = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
    Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
    my $eix = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($eix, $bit_precision);
    Math::MPFR::Rmpfr_eint($eix, $rx, $rnd);
    my $strval = Math::MPFR::Rmpfr_get_str($eix, 10, 0, $rnd);
    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + $strval;
  }

  $x = Math::BigFloat->new("$x") if defined $bignum::VERSION && ref($x) ne 'Math::BigFloat';

  my $tol = 1e-16;
  my $sum = 0.0;
  my($y, $t);
  my $c = 0.0;
  my $val; # The result from one of the four methods

  if ($x < -1) {
    # Continued fraction
    my $lc = 0;
    my $ld = 1 / (1 - $x);
    $val = $ld * (-exp($x));
    for my $n (1 .. 100000) {
      $lc = 1 / (2*$n + 1 - $x - $n*$n*$lc);
      $ld = 1 / (2*$n + 1 - $x - $n*$n*$ld);
      my $old = $val;
      $val *= $ld/$lc;
      last if abs($val - $old) <= ($tol * abs($val));
    }
  } elsif ($x < 0) {
    # Rational Chebyshev approximation
    my @C6p = ( -148151.02102575750838086,
                 150260.59476436982420737,
                  89904.972007457256553251,
                  15924.175980637303639884,
                   2150.0672908092918123209,
                    116.69552669734461083368,
                      5.0196785185439843791020);
    my @C6q = (  256664.93484897117319268,
                 184340.70063353677359298,
                  52440.529172056355429883,
                   8125.8035174768735759866,
                    750.43163907103936624165,
                     40.205465640027706061433,
                      1.0000000000000000000000);
    my $sumn = $C6p[0]-$x*($C6p[1]-$x*($C6p[2]-$x*($C6p[3]-$x*($C6p[4]-$x*($C6p[5]-$x*$C6p[6])))));
    my $sumd = $C6q[0]-$x*($C6q[1]-$x*($C6q[2]-$x*($C6q[3]-$x*($C6q[4]-$x*($C6q[5]-$x*$C6q[6])))));
    $val = log(-$x) - ($sumn / $sumd);
  } elsif ($x < -log($tol)) {
    # Convergent series
    my $fact_n = 1;
    $y = CONST_EULER-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
    $y = log($x)-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
    for my $n (1 .. 200) {
      $fact_n *= $x/$n;
      my $term = $fact_n / $n;
      $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
      last if $term < $tol;
    }
    $val = $sum;
  } else {
    # Asymptotic divergent series
    my $invx = 1.0 / $x;
    my $term = $invx;
    $sum = 1.0 + $term;
    for my $n (2 .. 200) {
      my $last_term = $term;
      $term *= $n * $invx;
      last if $term < $tol;
      if ($term < $last_term) {
        $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
      } else {
        $y = (-$last_term/3)-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
        last;
      }
    }
    $val = exp($x) * $invx * $sum;
  }
  $val;
}

sub LogarithmicIntegral {
  my($x,$opt) = @_;
  return 0              if $x == 0;
  return - MPU_INFINITY if $x == 1;
  return MPU_INFINITY   if $x == MPU_INFINITY;
  croak "Invalid input to LogarithmicIntegral:  x must be > 0" if $x <= 0;
  $opt = 0 unless defined $opt;

  # Remember MPFR eint doesn't handle negative inputs
  if ($x >= 1 && _MPFR_available()) {
    my $wantbf = 0;
    my $xdigits = 18;
    if ($opt) {
      $wantbf = length($x);
      $xdigits = $wantbf;
    } elsif (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
      $wantbf = _find_big_acc($x);
      $xdigits = $wantbf;
    }
    $xdigits += length(int(log(0.0+"$x"))) + 1;
    my $rnd = 0;  # MPFR_RNDN;
    my $bit_precision = int($xdigits * 3.322) + 4;
    my $rx = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
    Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
    Math::MPFR::Rmpfr_log($rx, $rx, $rnd);
    my $lix = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($lix, $bit_precision);
    Math::MPFR::Rmpfr_eint($lix, $rx, $rnd);
    my $strval = Math::MPFR::Rmpfr_get_str($lix, 10, 0, $rnd);
    return ($wantbf)  ?  _upgrade_to_float($strval,$wantbf)  :  0.0 + $strval;
  }

  if ($x == 2) {
    my $li2const = (ref($x) eq 'Math::BigFloat') ? Math::BigFloat->new(CONST_LI2) : 0.0+CONST_LI2;
    return $li2const;
  }

  $x = _bigint_to_int($x) if ref($x) && !defined $bignum::VERSION && $x <= 1e16;
  $x = Math::BigFloat->new("$x") if defined $bignum::VERSION && ref($x) ne 'Math::BigFloat';
  $x = _upgrade_to_float($x) if ref($x) && ref($x) ne 'Math::BigFloat' && $x > 1e16;

  # Do divergent series here for big inputs.  Common for big pc approximations.
  # Why is this here?
  #   1) exp(log(x)) results in a lot of lost precision
  #   2) exp(x) with lots of precision turns out to be really slow, and in
  #      this case it was unnecessary.
  my $tol = 1e-16;
  my $xdigits = 0;
  my $finalacc = 0;
  if (ref($x) =~ /^Math::Big/) {
    $xdigits = _find_big_acc($x);
    my $xlen = length($x->copy->bfloor->bstr());
    $xdigits = $xlen if $xdigits < $xlen;
    $finalacc = $xdigits;
    $xdigits += length(int(log(0.0+"$x"))) + 1;
    $tol = Math::BigFloat->new(10)->bpow(-$xdigits);
    $x->accuracy($xdigits);
  }
  my $logx = $xdigits ? $x->copy->blog(undef,$xdigits) : log($x);

  if ($x > 1e16) {
    my $invx = ref($logx) ? Math::BigFloat->bone / $logx : 1.0/$logx;
    # n = 0  =>  0!/(logx)^0 = 1/1 = 1
    # n = 1  =>  1!/(logx)^1 = 1/logx
    my $term = $invx;
    my $sum = 1.0 + $term;
    for my $n (2 .. 200) {
      my $last_term = $term;
      $term *= $n * $invx;
      last if $term < $tol;
      if ($term < $last_term) {
        $sum += $term;
      } else {
        $sum -= ($last_term/3);
        last;
      }
      $term->bround($xdigits) if $xdigits;
    }
    my $val = $x * $invx * $sum;
    $val->accuracy($finalacc) if $xdigits;
    return $val;
  }
  # Convergent series.
  if ($x >= 1) {
    my $fact_n = 1.0;
    my $nfac = 1.0;
    my $sum  = 0.0;
    for my $n (1 .. 200) {
      $fact_n *= $logx/$n;
      my $term = $fact_n / $n;
      $sum += $term;
      last if $term < $tol;
      $term->bround($xdigits) if $xdigits;
    }
    my $eulerconst = (ref($x) eq 'Math::BigFloat') ? Math::BigFloat->new(CONST_EULER) : 0.0+CONST_EULER;
    my $val = $eulerconst + log($logx) + $sum;
    $val->accuracy($finalacc) if $xdigits;
    return $val;
  }

  ExponentialIntegral($logx);
}

# Riemann Zeta function for native integers.
my @_Riemann_Zeta_Table = (
  0.6449340668482264364724151666460251892,  # zeta(2) - 1
  0.2020569031595942853997381615114499908,
  0.0823232337111381915160036965411679028,
  0.0369277551433699263313654864570341681,
  0.0173430619844491397145179297909205279,
  0.0083492773819228268397975498497967596,
  0.0040773561979443393786852385086524653,
  0.0020083928260822144178527692324120605,
  0.0009945751278180853371459589003190170,
  0.0004941886041194645587022825264699365,
  0.0002460865533080482986379980477396710,
  0.0001227133475784891467518365263573957,
  0.0000612481350587048292585451051353337,
  0.0000305882363070204935517285106450626,
  0.0000152822594086518717325714876367220,
  0.0000076371976378997622736002935630292,
  0.0000038172932649998398564616446219397,
  0.0000019082127165539389256569577951013,
  0.0000009539620338727961131520386834493,
  0.0000004769329867878064631167196043730,
  0.0000002384505027277329900036481867530,
  0.0000001192199259653110730677887188823,
  0.0000000596081890512594796124402079358,
  0.0000000298035035146522801860637050694,
  0.0000000149015548283650412346585066307,
  0.0000000074507117898354294919810041706,
  0.0000000037253340247884570548192040184,
  0.0000000018626597235130490064039099454,
  0.0000000009313274324196681828717647350,
  0.0000000004656629065033784072989233251,
  0.0000000002328311833676505492001455976,
  0.0000000001164155017270051977592973835,
  0.0000000000582077208790270088924368599,
  0.0000000000291038504449709968692942523,
  0.0000000000145519218910419842359296322,
  0.0000000000072759598350574810145208690,
  0.0000000000036379795473786511902372363,
  0.0000000000018189896503070659475848321,
  0.0000000000009094947840263889282533118,
);


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

  my $ix = ($x == int($x))  ?  "" . Math::BigInt->new($x)  :  0;

  # Try MPFR
  if (_MPFR_available()) {
    my($wantbf,$xdigits) = _bfdigits($x);
    my $rnd = 0;  # MPFR_RNDN;
    my $bit_precision = int($xdigits * 3.322) + 8;
    # Add more bits to account for the leading zeros.
    my $extra_bits = int((int(abs($x)/3)-1) * 3.322 + 0.5);

    my $zetax = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($zetax, $bit_precision + $extra_bits);

    if ($ix) {
      Math::MPFR::Rmpfr_zeta_ui($zetax, $ix, $rnd);
    } else {
      my $rx = Math::MPFR->new();
      Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
      Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
      Math::MPFR::Rmpfr_zeta($zetax, $rx, $rnd);
    }
    Math::MPFR::Rmpfr_sub_ui($zetax, $zetax, 1, $rnd);
    my $strval = Math::MPFR::Rmpfr_get_str($zetax, 10, $xdigits, $rnd);
    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + $strval;
  }

  # Try our GMP code if possible.
  if ($Math::Prime::Util::_GMPfunc{"zeta"}) {
    my($wantbf,$xdigits) = _bfdigits($x);
    # If we knew the *exact* number of zero digits, we could let GMP zeta
    # handle the correct rounding.  But we don't, so we have to go over.
    my $zero_dig = "".int($x / 3) - 1;
    my $strval = Math::Prime::Util::GMP::zeta($x, $xdigits + 8 + $zero_dig);
    if ($strval =~ s/^(1\.0*)/./) {
      $strval .= "e-".(length($1)-2) if length($1) > 2;
    } else {
      $strval =~ s/^(\d+)/$1-1/e;
    }

    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  : 0.0 + $strval;
  }

  # If we need a bigfloat result, then call our PP routine.
  if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
    require Math::Prime::Util::ZetaBigFloat;
    return Math::Prime::Util::ZetaBigFloat::RiemannZeta($x);
  }

  # No MPFR, no BigFloat.
  return 0.0 + $_Riemann_Zeta_Table[int($x)-2]
    if $x == int($x) && defined $_Riemann_Zeta_Table[int($x)-2];
  my $tol = 1.11e-16;

  # Series based on (2n)! / B_2n.
  # This is a simplification of the Cephes zeta function.
  my @A = (
      12.0,
     -720.0,
      30240.0,
     -1209600.0,
      47900160.0,
     -1892437580.3183791606367583212735166426,
      74724249600.0,
     -2950130727918.1642244954382084600497650,
      116467828143500.67248729113000661089202,
     -4597978722407472.6105457273596737891657,
      181521054019435467.73425331153534235290,
     -7166165256175667011.3346447367083352776,
      282908877253042996618.18640556532523927,
  );
  my $s = 0.0;
  my $rb = 0.0;
  foreach my $i (2 .. 10) {
    $rb = $i ** -$x;
    $s += $rb;
    return $s if abs($rb/$s) < $tol;
  }
  my $w = 10.0;
  $s = $s  +  $rb*$w/($x-1.0)  -  0.5*$rb;
  my $ra = 1.0;
  foreach my $i (0 .. 12) {
    my $k = 2*$i;
    $ra *= $x + $k;
    $rb /= $w;
    my $t = $ra*$rb/$A[$i];
    $s += $t;
    $t = abs($t/$s);
    last if $t < $tol;
    $ra *= $x + $k + 1.0;
    $rb /= $w;
  }
  return $s;
}

# Riemann R function
sub RiemannR {
  my($x) = @_;

  croak "Invalid input to ReimannR:  x must be > 0" if $x <= 0;

  # Use MPFR if possible.
  if (_MPFR_available()) {
    my($wantbf,$xdigits) = _bfdigits($x);
    my $rnd = 0;  # MPFR_RNDN;
    my $bit_precision = int($xdigits * 3.322) + 8;  # Add some extra

    my $rlogx = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rlogx, $bit_precision);
    Math::MPFR::Rmpfr_set_str($rlogx, "$x", 10, $rnd);
    Math::MPFR::Rmpfr_log($rlogx, $rlogx, $rnd);

    my $rpart_term = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rpart_term, $bit_precision);
    Math::MPFR::Rmpfr_set_str($rpart_term, "1", 10, $rnd);

    my $rzeta = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rzeta, $bit_precision);
    my $rterm = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rterm, $bit_precision);

    my $rsum = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rsum, $bit_precision);
    Math::MPFR::Rmpfr_set_str($rsum, "1", 10, $rnd);

    my $rstop = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($rstop, $bit_precision);
    Math::MPFR::Rmpfr_set_str($rstop, "1e-$xdigits", 10, $rnd);

    for my $k (1 .. 100000) {
      Math::MPFR::Rmpfr_mul($rpart_term, $rpart_term, $rlogx, $rnd);
      Math::MPFR::Rmpfr_div_ui($rpart_term, $rpart_term, $k, $rnd);

      Math::MPFR::Rmpfr_zeta_ui($rzeta, $k+1, $rnd);
      Math::MPFR::Rmpfr_sub_ui($rzeta, $rzeta, 1, $rnd);
      Math::MPFR::Rmpfr_mul_ui($rzeta, $rzeta, $k, $rnd);
      Math::MPFR::Rmpfr_add_ui($rzeta, $rzeta, $k, $rnd);
      Math::MPFR::Rmpfr_div($rterm, $rpart_term, $rzeta, $rnd);

      last if Math::MPFR::Rmpfr_less_p($rterm, $rstop);
      Math::MPFR::Rmpfr_add($rsum, $rsum, $rterm, $rnd);
    }
    my $strval = Math::MPFR::Rmpfr_get_str($rsum, 10, $xdigits, $rnd);
    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + $strval;
  }

  if ($Math::Prime::Util::_GMPfunc{"riemannr"}) {
    my($wantbf,$xdigits) = _bfdigits($x);
    my $strval = Math::Prime::Util::GMP::riemannr($x, $xdigits);
    return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + $strval;
  }

# TODO: look into this as a generic solution
if (0 && $Math::Prime::Util::_GMPfunc{"zeta"}) {
  my($wantbf,$xdigits) = _bfdigits($x);
  $x = _upgrade_to_float($x);

  my $extra_acc = 4;
  $xdigits += $extra_acc;
  $x->accuracy($xdigits);

  my $logx = log($x);
  my $part_term = $x->copy->bone;
  my $sum = $x->copy->bone;
  my $tol = $x->copy->bone->brsft($xdigits-1, 10);
  my $bigk = $x->copy->bone;
  my $term;
  for my $k (1 .. 10000) {
    $part_term *= $logx / $bigk;
    my $zarg = $bigk->copy->binc;
    my $zeta = (RiemannZeta($zarg) * $bigk) + $bigk;
    #my $strval = Math::Prime::Util::GMP::zeta($k+1, $xdigits + int(($k+1) / 3));
    #my $zeta = Math::BigFloat->new($strval)->bdec->bmul($bigk)->badd($bigk);
    $term = $part_term / $zeta;
    $sum += $term;
    last if $term < ($tol * $sum);
    $bigk->binc;
  }
  $sum->bround($xdigits-$extra_acc);
  my $strval = "$sum";
  return ($wantbf)  ?  Math::BigFloat->new($strval,$wantbf)  :  0.0 + $strval;
}

  if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
    require Math::Prime::Util::ZetaBigFloat;
    return Math::Prime::Util::ZetaBigFloat::RiemannR($x);
  }

  my $sum = 0.0;
  my $tol = 1e-18;
  my($c, $y, $t) = (0.0);
  if ($x > 10**17) {
    my @mob = Math::Prime::Util::moebius(0,300);
    for my $k (1 .. 300) {
      next if $mob[$k] == 0;
      my $term = $mob[$k] / $k *
                 Math::Prime::Util::LogarithmicIntegral($x**(1.0/$k));
      $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
      last if abs($term) < ($tol * abs($sum));
    }
  } else {
    $y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
    my $flogx = log($x);
    my $part_term = 1.0;
    for my $k (1 .. 10000) {
      my $zeta = ($k <= $#_Riemann_Zeta_Table)
                 ? $_Riemann_Zeta_Table[$k+1-2]    # Small k from table
                 : RiemannZeta($k+1);              # Large k from function
      $part_term *= $flogx / $k;
      my $term = $part_term / ($k + $k * $zeta);
      $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
      last if $term < ($tol * $sum);
    }
  }
  return $sum;
}

sub LambertW {
  my $x = shift;
  croak "Invalid input to LambertW:  x must be >= -1/e" if $x < -0.36787944118;
  $x = _upgrade_to_float($x) if ref($x) eq 'Math::BigInt';
  my $xacc = ref($x) ? _find_big_acc($x) : 0;
  my $w;

  if ($Math::Prime::Util::_GMPfunc{"lambertw"}) {
    my $w = (!$xacc)
          ? 0.0 + Math::Prime::Util::GMP::lambertw($x)
          : $x->copy->bzero->badd(Math::Prime::Util::GMP::lambertw($x, $xacc));
    return $w;
  }

  # Approximation
  if ($x < -0.06) {
    my $ti = $x * 2 * exp($x-$x+1) + 2;
    return -1 if $ti <= 0;
    my $t  = sqrt($ti);
    $w = (-1 + 1/6*$t + (257/720)*$t*$t + (13/720)*$t*$t*$t) / (1 + (5/6)*$t + (103/720)*$t*$t);
  } elsif ($x < 1.363) {
    my $l1 = log($x + 1);
    $w = $l1 * (1 - log(1+$l1) / (2+$l1));
  } elsif ($x < 3.7) {
    my $l1 = log($x);
    my $l2 = log($l1);
    $w = $l1 - $l2 - log(1 - $l2/$l1)/2.0;
  } else {
    my $l1 = log($x);
    my $l2 = log($l1);
    my $d1 = 2 * $l1 * $l1;
    my $d2 = 3 * $l1 * $d1;
    my $d3 = 2 * $l1 * $d2;
    my $d4 = 5 * $l1 * $d3;
    $w = $l1 - $l2 + $l2/$l1 + $l2*($l2-2)/$d1
       + $l2*(6+$l2*(-9+2*$l2))/$d2
       + $l2*(-12+$l2*(36+$l2*(-22+3*$l2)))/$d3
       + $l2*(60+$l2*(-300+$l2*(350+$l2*(-125+12*$l2))))/$d4;
  }

  # Now iterate to get the answer
  #
  # Newton:
  #   $w = $w*(log($x) - log($w) + 1) / ($w+1);
  # Halley:
  #   my $e = exp($w);
  #   my $f = $w * $e - $x;
  #   $w -= $f / ($w*$e+$e - ($w+2)*$f/(2*$w+2));

  # Fritsch converges quadratically, so tolerance could be 4x smaller.  Use 2x.
  my $tol = ($xacc) ? 10**(-int(1+$xacc/2)) : 1e-16;
  $w->accuracy($xacc+10) if $xacc;
  for (1 .. 200) {
    last if $w == 0;
    my $w1 = $w + 1;
    my $zn = log($x/$w) - $w;
    my $qn = $w1 * 2 * ($w1+(2*$zn/3));
    my $en = ($zn/$w1) * ($qn-$zn)/($qn-$zn*2);
    my $wen = $w * $en;
    $w += $wen;
    last if abs($wen) < $tol;
  }
  $w->accuracy($xacc) if $xacc;

  $w;
}

my $_Pi = "3.14159265358979323846264338328";
sub Pi {
  my $digits = shift;
  return 0.0+$_Pi unless $digits;
  return 0.0+sprintf("%.*lf", $digits-1, $_Pi) if $digits < 15;
  return _upgrade_to_float($_Pi, $digits) if $digits < 30;

  # Performance ranking:
  #   MPFR             The first two are fastest by a wide margin
  #   MPU::GMP         Both use AGM.  MPFR is very slightly faster.
  #   Perl AGM w/GMP   also AGM, nice growth rate, but slower than above
  #   C pidigits       much worse than above, but faster than the others
  #   Perl AGM         without Math::BigInt::GMP, it's sluggish
  #   Math::BigFloat   much slower than AGM
  #
  # With a few thousand digits, any of the top 4 are fine.
  # At 10k digits, the first two are pulling away.
  # At 50k digits, the first three are 5-20x faster than C pidigits, and
  #   pray you're not having to the Perl BigFloat methods without GMP.
  # At 100k digits, the first two are 15x faster than the third, C pidigits
  #   is 200x slower, and the rest thousands of times slower.
  # At 1M digits, the first two are under 2 seconds, the third is over a
  #   minute, and C pixigits at 1.5 hours.
  #
  # Interestingly, Math::BigInt::Pari, while greatly faster than Calc, is
  # *much* slower than GMP for these operations (both AGM and Machin).  While
  # Perl AGM with the Math::BigInt::GMP backend will pull away from C pidigits,
  # using it with the other backends doesn't do so.
  #
  # The GMP program at https://gmplib.org/download/misc/gmp-chudnovsky.c
  # will run ~4x faster than the MPFR code.

  my $have_bigint_gmp = Math::BigInt->config()->{lib} =~ /GMP/;
  my $have_xdigits    = Math::Prime::Util::prime_get_config()->{'xs'};
  my $_verbose = Math::Prime::Util::prime_get_config()->{'verbose'};

  # Uses AGM to get performance almost as good as MPFR
  if ($Math::Prime::Util::_GMPfunc{"Pi"}) {
    print "  using MPUGMP for Pi($digits)\n" if $_verbose;
    return _upgrade_to_float( Math::Prime::Util::GMP::Pi($digits) );
  }

  # MPFR is a bit faster than MPU-GMP's AGM.  Both are much faster than others.
  if ( (!$have_xdigits || $digits > 60) && _MPFR_available()) {
    print "  using MPFR for Pi($digits)\n" if $_verbose;
    my $rnd = 0;  # MPFR_RNDN;
    my $bit_precision = int($digits * 3.322) + 40;
    my $pi = Math::MPFR->new();
    Math::MPFR::Rmpfr_set_prec($pi, $bit_precision);
    Math::MPFR::Rmpfr_const_pi($pi, $rnd);
    my $strval = Math::MPFR::Rmpfr_get_str($pi, 10, $digits, $rnd);
    return _upgrade_to_float($strval);
  }

  # We could consider looking for Pari

  # This has a *much* better growth rate than the later solutions.
  if ( !$have_xdigits || ($have_bigint_gmp && $digits > 100) ) {
    print "  using Perl AGM for Pi($digits)\n" if $_verbose;
    # Brent-Salamin (aka AGM or Gauss-Legendre)
    $digits += 8;
    my $HALF = _upgrade_to_float(0.5);
    my ($an, $bn, $tn, $pn) = ($HALF->copy->bone, $HALF->copy->bsqrt($digits),
                               $HALF->copy->bmul($HALF), $HALF->copy->bone);
    while ($pn < $digits) {
      my $prev_an = $an->copy;
      $an->badd($bn)->bmul($HALF, $digits);
      $bn->bmul($prev_an)->bsqrt($digits);
      $prev_an->bsub($an);
      $tn->bsub($pn * $prev_an * $prev_an);
      $pn->badd($pn);
    }
    $an->badd($bn);
    $an->bmul($an,$digits)->bdiv(4*$tn, $digits-8);
    return $an;
  }

  # Spigot method in C.  Low overhead but not good growth rate.
  if ($have_xdigits) {
    print "  using XS spigot for Pi($digits)\n" if $_verbose;
    return _upgrade_to_float(Math::Prime::Util::_pidigits($digits));
  }

  # We're going to have to use the Math::BigFloat code.
  # 1) it rounds incorrectly (e.g. 761, 1372, 1509,...).
  #    Fix by adding some digits and rounding.
  # 2) AGM is *much* faster once past ~2000 digits
  # 3) It is very slow without the GMP backend.  The Pari backend helps
  #    but it still pretty bad.  With Calc it's glacial for large inputs.

  #           Math::BigFloat                AGM              spigot   AGM
  # Size     GMP    Pari  Calc        GMP    Pari  Calc        C      C+GMP
  #   500   0.04    0.60   0.30      0.08    0.10   0.47      0.09    0.06
  #  1000   0.04    0.11   1.82      0.09    0.14   1.82      0.09    0.06
  #  2000   0.07    0.37  13.5       0.09    0.34   9.16      0.10    0.06
  #  4000   0.14    2.17 107.8       0.12    1.14  39.7       0.20    0.06
  #  8000   0.52   15.7              0.22    4.63 186.2       0.56    0.08
  # 16000   2.73  121.8              0.52   19.2              2.00    0.08
  # 32000  15.4                      1.42                     7.78    0.12
  #                                   ^                        ^       ^
  #                                   |      use this THIRD ---+       |
  #                use this SECOND ---+                                |
  #                                                  use this FIRST ---+
  # approx
  # growth  5.6x    7.6x   8.0x      2.7x    4.1x   4.7x      3.9x    2.0x

  print "  using BigFloat for Pi($digits)\n" if $_verbose;
  _upgrade_to_float(0);
  return Math::BigFloat::bpi($digits+10)->round($digits);
}

sub forpart {
  my($sub, $n, $rhash) = @_;
  _forcompositions(1, $sub, $n, $rhash);
}
sub forcomp {
  my($sub, $n, $rhash) = @_;
  _forcompositions(0, $sub, $n, $rhash);
}
sub _forcompositions {
  my($ispart, $sub, $n, $rhash) = @_;
  _validate_positive_integer($n);
  my($mina, $maxa, $minn, $maxn, $primeq) = (1,$n,1,$n,-1);
  if (defined $rhash) {
    croak "forpart second argument must be a hash reference"
      unless ref($rhash) eq 'HASH';
    if (defined $rhash->{amin}) {
      $mina = $rhash->{amin};
      _validate_positive_integer($mina);
    }
    if (defined $rhash->{amax}) {
      $maxa = $rhash->{amax};
      _validate_positive_integer($maxa);
    }
    $minn = $maxn = $rhash->{n} if defined $rhash->{n};
    $minn = $rhash->{nmin} if defined $rhash->{nmin};
    $maxn = $rhash->{nmax} if defined $rhash->{nmax};
    _validate_positive_integer($minn);
    _validate_positive_integer($maxn);
    if (defined $rhash->{prime}) {
      $primeq = $rhash->{prime};
      _validate_positive_integer($primeq);
    }
   $mina = 1 if $mina < 1;
   $maxa = $n if $maxa > $n;
   $minn = 1 if $minn < 1;
   $maxn = $n if $maxn > $n;
   $primeq = 2 if $primeq != -1 && $primeq != 0;
  }

  $sub->() if $n == 0 && $minn <= 1;
  return if $n < $minn || $minn > $maxn || $mina > $maxa || $maxn <= 0 || $maxa <= 0;

  my $oldforexit = Math::Prime::Util::_start_for_loop();
  my ($x, $y, $r, $k);
  my @a = (0) x ($n);
  $k = 1;
  $a[0] = $mina - 1;
  $a[1] = $n - $mina + 1;
  while ($k != 0) {
    $x = $a[$k-1]+1;
    $y = $a[$k]-1;
    $k--;
    $r = $ispart ? $x : 1;
    while ($r <= $y) {
      $a[$k] = $x;
      $x = $r;
      $y -= $x;
      $k++;
    }
    $a[$k] = $x + $y;
    # Restrict size
    while ($k+1 > $maxn) {
      $a[$k-1] += $a[$k];
      $k--;
    }
    next if $k+1 < $minn;
    # Restrict values
    if ($mina > 1 || $maxa < $n) {
      last if $a[0] > $maxa;
      if ($ispart) {
        next if $a[$k] > $maxa;
      } else {
        next if Math::Prime::Util::vecany(sub{ $_ < $mina || $_ > $maxa }, @a[0..$k]);
      }
    }
    next if $primeq == 0 && Math::Prime::Util::vecany(sub{ is_prime($_) }, @a[0..$k]);
    next if $primeq == 2 && Math::Prime::Util::vecany(sub{ !is_prime($_) }, @a[0..$k]);
    last if Math::Prime::Util::_get_forexit();
    $sub->(@a[0 .. $k]);
  }
  Math::Prime::Util::_end_for_loop($oldforexit);
}
sub forcomb {
  my($sub, $n, $k) = @_;
  _validate_positive_integer($n);

  my($begk, $endk);
  if (defined $k) {
    _validate_positive_integer($k);
    return if $k > $n;
    $begk = $endk = $k;
  } else {
    $begk = 0;
    $endk = $n;
  }

  my $oldforexit = Math::Prime::Util::_start_for_loop();
  for my $k ($begk .. $endk) {
    if ($k == 0) {
      $sub->();
    } else {
      my @c = 0 .. $k-1;
      while (1) {
        $sub->(@c);
        last if Math::Prime::Util::_get_forexit();
        next if $c[-1]++ < $n-1;
        my $i = $k-2;
        $i-- while $i >= 0 && $c[$i] >= $n-($k-$i);
        last if $i < 0;
        $c[$i]++;
        while (++$i < $k) { $c[$i] = $c[$i-1] + 1; }
      }
    }
    last if Math::Prime::Util::_get_forexit();
  }
  Math::Prime::Util::_end_for_loop($oldforexit);
}
sub _forperm {
  my($sub, $n, $all_perm) = @_;
  my $k = $n;
  my @c = reverse 0 .. $k-1;
  my $inc = 0;
  my $send = 1;
  my $oldforexit = Math::Prime::Util::_start_for_loop();
  while (1) {
    if (!$all_perm) {   # Derangements via simple filtering.
      $send = 1;
      for my $p (0 .. $#c) {
        if ($c[$p] == $k-$p-1) {
          $send = 0;
          last;
        }
      }
    }
    if ($send) {
      $sub->(reverse @c);
      last if Math::Prime::Util::_get_forexit();
    }
    if (++$inc & 1) {
      @c[0,1] = @c[1,0];
      next;
    }
    my $j = 2;
    $j++ while $j < $k && $c[$j] > $c[$j-1];
    last if $j >= $k;
    my $m = 0;
    $m++ while $c[$j] > $c[$m];
    @c[$j,$m] = @c[$m,$j];
    @c[0..$j-1] = reverse @c[0..$j-1];
  }
  Math::Prime::Util::_end_for_loop($oldforexit);
}
sub forperm {
  my($sub, $n, $k) = @_;
  _validate_positive_integer($n);
  croak "Too many arguments for forperm" if defined $k;
  return $sub->() if $n == 0;
  return $sub->(0) if $n == 1;
  _forperm($sub, $n, 1);
}
sub forderange {
  my($sub, $n, $k) = @_;
  _validate_positive_integer($n);
  croak "Too many arguments for forderange" if defined $k;
  return $sub->() if $n == 0;
  return if $n == 1;
  _forperm($sub, $n, 0);
}

sub _multiset_permutations {
  my($sub, $prefix, $ar, $sum) = @_;

  return if $sum == 0;

  # Remove any values with 0 occurances
  my @n = grep { $_->[1] > 0 } @$ar;

  if ($sum == 1) {                       # A single value
    $sub->(@$prefix, $n[0]->[0]);
  } elsif ($sum == 2) {                  # Optimize the leaf case
    my($n0,$n1) = map { $_->[0] } @n;
    if (@n == 1) {
      $sub->(@$prefix, $n0, $n0);
    } else {
      $sub->(@$prefix, $n0, $n1);
      $sub->(@$prefix, $n1, $n0) unless Math::Prime::Util::_get_forexit();
    }
  } elsif (0 && $sum == scalar(@n)) {         # All entries have 1 occurance
    # TODO:  Figure out a way to use this safely.  We need to capture any
    #        lastfor that was seen in the forperm.
    my @i = map { $_->[0] } @n;
    Math::Prime::Util::forperm(sub { $sub->(@$prefix, @i[@_]) }, 1+$#i);
  } else {                               # Recurse over each leading value
    for my $v (@n) {
      $v->[1]--;
      push @$prefix, $v->[0];
      no warnings 'recursion';
      _multiset_permutations($sub, $prefix, \@n, $sum-1);
      pop @$prefix;
      $v->[1]++;
      last if Math::Prime::Util::_get_forexit();
    }
  }
}

sub numtoperm {
  my($n,$k) = @_;
  _validate_positive_integer($n);
  _validate_positive_integer($k);
  return () if $n == 0;
  return (0) if $n == 1;
  my $f = factorial($n-1);
  $k %= vecprod($f,$n) if int($k/$f) >= $n;
  my @S = map { $_ } 0 .. $n-1;
  my @V;
  while ($n-- > 0) {
    my $i = int($k/$f);
    push @V, splice(@S,$i,1);
    last if $n == 0;
    $k -= $i*$f;
    $f /= $n;
  }
  @V;
}

sub permtonum {
  my $A = shift;
  croak "permtonum argument must be an array reference"
    unless ref($A) eq 'ARRAY';
  my $n = scalar(@$A);
  return 0 if $n == 0;
  {
    my %S;
    for my $v (@$A) {
      croak "permtonum invalid permutation array"
        if !defined $v || $v < 0 || $v >= $n || $S{$v}++;
    }
  }
  my $f = factorial($n-1);
  my $rank = 0;
  for my $i (0 .. $n-2) {
    my $k = 0;
    for my $j ($i+1 .. $n-1) {
      $k++ if $A->[$j] < $A->[$i];
    }
    $rank = Math::Prime::Util::vecsum($rank, Math::Prime::Util::vecprod($k,$f));
    $f /= $n-$i-1;
  }
  $rank;
}

sub randperm {
  my($n,$k) = @_;
  _validate_positive_integer($n);
  if (defined $k) {
    _validate_positive_integer($k);
  }
  $k = $n if !defined($k) || $k > $n;
  return () if $k == 0;

  my @S;
  if ("$k"/"$n" <= 0.30) {
    my %seen;
    my $v;
    for my $i (1 .. $k) {
      do { $v = Math::Prime::Util::urandomm($n); } while $seen{$v}++;
      push @S,$v;
    }
  } else {
    @S = map { $_ } 0..$n-1;
    for my $i (0 .. $n-2) {
      last if $i >= $k;
      my $j = Math::Prime::Util::urandomm($n-$i);
      @S[$i,$i+$j] = @S[$i+$j,$i];
    }
    $#S = $k-1;
  }
  return @S;
}

sub shuffle {
  my @S=@_;
  # Note: almost all the time is spent in urandomm.
  for (my $i = $#S; $i >= 1; $i--) {
    my $j = Math::Prime::Util::urandomm($i+1);
    @S[$i,$j] = @S[$j,$i];
  }
  @S;
}

###############################################################################
#       Random numbers
###############################################################################

# PPFE:  irand irand64 drand random_bytes csrand srand _is_csprng_well_seeded
sub urandomb {
  my($n) = @_;
  return 0 if $n <= 0;
  return ( Math::Prime::Util::irand() >> (32-$n) ) if $n <= 32;
  return ( Math::Prime::Util::irand64() >> (64-$n) ) if MPU_MAXBITS >= 64 && $n <= 64;
  my $bytes = Math::Prime::Util::random_bytes(($n+7)>>3);
  my $binary = substr(unpack("B*",$bytes),0,$n);
  return Math::BigInt->new("0b$binary");
}
sub urandomm {
  my($n) = @_;
  # _validate_positive_integer($n);
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::urandomm($n))
    if $Math::Prime::Util::_GMPfunc{"urandomm"};
  return 0 if $n <= 1;
  my $r;
  if ($n <= 4294967295) {
    my $rmax = int(4294967295 / $n) * $n;
    do { $r = Math::Prime::Util::irand() } while $r >= $rmax;
  } elsif (!ref($n)) {
    my $rmax = int(~0 / $n) * $n;
    do { $r = Math::Prime::Util::irand64() } while $r >= $rmax;
  } else {
    # TODO: verify and try to optimize this
    my $bits = length($n->as_bin) - 2;
    my $bytes = 1 + (($bits+7)>>3);
    my $rmax = Math::BigInt->bone->blsft($bytes*8)->bdec;
    my $overflow = $rmax - ($rmax % $n);
    do { $r = Math::Prime::Util::urandomb($bytes*8); } while $r >= $overflow;
  }
  return $r % $n;
}

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

  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::random_prime($low, $high))
    if $Math::Prime::Util::_GMPfunc{"random_prime"};

  require Math::Prime::Util::RandomPrimes;
  return Math::Prime::Util::RandomPrimes::random_prime($low,$high);
}

sub random_ndigit_prime {
  my($digits) = @_;
  _validate_positive_integer($digits, 1);
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::random_ndigit_prime($digits))
    if $Math::Prime::Util::_GMPfunc{"random_ndigit_prime"};
  require Math::Prime::Util::RandomPrimes;
  return Math::Prime::Util::RandomPrimes::random_ndigit_prime($digits);
}
sub random_nbit_prime {
  my($bits) = @_;
  _validate_positive_integer($bits, 2);
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::random_nbit_prime($bits))
    if $Math::Prime::Util::_GMPfunc{"random_nbit_prime"};
  require Math::Prime::Util::RandomPrimes;
  return Math::Prime::Util::RandomPrimes::random_nbit_prime($bits);
}
sub random_strong_prime {
  my($bits) = @_;
  _validate_positive_integer($bits, 128);
  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::random_strong_prime($bits))
    if $Math::Prime::Util::_GMPfunc{"random_strong_prime"};
  require Math::Prime::Util::RandomPrimes;
  return Math::Prime::Util::RandomPrimes::random_strong_prime($bits);
}

sub random_maurer_prime {
  my($bits) = @_;
  _validate_positive_integer($bits, 2);

  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::random_maurer_prime($bits))
    if $Math::Prime::Util::_GMPfunc{"random_maurer_prime"};

  require Math::Prime::Util::RandomPrimes;
  my ($n, $cert) = Math::Prime::Util::RandomPrimes::random_maurer_prime_with_cert($bits);
  croak "maurer prime $n failed certificate verification!"
        unless Math::Prime::Util::verify_prime($cert);

  return $n;
}

sub random_shawe_taylor_prime {
  my($bits) = @_;
  _validate_positive_integer($bits, 2);

  return Math::Prime::Util::_reftyped($_[0], Math::Prime::Util::GMP::random_shawe_taylor_prime($bits))
    if $Math::Prime::Util::_GMPfunc{"random_shawe_taylor_prime"};

  require Math::Prime::Util::RandomPrimes;
  my ($n, $cert) = Math::Prime::Util::RandomPrimes::random_shawe_taylor_prime_with_cert($bits);
  croak "shawe-taylor prime $n failed certificate verification!"
        unless Math::Prime::Util::verify_prime($cert);

  return $n;
}

sub miller_rabin_random {
  my($n, $k, $seed) = @_;
  _validate_positive_integer($n);
  if (scalar(@_) == 1 ) { $k = 1; } else { _validate_positive_integer($k); }

  return 1 if $k <= 0;

  if ($Math::Prime::Util::_GMPfunc{"miller_rabin_random"}) {
    return Math::Prime::Util::GMP::miller_rabin_random($n, $k, $seed) if defined $seed;
    return Math::Prime::Util::GMP::miller_rabin_random($n, $k);
  }

  # Math::Prime::Util::prime_get_config()->{'assume_rh'})  ==>  2*log(n)^2
  if ($k >= int(3*$n/4) ) {
    for (2 .. int(3*$n/4)+2) {
      return 0 unless Math::Prime::Util::is_strong_pseudoprime($n, $_);
    }
    return 1;
  }
  my $brange = $n-2;
  return 0 unless Math::Prime::Util::is_strong_pseudoprime($n, Math::Prime::Util::urandomm($brange)+2 );
  $k--;
  while ($k > 0) {
    my $nbases = ($k >= 20) ? 20 : $k;
    return 0 unless is_strong_pseudoprime($n, map { urandomm($brange)+2 } 1 .. $nbases);
    $k -= $nbases;
  }
  1;
}

sub random_semiprime {
  my($b) = @_;
  return 0 if defined $b && int($b) < 0;
  _validate_positive_integer($b,4);

  my $n;
  my $min = ($b <= MPU_MAXBITS)  ?  (1 << ($b-1))  :  BTWO->copy->bpow($b-1);
  my $max = $min + ($min - 1);
  my $L = $b >> 1;
  my $N = $b - $L;
  my $one = ($b <= MPU_MAXBITS) ? 1 : BONE;
  do {
    $n = $one * random_nbit_prime($L) * random_nbit_prime($N);
  } while $n < $min || $n > $max;
  $n = _bigint_to_int($n) if ref($n) && $n->bacmp(BMAX) <= 0;
  $n;
}

sub random_unrestricted_semiprime {
  my($b) = @_;
  return 0 if defined $b && int($b) < 0;
  _validate_positive_integer($b,3);

  my $n;
  my $min = ($b <= MPU_MAXBITS)  ?  (1 << ($b-1))  :  BTWO->copy->bpow($b-1);
  my $max = $min + ($min - 1);

  if ($b <= 64) {
    do {
      $n = $min + urandomb($b-1);
    } while !Math::Prime::Util::is_semiprime($n);
  } else {
    # Try to get probabilities right for small divisors
    my %M = (
      2 => 1.91218397452243,
      3 => 1.33954826555021,
      5 => 0.854756717114822,
      7 => 0.635492301836862,
      11 => 0.426616792046787,
      13 => 0.368193843118344,
      17 => 0.290512701603111,
      19 => 0.263359264658156,
      23 => 0.222406328935102,
      29 => 0.181229250520242,
      31 => 0.170874199059434,
      37 => 0.146112155735473,
      41 => 0.133427839963585,
      43 => 0.127929010905662,
      47 => 0.118254609086782,
      53 => 0.106316418106489,
      59 => 0.0966989675438643,
      61 => 0.0938833658008547,
      67 => 0.0864151823151671,
      71 => 0.0820822953188297,
      73 => 0.0800964416340746,
      79 => 0.0747060914833344,
      83 => 0.0714973706654851,
      89 => 0.0672115468436284,
      97 => 0.0622818892486191,
      101 => 0.0600855891549939,
      103 => 0.0590613570015407,
      107 => 0.0570921135626976,
      109 => 0.0561691667641485,
      113 => 0.0544330141081874,
      127 => 0.0490620204315701,
    );
    my ($p,$r);
    $r = Math::Prime::Util::drand();
    for my $prime (2..127) {
      next unless defined $M{$prime};
      my $PR = $M{$prime} / $b  +  0.19556 / $prime;
      if ($r <= $PR) {
        $p = $prime;
        last;
      }
      $r -= $PR;
    }
    if (!defined $p) {
      # Idea from Charles Greathouse IV, 2010.  The distribution is right
      # at the high level (small primes weighted more and not far off what
      # we get with the uniform selection), but there is a noticeable skew
      # toward primes with a large gap after them.  For instance 3 ends up
      # being weighted as much as 2, and 7 more than 5.
      #
      # Since we handled small divisors earlier, this is less bothersome.
      my $M = 0.26149721284764278375542683860869585905;
      my $weight = $M + log($b * log(2)/2);
      my $minr = log(log(131));
      do {
        $r  = Math::Prime::Util::drand($weight) - $M;
      } while $r < $minr;
      # Using Math::BigFloat::bexp is ungodly slow, so avoid at all costs.
      my $re = exp($r);
      my $a = ($re < log(~0)) ? int(exp($re)+0.5)
                              : _upgrade_to_float($re)->bexp->bround->as_int;
      $p = $a < 2 ? 2 : Math::Prime::Util::prev_prime($a+1);
    }
    my $ranmin = ref($min) ? $min->badd($p-1)->bdiv($p)->as_int : int(($min+$p-1)/$p);
    my $ranmax = ref($max) ? $max->bdiv($p)->as_int : int($max/$p);
    my $q = random_prime($ranmin, $ranmax);
    $n = Math::Prime::Util::vecprod($p,$q);
  }
  $n = _bigint_to_int($n) if ref($n) && $n->bacmp(BMAX) <= 0;
  $n;
}

1;

__END__