[libmath-prime-util-perl] 15/59: Lots of bignum support
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:53 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.10
in repository libmath-prime-util-perl.
commit d94fc1102872ebb07e1e321aafc8c296ad2d5e25
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Jun 29 21:33:30 2012 -0600
Lots of bignum support
---
Makefile.PL | 3 +
TODO | 10 ++-
XS.xs | 8 +-
lib/Math/Prime/Util.pm | 225 +++++++++++++++++++++++++++++++++++++---------
lib/Math/Prime/Util/PP.pm | 93 +++++++++++++------
t/15-probprime.t | 2 +-
t/19-moebius.t | 55 ++++++++++++
util.c | 22 ++---
util.h | 6 +-
9 files changed, 335 insertions(+), 89 deletions(-)
diff --git a/Makefile.PL b/Makefile.PL
index a9afe32..b2f5988 100644
--- a/Makefile.PL
+++ b/Makefile.PL
@@ -23,6 +23,9 @@ WriteMakefile1(
'Carp' => 0,
'Tie::Array' => 0,
'base' => 0,
+ 'bignum' => 0,
+ 'Math::BigInt' => 0,
+ 'Math::BigFloat' => 0,
},
MIN_PERL_VERSION => 5.006002,
diff --git a/TODO b/TODO
index 828b213..f04f361 100644
--- a/TODO
+++ b/TODO
@@ -14,7 +14,7 @@
- Is the current PP.pm setup the way we want to do it?
- input validation (in XS, or do we need to make Perl wrappers for everything?)
- We can do inpuut val in XS by looking at the NV. But I think long term we'll
+ We can do input val in XS by looking at the NV. But I think long term we'll
have a little Perl front end for everything to route to bignum routines or
regular routines.
@@ -26,3 +26,11 @@
version does it in a painful way. Something simpler to be had?
- In PP nth_prime* and prime_count*, turn on bignum if bigint is on.
+
+- Add euler_phi documentation
+
+- More testing for bignum, including a test suite file
+
+- Move the approx / bounds functions into main, to support bignum properly.
+
+- need next_prime and prev_prime bignum support. random_ndigit_prime needs it.
diff --git a/XS.xs b/XS.xs
index d7f2f32..842e2b3 100644
--- a/XS.xs
+++ b/XS.xs
@@ -212,7 +212,7 @@ erat_primes(IN UV low, IN UV high)
void
-factor(IN UV n)
+XS_factor(IN UV n)
PPCODE:
if (n < 4) {
XPUSHs(sv_2mortal(newSVuv( n ))); /* If n is 0-3, we're done. */
@@ -371,10 +371,10 @@ int
is_prob_prime(IN UV n)
double
-ExponentialIntegral(double x)
+XS_ExponentialIntegral(double x)
double
-LogarithmicIntegral(double x)
+XS_LogarithmicIntegral(double x)
double
-RiemannR(double x)
+XS_RiemannR(double x)
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d81d3b5..f7511d6 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -12,6 +12,7 @@ BEGIN {
# use parent qw( Exporter );
use base qw( Exporter );
our @EXPORT_OK = qw(
+ prime_get_config
prime_precalc prime_memfree
is_prime is_prob_prime miller_rabin
primes
@@ -19,23 +20,29 @@ our @EXPORT_OK = qw(
prime_count prime_count_lower prime_count_upper prime_count_approx
nth_prime nth_prime_lower nth_prime_upper nth_prime_approx
random_prime random_ndigit_prime
- factor all_factors
+ factor all_factors moebius euler_phi
ExponentialIntegral LogarithmicIntegral RiemannR
);
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
-my $_pure_perl;
+my %_Config;
BEGIN {
+
+ # Load PP code. Nothing exported.
+ require Math::Prime::Util::PP;
+
eval {
require XSLoader;
XSLoader::load(__PACKAGE__, $Math::Prime::Util::VERSION);
prime_precalc(0);
- $_pure_perl = 0;
+ $_Config{'xs'} = 1;
+ $_Config{'maxbits'} = _maxbits();
1;
} or do {
- require Math::Prime::Util::PP;
- $_pure_perl = 1;
+die $@;
+ $_Config{'xs'} = 0;
+ $_Config{'maxbits'} = Math::Prime::Util::PP::_maxbits();
#carp "Using Pure Perl implementation";
*_get_prime_cache_size = \&Math::Prime::Util::PP::_get_prime_cache_size;
*_maxbits = \&Math::Prime::Util::PP::_maxbits;
@@ -72,18 +79,39 @@ BEGIN {
*RiemannR = \&Math::Prime::Util::PP::RiemannR;
*LogarithmicIntegral = \&Math::Prime::Util::PP::LogarithmicIntegral;
*ExponentialIntegral = \&Math::Prime::Util::PP::ExponentialIntegral;
- # TODO: We should have some tests to verify XS vs. PP.
}
}
END {
_prime_memfreeall;
}
-my $_maxparam = (_maxbits == 32) ? 4294967295 : 18446744073709551615;
-my $_maxdigits = (_maxbits == 32) ? 10 : 20;
-my $_maxprime = (_maxbits == 32) ? 4294967291 : 18446744073709551557;
-my $_maxprimeidx=(_maxbits == 32) ? 203280221 : 425656284035217743;
+if ($_Config{'maxbits'} == 32) {
+ $_Config{'maxparam'} = 4294967295;
+ $_Config{'maxdigits'} = 10;
+ $_Config{'maxprime'} = 4294967291;
+ $_Config{'maxprimeidx'} = 203280221;
+} else {
+ $_Config{'maxparam'} = 18446744073709551615;
+ $_Config{'maxdigits'} = 20;
+ $_Config{'maxprime'} = 18446744073709551557;
+ $_Config{'maxprimeidx'} = 425656284035217743;
+}
+sub prime_get_config {
+ my %config = %_Config;
+ return \%config;
+}
+
+sub _validate_positive_integer {
+ my($n, $min, $max) = @_;
+ croak "Parameter must be defined" if !defined $n;
+ croak "Parameter '$n' must be a positive integer" if $n =~ tr/0123456789//c;
+ croak "Parameter '$n' must be >= $min" if defined $min && $n < $min;
+ croak "Parameter '$n' must be <= $max" if defined $max && $n > $max;
+ croak "Parameter '$n' outside of integer range" if $n > $_Config{'maxparam'}
+ && ref($n) !~ /^Math::Big/;
+ 1;
+}
sub primes {
my $optref = (ref $_[0] eq 'HASH') ? shift : {};
@@ -91,25 +119,16 @@ sub primes {
croak "too many parameters to primes" unless scalar @_ <= 2;
my $low = (@_ == 2) ? shift : 2;
my $high = shift;
- my $sref = [];
- # Validate parameters
- if ( (!defined $low) || (!defined $high) ||
- ($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
- ) {
- $low = '<undef>' unless defined $low;
- $high = '<undef>' unless defined $high;
- croak "Parameters [ $low $high ] must be positive integers";
- }
-
- # Verify the parameters are in range.
- if (!$_pure_perl || !defined $bigint::VERSION) {
- croak "Parameters [ $low $high ] not in range 0-$_maxparam" unless $low <= $_maxparam && $high <= $_maxparam;
- }
+ _validate_positive_integer($low);
+ _validate_positive_integer($high);
+ my $sref = [];
return $sref if ($low > $high) || ($high < 2);
- return Math::Prime::Util::PP::primes($low,$high) if $_pure_perl;
+ if ( (!$_Config{'xs'}) || ($high > $_Config{'maxparam'}) ) {
+ return Math::Prime::Util::PP::primes($low,$high);
+ }
my $method = $optref->{'method'};
$method = 'Dynamic' unless defined $method;
@@ -162,19 +181,11 @@ sub primes {
sub random_prime {
my $low = (@_ == 2) ? shift : 2;
my $high = shift;
- if ( (!defined $low) || (!defined $high) ||
- ($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
- ) {
- $low = 'undef' unless defined $low;
- $high = 'undef' unless defined $high;
- croak "Parameters [ $low $high ] must be positive integers";
- }
- if (!$_pure_perl || !defined $bigint::VERSION) {
- croak "Parameters [ $low $high ] not in range 0-$_maxparam" unless $low <= $_maxparam && $high <= $_maxparam;
- }
- $low = 2 if $low < 2;
+ _validate_positive_integer($low);
+ _validate_positive_integer($high);
- # Make sure we have a valid range.
+ # Tighten the range to the nearest prime.
+ $low = 2 if $low < 2;
$low = next_prime($low - 1);
$high = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high);
return $low if ($low == $high) && is_prime($low);
@@ -220,15 +231,12 @@ sub random_prime {
my @_random_ndigit_ranges;
sub random_ndigit_prime {
- my $digits = shift;
- # TODO: bigint with many digits
- if ((!defined $digits) || ($digits > $_maxdigits) || ($digits < 1)) {
- croak "Digits must be between 1 and $_maxdigits";
- }
+ my($digits) = @_;
+ _validate_positive_integer($digits, 1, (defined $bigint::VERSION) ? 10000 : $_Config{'maxdigits'} );
if (!defined $_random_ndigit_ranges[$digits]) {
my $low = ($digits == 1) ? 0 : int(10 ** ($digits-1));
my $high = int(10 ** $digits);
- $high = ~0 if $high > ~0;
+ $high = ~0 if $high > ~0 && !defined $bigint::VERSION;
$_random_ndigit_ranges[$digits] = [next_prime($low), prev_prime($high)];
}
my ($low, $high) = @{$_random_ndigit_ranges[$digits]};
@@ -237,7 +245,6 @@ sub random_ndigit_prime {
# Perhaps a random_nbit_prime ? Definition?
-
sub all_factors {
my $n = shift;
my @factors = factor($n);
@@ -255,6 +262,126 @@ sub all_factors {
return @factors;
}
+# A008683 Moebius function mu(n)
+# A030059, A013929, A030229, A002321, A005117, A013929 all relate.
+
+# One can argue for the Omega function (A001221), Euler Phi (A000010), and
+# Merten's functions also.
+
+sub moebius {
+ my($n) = @_;
+ _validate_positive_integer($n, 1);
+ return if $n <= 0;
+ return 1 if $n == 1;
+
+ # Quick check for small replicated factors
+ return 0 if ($n >= 25) && (($n % 4) == 0 || ($n % 9) == 0 || ($n % 25) == 0);
+
+ my @factors = factor($n);
+ my %all_factors;
+ foreach my $factor (@factors) {
+ return 0 if $all_factors{$factor}++;
+ }
+ return (((scalar @factors) % 2) == 0) ? 1 : -1;
+}
+
+sub euler_phi {
+ my($n) = @_;
+ _validate_positive_integer($n);
+ return 1 if $n <= 1; # phi(0) is disputed
+
+ my %factor_mult;
+ my @factors = grep { !$factor_mult{$_}++ } factor($n);
+ my $totient = $n;
+ foreach my $factor (@factors) {
+ # These divisions should always be exact
+ $totient = int($totient/$factor) * ($factor-1);
+ }
+ $totient;
+}
+
+
+#############################################################################
+# Front ends to functions.
+#
+# These will do input validation, then call the appropriate internal function
+# based on the input (XS, GMP, PP).
+#############################################################################
+
+#
+# This works, but sadly has a lot of overhead -- 7x more overhead than the
+# entire is_prime C function for n under 10M or so.
+#
+#sub is_prime {
+# my($n) = @_;
+# _validate_positive_integer($n);
+#
+# return XS_is_prime($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
+# return Math::Prime::Util::PP::is_prime($n);
+#}
+
+sub factor {
+ my($n) = @_;
+ _validate_positive_integer($n);
+
+ return XS_factor($n) if $_Config{'xs'} && $n <= $_Config{'maxparam'};
+
+ $n = $n->as_number() if ref($n) =~ /^Math::BigFloat/;
+ $n = $n->numify if $n <= ~0;
+ Math::Prime::Util::PP::factor($n);
+}
+
+
+sub RiemannR {
+ my($n) = @_;
+ croak("Invalid input to ReimannR: x must be > 0") if $n <= 0;
+
+ if (ref($n) eq 'Math::BigInt' && !defined $Math::BigFloat::VERSION) {
+ require Math::BigFloat;
+ $n = new Math::BigFloat "$n";
+ }
+
+ return Math::Prime::Util::PP::RiemannR($n, 1e-30) if ref($n) =~ /^Math::Big/;
+ return Math::Prime::Util::PP::RiemannR($n) if !$_Config{'xs'};
+ return XS_RiemannR($n);
+
+ # We could make a new object, like:
+ # require Math::BigFloat;
+ # my $bign = new Math::BigFloat "$n";
+ # my $result = Math::Prime::Util::PP::RiemannR($bign);
+ # return $result;
+}
+
+sub ExponentialIntegral {
+ my($n) = @_;
+ croak "Invalid input to ExponentialIntegral: x must be != 0" if $n == 0;
+
+ if (ref($n) eq 'Math::BigInt' && !defined $Math::BigFloat::VERSION) {
+ require Math::BigFloat;
+ $n = new Math::BigFloat "$n";
+ }
+
+ return Math::Prime::Util::PP::ExponentialIntegral($n, 1e-30) if ref($n) =~ /^Math::Big/;
+ return Math::Prime::Util::PP::ExponentialIntegral($n) if !$_Config{'xs'};
+ return XS_ExponentialIntegral($n);
+}
+
+sub LogarithmicIntegral {
+ my($n) = @_;
+ return 0 if $n == 0;
+ croak("Invalid input to LogarithmicIntegral: x must be >= 0") if $n <= 0;
+
+ if (defined $Math::BigFloat::VERSION) {
+ return Math::BigFloat->binf('-') if $n == 1;
+ return Math::BigFloat->new('1.045163780117492784844588889194613136522615578151201575832909144075013205210359530172717405626383356306') if $n == 2;
+ } else {
+ return 0+'-inf' if $n == 1;
+ return 1.045163780117492784844588889194613136522615578151 if $n == 2;
+ }
+ ExponentialIntegral(log($n));
+}
+
+
use Math::Prime::Util::MemFree;
1;
@@ -619,6 +746,16 @@ digits between 1 and the maximum native type (10 for 32-bit, 20 for 64-bit).
One of the primes within that range (e.g. 1000 - 9999 for 4-digits) will be
uniformly selected using the L<rand> function.
+=head2 moebius
+
+ say "$n is square free" if moebius($n) != 0;
+ $sum += moebius($_) for (1..200); say "Mertens(200) = $sum";
+
+Returns the Möbius function (also called the Moebius, Mobius, or MoebiusMu
+function) for a positive non-zero integer input. This function is 1 if
+C<n = 1>, 0 if C<n> is not square free (i.e. C<n> has a repeated factor),
+and C<-1^t> if C<n> is a product of C<t> distinct primes. This is an
+important function in prime number theory.
=head1 UTILITY FUNCTIONS
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 84613d6..74047d9 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -584,6 +584,18 @@ sub _mulmod {
$r;
}
+sub _native_powmod {
+ my($n, $power, $m) = @_;
+ my $t = 1;
+ $n = $n % $m;
+ while ($power) {
+ $t = ($t * $n) % $m if ($power & 1) != 0;
+ $n = ($n * $n) % $m;
+ $power >>= 1;
+ }
+ $t;
+}
+
sub _powmod {
my($n, $power, $m) = @_;
my $t = 1;
@@ -639,20 +651,41 @@ sub miller_rabin {
$d >>= 1;
}
- foreach my $a (@bases) {
- croak "Base $a is invalid" if $a < 2;
- my $x = _powmod($a, $d, $n);
- next if ($x == 1) || ($x == ($n-1));
+ if ( ($n < $_half_word) || (ref($n) =~ /^Math::Big/) ) {
+
+ foreach my $a (@bases) {
+ croak "Base $a is invalid" if $a < 2;
+ my $x = _native_powmod($a, $d, $n);
+ next if ($x == 1) || ($x == ($n-1));
+ foreach my $r (1 .. $s) {
+ $x = ($x*$x) % $n;
+ return 0 if $x == 1;
+ if ($x == ($n-1)) {
+ $a = 0;
+ last;
+ }
+ }
+ return 0 if $a != 0;
+ }
- foreach my $r (1 .. $s) {
- $x = ($x < $_half_word) ? ($x*$x) % $n : _mulmod($x, $x, $n);
- return 0 if $x == 1;
- if ($x == ($n-1)) {
- $a = 0;
- last;
+ } else {
+
+ foreach my $a (@bases) {
+ croak "Base $a is invalid" if $a < 2;
+ my $x = _powmod($a, $d, $n);
+ next if ($x == 1) || ($x == ($n-1));
+
+ foreach my $r (1 .. $s) {
+ $x = ($x < $_half_word) ? ($x*$x) % $n : _mulmod($x, $x, $n);
+ return 0 if $x == 1;
+ if ($x == ($n-1)) {
+ $a = 0;
+ last;
+ }
}
+ return 0 if $a != 0;
}
- return 0 if $a != 0;
+
}
1;
}
@@ -682,10 +715,15 @@ sub is_prob_prime {
sub _basic_factor {
# MODIFIES INPUT SCALAR
return ($_[0]) if $_[0] < 4;
+
my @factors;
- while ( ($_[0] % 2) == 0 ) { push @factors, 2; $_[0] /= 2; }
- while ( ($_[0] % 3) == 0 ) { push @factors, 3; $_[0] /= 3; }
- while ( ($_[0] % 5) == 0 ) { push @factors, 5; $_[0] /= 5; }
+ while ( ($_[0] % 2) == 0 ) { push @factors, 2; $_[0] = int($_[0] / 2); }
+ while ( ($_[0] % 3) == 0 ) { push @factors, 3; $_[0] = int($_[0] / 3); }
+ while ( ($_[0] % 5) == 0 ) { push @factors, 5; $_[0] = int($_[0] / 5); }
+
+ # Stop using bignum if we can
+ $_[0] = $_[0]->numify if ref($_[0]) =~ /^Math::Big/ && $_[0] <= ~0;
+
if ( ($_[0] > 1) && _is_prime7($_[0]) ) {
push @factors, $_[0];
$_[0] = 1;
@@ -706,7 +744,7 @@ sub trial_factor {
if ( ($n % $f) == 0) {
while ( ($n % $f) == 0) {
push @factors, $f;
- $n /= $f;
+ $n = int($n/$f);
}
$limit = int( sqrt($n) + 0.001);
}
@@ -725,13 +763,16 @@ sub factor {
my @factors = _basic_factor($n);
return @factors if $n < 4;
- while (($n % 7) == 0) { push @factors, 7; $n /= 7; }
- while (($n % 11) == 0) { push @factors, 11; $n /= 11; }
- while (($n % 13) == 0) { push @factors, 13; $n /= 13; }
- while (($n % 17) == 0) { push @factors, 17; $n /= 17; }
- while (($n % 19) == 0) { push @factors, 19; $n /= 19; }
- while (($n % 23) == 0) { push @factors, 23; $n /= 23; }
- while (($n % 29) == 0) { push @factors, 29; $n /= 29; }
+ while (($n % 7) == 0) { push @factors, 7; $n = int($n / 7); }
+ while (($n % 11) == 0) { push @factors, 11; $n = int($n / 11); }
+ while (($n % 13) == 0) { push @factors, 13; $n = int($n / 13); }
+ while (($n % 17) == 0) { push @factors, 17; $n = int($n / 17); }
+ while (($n % 19) == 0) { push @factors, 19; $n = int($n / 19); }
+ while (($n % 23) == 0) { push @factors, 23; $n = int($n / 23); }
+ while (($n % 29) == 0) { push @factors, 29; $n = int($n / 29); }
+
+ # Stop using bignum if possible
+ $n = $n->numify if ref($n) =~ /^Math::Big/ && $n <= ~0;
my @nstack = ($n);
while (@nstack) {
@@ -885,8 +926,8 @@ my $_const_euler = 0.57721566490153286060651209008240243104215933593992;
my $_const_li2 = 1.045163780117492784844588889194613136522615578151;
sub ExponentialIntegral {
- my($x) = @_;
- my $tol = 1e-16;
+ my($x, $tol) = @_;
+ $tol = 1e-16 unless defined $tol;
my $sum = 0.0;
my($y, $t);
my $c = 0.0;
@@ -1032,8 +1073,8 @@ sub _evaluate_zeta {
# Riemann R function
sub RiemannR {
- my($x) = @_;
- my $tol = 1e-16;
+ my($x, $tol) = @_;
+ $tol = 1e-16 unless defined $tol;
my $sum = 0.0;
my($y, $t);
my $c = 0.0;
diff --git a/t/15-probprime.t b/t/15-probprime.t
index 79a8779..67c14ee 100644
--- a/t/15-probprime.t
+++ b/t/15-probprime.t
@@ -5,7 +5,7 @@ use warnings;
use Test::More;
use Math::Prime::Util qw/is_prime is_prob_prime miller_rabin/;
-my $use64 = Math::Prime::Util::_maxbits > 32;
+my $use64 = Math::Prime::Util::prime_get_config->{'maxbits'} > 32;
my $extra = defined $ENV{RELEASE_TESTING} && $ENV{RELEASE_TESTING};
plan tests => 6 + 19 + 3573 + (5 + 29 + 22 + 23 + 16) + 15 + 27
diff --git a/t/19-moebius.t b/t/19-moebius.t
new file mode 100644
index 0000000..9ab422c
--- /dev/null
+++ b/t/19-moebius.t
@@ -0,0 +1,55 @@
+#!/usr/bin/env perl
+use strict;
+use warnings;
+
+use Test::More;
+use Math::Prime::Util qw/moebius euler_phi/;
+
+my @moeb_vals = (qw/ 1 -1 -1 0 -1 1 -1 0 0 1 -1 0 -1 1 1 0 -1 0 -1 0 /);
+my %mertens = (
+ 1 => 1,
+ 2 => 0,
+ 3 => -1,
+ 4 => -1,
+ 5 => -2,
+ 10 => -1,
+ 100 => 1,
+ 1000 => 2,
+ 10000 => -23,
+ 100000 => -48,
+# 1000000 => 212,
+# 10000000 => 1037,
+);
+
+my %totients = (
+ 1 => 1,
+ 2 => 1,
+ 20 => 8,
+ 36 => 12,
+ 123456 => 41088,
+ 123457 => 123456,
+ 123456789 => 82260072,
+);
+
+plan tests => 0 + 1
+ + scalar @moeb_vals
+ + scalar(keys %mertens)
+ + scalar(keys %totients);
+
+ok(!eval { moebius(0); }, "moebius(0)");
+
+my $i = 1;
+foreach my $m (@moeb_vals) {
+ is( moebius($i), $m, "moebius($i) == $m" );
+ $i++;
+}
+
+while (my($n, $mertens) = each (%mertens)) {
+ my $M = 0;
+ $M += moebius($_) for (1 .. $n);
+ is( $M, $mertens, "Mertens($n) == $mertens" );
+}
+
+while (my($n, $phi) = each (%totients)) {
+ is( euler_phi($n), $phi, "euler_phi($n) == $phi" );
+}
diff --git a/util.c b/util.c
index 749ddba..174d68f 100644
--- a/util.c
+++ b/util.c
@@ -508,11 +508,12 @@ UV prime_count_approx(UV x)
* Riemann's R function works even better, with errors of ~1828 at 10^10
* and 24M at 10^19.
*
- * Method 10^10 %error 10^19 %error
- * --------------- ------------ ------------
- * average bounds .01% .0002%
- * li(n) .0007% .00000004%
- * R(n) .0004% .00000001%
+ * Method 10^10 %error 10^19 %error
+ * ----------------- ------------ ------------
+ * average bounds .01% .0002%
+ * li(n) .0007% .00000004%
+ * li(n)-li(n^.5)/2 .0004% .00000001%
+ * R(n) .0004% .00000001%
*
* Getting fancier, one try using Riemann's pi formula:
* http://trac.sagemath.org/sage_trac/ticket/8135
@@ -521,7 +522,8 @@ UV prime_count_approx(UV x)
if (x < NPRIME_COUNT_SMALL)
return prime_count_small[x];
- R = RiemannR(x);
+ //R = XS_LogarithmicIntegral(x) - XS_LogarithmicIntegral(sqrt(x))/2;
+ R = XS_RiemannR(x);
/* We could add the additional factor:
* R = R - (1.0 / log(x)) + (M_1_PI * atan(M_PI/log(x)))
@@ -833,7 +835,7 @@ UV nth_prime(UV n)
static double const euler_mascheroni = 0.57721566490153286060651209008240243104215933593992;
static double const li2 = 1.045163780117492784844588889194613136522615578151;
-double ExponentialIntegral(double x) {
+double XS_ExponentialIntegral(double x) {
double const tol = 1e-16;
double val, term, fact_n;
double y, t;
@@ -925,12 +927,12 @@ double ExponentialIntegral(double x) {
return val;
}
-double LogarithmicIntegral(double x) {
+double XS_LogarithmicIntegral(double x) {
if (x == 0) return 0;
if (x == 1) return -INFINITY;
if (x == 2) return li2;
if (x <= 0) croak("Invalid input to LogarithmicIntegral: x must be > 0");
- return ExponentialIntegral(log(x));
+ return XS_ExponentialIntegral(log(x));
}
/*
@@ -1006,7 +1008,7 @@ static double evaluate_zeta(double x) {
return sum;
}
-double RiemannR(double x) {
+double XS_RiemannR(double x) {
double const tol = 1e-16;
double y, t, part_term, term, flogx, zeta;
double sum = 0.0;
diff --git a/util.h b/util.h
index a506ec2..b93c075 100644
--- a/util.h
+++ b/util.h
@@ -19,9 +19,9 @@ extern UV nth_prime_upper(UV x);
extern UV nth_prime_approx(UV x);
extern UV nth_prime(UV x);
-extern double ExponentialIntegral(double x);
-extern double LogarithmicIntegral(double x);
-extern double RiemannR(double x);
+extern double XS_ExponentialIntegral(double x);
+extern double XS_LogarithmicIntegral(double x);
+extern double XS_RiemannR(double x);
/* Above this value, is_prime will do deterministic Miller-Rabin */
/* With 64-bit math, we can do much faster mulmods from 2^16-2^32 */
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list