This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.37
in repository libmath-prime-util-perl.

commit f2d43c34330f108b3330e09a6688406bf7b35f0e
Author: Dana Jacobsen <d...@acm.org>
Date:   Tue Jan 21 15:40:27 2014 -0800

    Move some functions from Util.pm to util.c.  Faster, and reduces startup 
bloat.
---
 Changes                           |  14 +-
 XS.xs                             |  33 +++-
 bin/factor.pl                     |   5 +-
 examples/test-factor-gnufactor.pl |   2 +-
 lib/Math/Prime/Util.pm            | 337 +-------------------------------------
 lib/Math/Prime/Util/PP.pm         | 317 ++++++++++++++++++++++++++++++++++-
 lmo.c                             |   6 +-
 util.c                            | 161 ++++++++++++++++--
 util.h                            |   8 +-
 9 files changed, 522 insertions(+), 361 deletions(-)

diff --git a/Changes b/Changes
index d81e0b2..f0b0955 100644
--- a/Changes
+++ b/Changes
@@ -10,11 +10,21 @@ Revision history for Perl module Math::Prime::Util
     - Dynamically loads the PP code and Math::BigInt only when needed.  This
       removes a lot of bloat for the usual cases:
         2.0 MB  perl -E 'say 1'
-        4.2 MB  Math::Prime::XS
-        4.6 MB  MPU 0.37
+        4.4 MB  MPU 0.37
+        4.5 MB  Math::Prime::XS + Math::Factor::XS
+        5.3 MB  Math::Pari
         9.6 MB  MPU 0.36
        10.0 MB  MPU 0.35
 
+    - Together with the previous, adjusted factor script.  Over 2x faster
+      with native integer (non-expression) arguments.  This is just not
+      loading thousands of lines of Perl code that aren't used, which
+      was more time-consuming than the actual factoring.
+
+    - nth_prime_{lower,upper,approx} and prime_count_{lower,upper,approx}
+      moved to XS->PP.  This helps us slim down and cut startup overhead.
+
+
 0.36  2014-01-13
 
     [API Changes]
diff --git a/XS.xs b/XS.xs
index cb234f6..d148121 100644
--- a/XS.xs
+++ b/XS.xs
@@ -589,27 +589,46 @@ next_prime(IN SV* svn)
   ALIAS:
     prev_prime = 1
     nth_prime = 2
+    nth_prime_upper = 3
+    nth_prime_lower = 4
+    nth_prime_approx = 5
+    prime_count_upper = 6
+    prime_count_lower = 7
+    prime_count_approx = 8
   PPCODE:
     if (_validate_int(aTHX_ svn, 0)) {
       UV n = my_svuv(svn);
-      if ( ((ix == 0) && (n >= MPU_MAX_PRIME))    ||
-           ((ix == 2) && (n >= MPU_MAX_PRIME_IDX)) ) {
+      if ( (n >= MPU_MAX_PRIME     && ix == 0) ||
+           (n >= MPU_MAX_PRIME_IDX && (ix==2 || ix==3 || ix==4 || ix==5)) ) {
         /* Out of range.  Fall through to Perl. */
       } else {
         UV ret;
         switch (ix) {
           case 0: ret = next_prime(n);  break;
           case 1: ret = (n < 3) ? 0 : prev_prime(n);  break;
-          case 2:
-          default:ret = _XS_nth_prime(n);  break;
+          case 2: ret = nth_prime(n); break;
+          case 3: ret = nth_prime_upper(n); break;
+          case 4: ret = nth_prime_lower(n); break;
+          case 5: ret = nth_prime_approx(n); break;
+          case 6: ret = prime_count_upper(n); break;
+          case 7: ret = prime_count_lower(n); break;
+          case 8:
+          default:ret = prime_count_approx(n); break;
         }
         XSRETURN_UV(ret);
       }
     }
     switch (ix) {
-      case 0:  _vcallsub("_generic_next_prime");     break;
-      case 1:  _vcallsub("_generic_prev_prime");     break;
-      default: _vcallsub_with_pp("nth_prime");           break;
+      case 0:  _vcallsub("_generic_next_prime");        break;
+      case 1:  _vcallsub("_generic_prev_prime");        break;
+      case 2:  _vcallsub_with_pp("nth_prime");          break;
+      case 3:  _vcallsub_with_pp("nth_prime_upper");    break;
+      case 4:  _vcallsub_with_pp("nth_prime_lower");    break;
+      case 5:  _vcallsub_with_pp("nth_prime_approx");   break;
+      case 6:  _vcallsub_with_pp("prime_count_upper");  break;
+      case 7:  _vcallsub_with_pp("prime_count_lower");  break;
+      case 8:
+      default: _vcallsub_with_pp("prime_count_approx"); break;
     }
     return; /* skip implicit PUTBACK */
 
diff --git a/bin/factor.pl b/bin/factor.pl
index 6fca37b..779c2db 100755
--- a/bin/factor.pl
+++ b/bin/factor.pl
@@ -2,10 +2,8 @@
 use strict;
 use warnings;
 use Getopt::Long;
-use bigint try => 'GMP';
 use Math::Prime::Util qw/factor nth_prime prime_set_config/;
 $| = 1;
-no bigint;
 
 # Allow execution of any of these functions in the command line
 my @mpu_funcs = (qw/next_prime prev_prime prime_count nth_prime random_prime
@@ -23,7 +21,7 @@ GetOptions(\%opts,
           ) || die_usage();
 if (exists $opts{'version'}) {
   my $version_str =
-   "factor.pl version 1.1 using Math::Prime::Util $Math::Prime::Util::VERSION";
+   "factor.pl version 1.2 using Math::Prime::Util $Math::Prime::Util::VERSION";
   $version_str .= " and MPU::GMP $Math::Prime::Util::GMP::VERSION"
     if Math::Prime::Util::prime_get_config->{'gmp'};
   $version_str .= "\nWritten by Dana Jacobsen.\n";
@@ -69,6 +67,7 @@ sub eval_expr {
     $expr =~ s/:$mpu_func_map{$func}\(/Math::Prime::Util::$func(/g;
   }
   $expr =~ s/(\d+)/ Math::BigInt->new("$1") /g;
+  $expr = 'use Math::BigInt try=>"GMP"; ' . $expr;
   my $res = eval $expr; ## no critic
   die "Cannot eval: $expr\n" if !defined $res;
   $res = int($res->bstr) if ref($res) eq 'Math::BigInt' && $res <= ~0;
diff --git a/examples/test-factor-gnufactor.pl 
b/examples/test-factor-gnufactor.pl
index 229461e..f5614ba 100755
--- a/examples/test-factor-gnufactor.pl
+++ b/examples/test-factor-gnufactor.pl
@@ -150,7 +150,7 @@ sub mpu_factors {
     my @ns = @_;
     my $numpercommand = int( (4000-30)/(length($ns[-1])+1) );
     while (@ns) {
-      my $cs = join(" ", 'factor.pl', splice(@ns, 0, $numpercommand));
+      my $cs = join(" ", 'perl -Iblib/lib -Iblib/arch bin/factor.pl', 
splice(@ns, 0, $numpercommand));
       my $fout = qx{$cs};
       my @flines = split(/\n/, $fout);
       foreach my $fline (@flines) {
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d2f5ec6..bcd2488 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -114,6 +114,12 @@ BEGIN {
     *mertens       = \&Math::Prime::Util::_generic_mertens;
     *prime_count   = \&Math::Prime::Util::_generic_prime_count;
     *nth_prime     = \&Math::Prime::Util::PP::nth_prime;
+    *nth_prime_upper=\&Math::Prime::Util::PP::nth_prime_upper;
+    *nth_prime_lower=\&Math::Prime::Util::PP::nth_prime_lower;
+    *nth_prime_approx=\&Math::Prime::Util::PP::nth_prime_approx;
+    *prime_count_upper=\&Math::Prime::Util::PP::prime_count_upper;
+    *prime_count_lower=\&Math::Prime::Util::PP::prime_count_lower;
+    *prime_count_approx=\&Math::Prime::Util::PP::prime_count_approx;
     *carmichael_lambda = \&Math::Prime::Util::_generic_carmichael_lambda;
     *kronecker     = \&Math::Prime::Util::_generic_kronecker;
     *divisor_sum   = \&Math::Prime::Util::_generic_divisor_sum;
@@ -286,29 +292,6 @@ sub _upgrade_to_float {
   return Math::BigFloat->new($_[0]);
 }
 
-# TODO: Remove these, push more funcs into XS
-my @_primes_small = (
-   0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
-   101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
-   193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,
-   293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,
-   409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509);
-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;
-}
-
-
-
-
 
 #############################################################################
 
@@ -822,36 +805,8 @@ sub _generic_divisors {
   my($n) = @_;
   _validate_num($n) || _validate_positive_integer($n);
 
-  # In scalar context, returns sigma_0(n).  Very fast.
-  return divisor_sum($n,0) unless wantarray;
-  return ($n == 0) ? (0,1) : (1)  if $n <= 1;
-
-  my %all_factors;
-  my @factors = factor($n);
-  return (1,$n) if scalar @factors == 1;
-
-  if (ref($n) eq 'Math::BigInt') {
-    foreach my $f1 (@factors) {
-      my $big_f1 = Math::BigInt->new("$f1");
-      my @to_add = map { ($_ <= ''.~0) ? _bigint_to_int($_) : $_ }
-                   grep { $_ < $n }
-                   map { $big_f1 * $_ }
-                   keys %all_factors;
-      undef @all_factors{ $f1, @to_add };
-    }
-  } else {
-    foreach my $f1 (@factors) {
-      my @to_add = grep { $_ < $n }
-                   map { $f1 * $_ }
-                   keys %all_factors;
-      undef @all_factors{ $f1, @to_add };
-    }
-  }
-  # Add 1 and n
-  undef $all_factors{1};
-  undef $all_factors{$n};
-  my @divisors = sort {$a<=>$b} keys %all_factors;
-  return @divisors;
+  require Math::Prime::Util::PP;
+  return Math::Prime::Util::PP::divisors($n);
 }
 
 
@@ -1010,282 +965,6 @@ sub verify_prime {
 
 #############################################################################
 
-sub prime_count_approx {
-  my($x) = @_;
-  _validate_num($x) || _validate_positive_integer($x);
-
-  # With XS using 30k tables, this is super fast.
-  return prime_count($x) if $x < $_XS_MAXVAL && $x < 3_000_000;
-  # Give an exact answer for what we have in our little table.
-  return _tiny_prime_count($x) if $x < $_primes_small[-1];
-
-  # Below 2^58th or so, all differences between the high precision and C double
-  # precision result are less than 0.5.
-  if ($x <= $_XS_MAXVAL && $x <= 144115188075855872) {
-    return int(_XS_RiemannR($x) + 0.5);
-  }
-
-  # 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%          .06%
-  #    average bounds      .01%          .0002%
-  #    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.  If MPFR
-  # is available, then use it -- it rocks.  Otherwise, switch to LiCorr for
-  # very big values.  This is hacky and shouldn't be necessary.
-  my $result;
-  if ( $x < 1e36 || Math::Prime::Util::PP::_MPFR_available() ) {
-    if (ref($x) eq 'Math::BigFloat') {
-      # Make sure we get enough accuracy, and also not too much more than 
needed
-      $x->accuracy(length($x->bfloor->bstr())+2);
-    }
-    $result = RiemannR($x) + 0.5;
-  } else {
-    $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2);
-  }
-
-  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);
-
-  # With XS using 30k tables, this is super fast.
-  return prime_count($x) if $x < $_XS_MAXVAL && $x < 3_000_000;
-  # 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';
-
-  my $flogx = log($x);
-
-  # 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
-  # The Dusart (1999 or 2010) bounds are far, far better than the others.
-
-  my $result;
-  if ($x > 1000_000_000_000 && $_Config{'assume_rh'}) {
-    my $lix = LogarithmicIntegral($x);
-    my $sqx = sqrt($x);
-    # Schoenfeld bound:    (constant is 8 * Pi)
-    $result = $lix - (($sqx*$flogx) / 25.13274122871834590770114707);
-  } elsif ($x < 599) {
-    $result = $x / ($flogx - 0.7);   # For smaller numbers this works out well.
-  } else {
-    my $a;
-    # Hand tuned for small numbers (< 60_000M)
-    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 <     176000) { $a = 1.80; }
-    elsif ($x <     315000) { $a = 2.10; }
-    elsif ($x <    1100000) { $a = 2.20; }
-    elsif ($x <    4500000) { $a = 2.31; }
-    elsif ($x <  233000000) { $a = 2.36; }
-    elsif ($x < 5433800000) { $a = 2.32; }
-    elsif ($x <60000000000) { $a = 2.15; }
-    else                    { $a = 2.00; } # Dusart 2010, page 2
-    $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx));
-  }
-
-  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);
-
-  # With XS using 30k tables, this is super fast.
-  return prime_count($x) if $x < $_XS_MAXVAL && $x < 3_000_000;
-  # 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
-  # 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
-
-  # As with the lower bounds, Dusart bounds are best by far.
-
-  # Another possibility here for numbers under 3000M is to use Li(x)
-  # minus a correction.
-
-  my $flogx = log($x);
-
-  my $result;
-  if ($x > 10000_000_000_000 && $_Config{'assume_rh'}) {
-    my $lix = LogarithmicIntegral($x);
-    my $sqx = sqrt($x);
-    # Schoenfeld bound:    (constant is 8 * Pi)
-    $result = $lix + (($sqx*$flogx) / 25.13274122871834590770114707);
-  } elsif ($x <  1621) { $result = ($x / ($flogx - 1.048)) + 1.0; }
-    elsif ($x <  5000) { $result = ($x / ($flogx - 1.071)) + 1.0; }
-    elsif ($x < 15900) { $result = ($x / ($flogx - 1.098)) + 1.0; }
-    else {
-    my $a;
-    # Hand tuned for small numbers (< 60_000M)
-    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
-    #elsif ($x <60000000000) { $a = 2.362; }
-    #else                    { $a = 2.51;  } # Dusart 1999, page 14
-
-    # Old versions of Math::BigFloat will do the Wrong Thing with this.
-    $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0;
-  }
-
-  return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 
'Math::BigFloat';
-  return int($result);
-}
-
-#############################################################################
-
-sub nth_prime_approx {
-  my($n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
-
-  return $_primes_small[$n] if $n <= $#_primes_small;
-
-  $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 +=  7.52* $order; }
-  elsif ($n <       1271) { $approx +=  5.6 * $order; }
-  elsif ($n <       2000) { $approx +=  5.2 * $order; }
-  elsif ($n <       4000) { $approx +=  4.3 * $order; }
-  elsif ($n <      12000) { $approx +=  3.0 * $order; }
-  elsif ($n <     150000) { $approx +=  2.1 * $order; }
-  elsif ($n <  200000000) { $approx +=  0.0 * $order; }
-  else                    { $approx += -0.010 * $order; }
-  # $approx = -0.025 is better for the last, but it gives problems with some
-  # other code that always wants the asymptotic approximation to be >= actual.
-
-  return int($approx + 0.5);
-}
-
-# 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 $_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));
-
-  return int($lower);
-}
-
-# The nth prime will be less or equal to this number
-sub nth_prime_upper {
-  my($n) = @_;
-  _validate_num($n) || _validate_positive_integer($n);
-
-  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 >= 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);
-}
-
-
-#############################################################################
-
-
 #############################################################################
 
 sub RiemannZeta {
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index e72ae6d..c54d71e 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -76,6 +76,12 @@ sub _bigint_to_int {
                             : int("$_[0]");
 }
 
+sub _upgrade_to_float {
+  do { require Math::BigFloat; Math::BigFloat->import(); }
+    if !defined $Math::BigFloat::VERSION;
+  return Math::BigFloat->new($_[0]);
+}
+
 sub _validate_num {
   my($n, $min, $max) = @_;
   croak "Parameter must be defined" if !defined $n;
@@ -124,11 +130,7 @@ my @_primes_small = (
    101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,
    193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,
    293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,
-   409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499);
-my @_prime_count_small = (
-   0,0,1,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,8,8,8,8,9,9,9,9,9,9,10,10,
-   11,11,11,11,11,11,12,12,12,12,13,13,14,14,14,14,15,15,15,15,15,15,
-   16,16,16,16,16,16,17,17,18,18,18,18,18,18,19);
+   409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509);
 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,
@@ -139,6 +141,19 @@ 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);
 
+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) = @_;
 
@@ -966,6 +981,263 @@ sub nth_prime {
   $prime;
 }
 
+# The nth prime will be less or equal to this number
+sub nth_prime_upper {
+  my($n) = @_;
+  _validate_positive_integer($n);
+
+  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 >= 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 $_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));
+
+  return int($lower);
+}
+
+sub nth_prime_approx {
+  my($n) = @_;
+  _validate_num($n) || _validate_positive_integer($n);
+
+  return $_primes_small[$n] if $n <= $#_primes_small;
+
+  $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 +=  7.52* $order; }
+  elsif ($n <       1271) { $approx +=  5.6 * $order; }
+  elsif ($n <       2000) { $approx +=  5.2 * $order; }
+  elsif ($n <       4000) { $approx +=  4.3 * $order; }
+  elsif ($n <      12000) { $approx +=  3.0 * $order; }
+  elsif ($n <     150000) { $approx +=  2.1 * $order; }
+  elsif ($n <  200000000) { $approx +=  0.0 * $order; }
+  else                    { $approx += -0.010 * $order; }
+  # $approx = -0.025 is better for the last, but it gives problems with some
+  # other code that always wants the asymptotic approximation to be >= actual.
+
+  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%          .06%
+  #    average bounds      .01%          .0002%
+  #    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.  If MPFR
+  # is available, then use it -- it rocks.  Otherwise, switch to LiCorr for
+  # very big values.  This is hacky and shouldn't be necessary.
+  my $result;
+  if ( $x < 1e36 || _MPFR_available() ) {
+    if (ref($x) eq 'Math::BigFloat') {
+      # Make sure we get enough accuracy, and also not too much more than 
needed
+      $x->accuracy(length($x->bfloor->bstr())+2);
+    }
+    $result = RiemannR($x) + 0.5;
+  } else {
+    $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2);
+  }
+
+  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 $flogx = log($x);
+
+  # 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
+  # The Dusart (1999 or 2010) bounds are far, far better than the others.
+
+  my $result;
+  if ($x > 1000_000_000_000 && 
Math::Prime::Util::prime_get_config()->{'assume_rh'}) {
+    my $lix = LogarithmicIntegral($x);
+    my $sqx = sqrt($x);
+    # Schoenfeld bound:    (constant is 8 * Pi)
+    $result = $lix - (($sqx*$flogx) / 25.13274122871834590770114707);
+  } elsif ($x < 599) {
+    $result = $x / ($flogx - 0.7);   # For smaller numbers this works out well.
+  } else {
+    my $a;
+    # Hand tuned for small numbers (< 60_000M)
+    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 <     176000) { $a = 1.80; }
+    elsif ($x <     315000) { $a = 2.10; }
+    elsif ($x <    1100000) { $a = 2.20; }
+    elsif ($x <    4500000) { $a = 2.31; }
+    elsif ($x <  233000000) { $a = 2.36; }
+    elsif ($x < 5433800000) { $a = 2.32; }
+    elsif ($x <60000000000) { $a = 2.15; }
+    else                    { $a = 2.00; } # Dusart 2010, page 2
+    $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx));
+  }
+
+  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
+  # 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
+
+  # As with the lower bounds, Dusart bounds are best by far.
+
+  # Another possibility here for numbers under 3000M is to use Li(x)
+  # minus a correction.
+
+  my $flogx = log($x);
+
+  my $result;
+  if ($x > 10000_000_000_000 && 
Math::Prime::Util::prime_get_config()->{'assume_rh'}) {
+    my $lix = LogarithmicIntegral($x);
+    my $sqx = sqrt($x);
+    # Schoenfeld bound:    (constant is 8 * Pi)
+    $result = $lix + (($sqx*$flogx) / 25.13274122871834590770114707);
+  } elsif ($x <  1621) { $result = ($x / ($flogx - 1.048)) + 1.0; }
+    elsif ($x <  5000) { $result = ($x / ($flogx - 1.071)) + 1.0; }
+    elsif ($x < 15900) { $result = ($x / ($flogx - 1.098)) + 1.0; }
+    else {
+    my $a;
+    # Hand tuned for small numbers (< 60_000M)
+    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
+    #elsif ($x <60000000000) { $a = 2.362; }
+    #else                    { $a = 2.51;  } # Dusart 1999, page 14
+
+    # Old versions of Math::BigFloat will do the Wrong Thing with this.
+    $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0;
+  }
+
+  return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 
'Math::BigFloat';
+  return int($result);
+}
+
+
+#############################################################################
+
 sub _mulmod {
   my($x, $y, $n) = @_;
   return (($x * $y) % $n) if ($x|$y) < MPU_HALFWORD;
@@ -2553,6 +2825,41 @@ sub ecm_factor {
   @factors;
 }
 
+sub divisors {
+  my($n) = @_;
+  _validate_num($n) || _validate_positive_integer($n);
+
+  # 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;
+
+  my %all_factors;
+  my @factors = Math::Prime::Util::factor($n);
+  return (1,$n) if scalar @factors == 1;
+
+  if (ref($n) eq 'Math::BigInt') {
+    foreach my $f1 (@factors) {
+      my $big_f1 = Math::BigInt->new("$f1");
+      my @to_add = map { ($_ <= ''.~0) ? _bigint_to_int($_) : $_ }
+                   grep { $_ < $n }
+                   map { $big_f1 * $_ }
+                   keys %all_factors;
+      undef @all_factors{ $f1, @to_add };
+    }
+  } else {
+    foreach my $f1 (@factors) {
+      my @to_add = grep { $_ < $n }
+                   map { $f1 * $_ }
+                   keys %all_factors;
+      undef @all_factors{ $f1, @to_add };
+    }
+  }
+  # Add 1 and n
+  undef $all_factors{1};
+  undef $all_factors{$n};
+  my @divisors = sort {$a<=>$b} keys %all_factors;
+  return @divisors;
+}
 
 
 sub ExponentialIntegral {
diff --git a/lmo.c b/lmo.c
index 3d40af5..3db1ad1 100644
--- a/lmo.c
+++ b/lmo.c
@@ -293,8 +293,8 @@ static UV _phi_recurse(UV x, UV a) {
   UV i, c = (a > PHIC) ? PHIC : a;
   UV sum = tablephi(x, c);
   if (a > c) {
-    UV p  = _XS_nth_prime(c);
-    UV pa = _XS_nth_prime(a);
+    UV p  = nth_prime(c);
+    UV pa = nth_prime(a);
     for (i = c+1; i <= a; i++) {
       UV xp;
       p = next_prime(p);
@@ -344,7 +344,7 @@ UV legendre_phi(UV x, UV a) {
     uint32_t lastidx;
     UV res, max_cache_a = (a >= PHICACHEA) ? PHICACHEA : a+1;
     Newz(0, cache, PHICACHEX * max_cache_a, uint16_t);
-    primes = make_primelist(_XS_nth_prime(a+1), &lastidx);
+    primes = make_primelist(nth_prime(a+1), &lastidx);
     res = (UV) _phi(x, a, 1, primes, lastidx, cache);
     Safefree(primes);
     Safefree(cache);
diff --git a/util.c b/util.c
index 3a6dca8..e8e8e31 100644
--- a/util.c
+++ b/util.c
@@ -29,12 +29,14 @@
   extern long double logl(long double);
   extern long double fabsl(long double);
   extern long double floorl(long double);
+  extern long double ceill(long double);
 #else
   #define powl(x, y)  (long double) pow( (double) (x), (double) (y) )
   #define expl(x)     (long double) exp( (double) (x) )
   #define logl(x)     (long double) log( (double) (x) )
   #define fabsl(x)    (long double) fabs( (double) (x) )
   #define floorl(x)   (long double) floor( (double) (x) )
+  #define ceill(x)    (long double) ceil( (double) (x) )
 #endif
 
 #ifdef LDBL_INFINITY
@@ -570,7 +572,90 @@ UV _XS_prime_count(UV low, UV high)
   return count;
 }
 
+UV prime_count_approx(UV n)
+{
+  if (n < 3000000) return _XS_prime_count(2, n);
+  return (UV) (_XS_RiemannR( (long double) n ) + 0.5 );
+}
+
+UV prime_count_lower(UV n)
+{
+  long double fn, flogn, lower, a;
+
+  if (n < 33000) return _XS_prime_count(2, n);
+
+  fn     = (long double) n;
+  flogn  = logl(n);
+
+  if      (n <   176000)  a = 1.80;
+  else if (n <   315000)  a = 2.10;
+  else if (n <  1100000)  a = 2.20;
+  else if (n <  4500000)  a = 2.31;
+  else if (n <233000000)  a = 2.36;
+#if BITS_PER_WORD == 32
+  else a = 2.32;
+#else
+  else if (n < UVCONST( 5433800000)) a = 2.32;
+  else if (n < UVCONST(60000000000)) a = 2.15;
+  else a = 2.00;
+#endif
+
+  lower = fn/flogn * (1.0 + 1.0/flogn + a/(flogn*flogn));
+  return (UV) floorl(lower);
+}
+
+typedef struct {
+  UV thresh;
+  float aval;
+} thresh_t;
+
+static const thresh_t _upper_thresh[] = {
+  {     59000, 2.48 },
+  {    350000, 2.52 },
+  {    355991, 2.54 },
+  {    356000, 2.51 },
+  {   3550000, 2.50 },
+  {   3560000, 2.49 },
+  {   5000000, 2.48 },
+  {   8000000, 2.47 },
+  {  13000000, 2.46 },
+  {  18000000, 2.45 },
+  {  31000000, 2.44 },
+  {  41000000, 2.43 },
+  {  48000000, 2.42 },
+  { 119000000, 2.41 },
+  { 182000000, 2.40 },
+  { 192000000, 2.395 },
+  { 213000000, 2.390 },
+  { 271000000, 2.385 },
+  { 322000000, 2.380 },
+  { 400000000, 2.375 },
+  { 510000000, 2.370 },
+  { 682000000, 2.367 },
+  { UVCONST(2953652287), 2.362 }
+};
+#define NUPPER_THRESH (sizeof(_upper_thresh)/sizeof(_upper_thresh[0]))
+
+UV prime_count_upper(UV n)
+{
+  int i;
+  long double fn, flogn, upper, a;
+
+  if (n < 33000) return _XS_prime_count(2, n);
+
+  fn     = (long double) n;
+  flogn  = logl(n);
+
+  for (i = 0; i < NUPPER_THRESH; i++)
+    if (n < _upper_thresh[i].thresh)
+      break;
+
+  if (i < NUPPER_THRESH)   a = _upper_thresh[i].aval;
+  else                     a = 2.334;   /* Dusart 2010, page 2 */
 
+  upper = fn/flogn * (1.0 + 1.0/flogn + a/(flogn*flogn));
+  return (UV) ceill(upper);
+}
 
 static const unsigned short primes_small[] =
   {0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
@@ -580,18 +665,17 @@ static const unsigned short primes_small[] =
    409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499};
 #define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))
 
-/* Note: We're keeping this here because we use it for nth_prime */
 /* The nth prime will be less or equal to this number */
-static UV _XS_nth_prime_upper(UV n)
+UV nth_prime_upper(UV n)
 {
-  double fn, flogn, flog2n, upper;
+  long double fn, flogn, flog2n, upper;
 
   if (n < NPRIMES_SMALL)
     return primes_small[n];
 
-  fn     = (double) n;
-  flogn  = log(n);
-  flog2n = log(flogn);    /* Note distinction between log_2(n) and log^2(n) */
+  fn     = (long double) n;
+  flogn  = logl(n);
+  flog2n = logl(flogn);    /* Note distinction between log_2(n) and log^2(n) */
 
   if      (n >= 688383)    /* Dusart 2010 page 2 */
     upper = fn * (flogn + flog2n - 1.0 + ((flog2n-2.00)/flogn));
@@ -615,16 +699,73 @@ static UV _XS_nth_prime_upper(UV n)
    *    nth_prime_lower(n)  <=  nth_prime(n)  <=  nth_prime_upper(n)
    */
   /* Watch out for  overflow */
-  if (upper >= (double)UV_MAX) {
+  if (upper >= (long double)UV_MAX) {
     if (n <= MPU_MAX_PRIME_IDX) return MPU_MAX_PRIME;
     croak("nth_prime_upper(%"UVuf") overflow", n);
   }
 
-  return (UV) ceil(upper);
+  return (UV) ceill(upper);
+}
+
+/* The nth prime will be greater than or equal to this number */
+UV nth_prime_lower(UV n)
+{
+  long double fn, flogn, flog2n, lower;
+
+  if (n < NPRIMES_SMALL)
+    return primes_small[n];
+
+  fn     = (long double) n;
+  flogn  = logl(n);
+  flog2n = logl(flogn);    /* Note distinction between log_2(n) and log^2(n) */
+
+  /* Dusart 2010 page 2, for all n >= 3 */
+  lower = fn * (flogn + flog2n - 1.0 + ((flog2n-2.10)/flogn));
+
+  return (UV) floorl(lower);
+}
+
+UV nth_prime_approx(UV n)
+{
+  long double fn, flogn, flog2n, approx, order;
+
+  if (n < NPRIMES_SMALL)
+    return primes_small[n];
+
+  fn     = (long double) n;
+  flogn  = logl(n);
+  flog2n = logl(flogn);    /* Note distinction between log_2(n) and log^2(n) */
+
+  /* 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)
+   */
+
+  approx = fn * (  flogn + flog2n - 1.0
+                 + ((flog2n - 2.0) / flogn)
+                 - (((flog2n*flog2n) - 6.0*flog2n + 11.0) / (2*flogn*flogn))
+                );
+
+  /* Apply a correction */
+  order = flog2n / flogn;
+  order = order * order * order * fn;
+  if      (n <      259) { approx += 10.4  * order; }
+  else if (n <      775) { approx +=  7.52 * order; }
+  else if (n <     1271) { approx +=  5.6  * order; }
+  else if (n <     2000) { approx +=  5.2  * order; }
+  else if (n <     4000) { approx +=  4.3  * order; }
+  else if (n <    12000) { approx +=  3.0  * order; }
+  else if (n <   150000) { approx +=  2.1  * order; }
+  else if (n <200000000) {                          }
+  else                   { approx += -0.01 * order; } /* -0.25 is closer */
+
+  return (UV) floorl(approx + 0.5);
 }
 
 
-UV _XS_nth_prime(UV n)
+UV nth_prime(UV n)
 {
   const unsigned char* cache_sieve;
   unsigned char* segment;
@@ -638,7 +779,7 @@ UV _XS_nth_prime(UV n)
     return primes_small[n];
 
   /* Determine a bound on the nth prime.  We know it comes before this. */
-  upper_limit = _XS_nth_prime_upper(n);
+  upper_limit = nth_prime_upper(n);
   MPUassert(upper_limit > 0, "nth_prime got an upper limit of 0");
 
   /* For relatively small values, generate a sieve and count the results.
diff --git a/util.h b/util.h
index 715e4b7..76b9fb5 100644
--- a/util.h
+++ b/util.h
@@ -13,7 +13,13 @@ extern UV  next_prime(UV x);
 extern UV  prev_prime(UV x);
 
 extern UV  _XS_prime_count(UV low, UV high);
-extern UV  _XS_nth_prime(UV x);
+extern UV  nth_prime(UV x);
+extern UV  nth_prime_upper(UV x);
+extern UV  nth_prime_lower(UV x);
+extern UV  nth_prime_approx(UV x);
+extern UV  prime_count_lower(UV x);
+extern UV  prime_count_upper(UV x);
+extern UV  prime_count_approx(UV x);
 
 extern signed char* _moebius_range(UV low, UV high);
 extern UV*    _totient_range(UV low, UV high);

-- 
Alioth's /usr/local/bin/git-commit-notice on 
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to