This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.15
in repository libmath-prime-util-perl.
commit e04d2d3a57f6d993a8b6377a145d25be89708833
Author: Dana Jacobsen <d...@acm.org>
Date: Sun Dec 2 04:48:41 2012 -0800
Use Math::MPFR if possible for Ei/li/Zeta/R functions. Huge speedup and
accuracy gains for BigFloats.
---
Changes | 13 ++
TODO | 5 -
lib/Math/Prime/Util.pm | 110 +++++++++++-----
lib/Math/Prime/Util/PP.pm | 243 ++++++++++++++++++++++++++++++++----
lib/Math/Prime/Util/ZetaBigFloat.pm | 86 +++++++++----
util.c | 25 ++--
6 files changed, 387 insertions(+), 95 deletions(-)
diff --git a/Changes b/Changes
index f061020..c660517 100644
--- a/Changes
+++ b/Changes
@@ -1,5 +1,18 @@
Revision history for Perl extension Math::Prime::Util.
+0.15 ?? December 2012
+
+ - Lots of internal changes to Ei, li, Zeta, and R functions:
+ - Native Zeta and R have slightly more accurate results.
+ - For bignums, use Math::MPFR if possible. MUCH faster.
+ Also allows extended precision while still being fast.
+ - Better accuracy for standard bignums.
+ - All four functions do:
+ - XS if native input.
+ - MPFR to whatever accuracy is desired, if Math::MPFR installed.
+ - BigFloat versions if no MPFR and BigFloat input.
+ - standard version if no MPFR and not a BigFloat.
+
0.14 29 November 2012
- Compilation and test issues:
diff --git a/TODO b/TODO
index f01fe4b..a83b661 100644
--- a/TODO
+++ b/TODO
@@ -29,11 +29,6 @@
- Make proper pminus1 in PP code, like factor.c.
-- For bignums, RiemannZeta and RiemmannR are slow and give questionable
- precision. We should be able to do better. One problem is the accuracy
- bug in Math::BigFloat. Perhaps check for Math::MPFR installed and use it?
- Alternately write real-number pow function using GMP.
-
- An assembler version of mulmod for i386 would be _really_ helpful for
all the non-x86-64 Intel machines.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 0293972..855eec5 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1172,9 +1172,11 @@ sub prime_count_approx {
# my $result = int(LogarithmicIntegral($x) -
LogarithmicIntegral(sqrt($x))/2);
- my $xlen = (ref($x) eq 'Math::BigFloat') ? length($x->bfloor->bstr()) :
length(int($x));
- my $tol = 10**-$xlen;
- my $result = RiemannR($x, $tol) + 0.5;
+ 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);
+ }
+ my $result = RiemannR($x) + 0.5;
return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq
'Math::BigFloat';
return int($result);
@@ -1418,18 +1420,18 @@ sub RiemannZeta {
my($n) = @_;
croak("Invalid input to ReimannZeta: x must be > 0") if $n <= 0;
- return Math::Prime::Util::PP::RiemannZeta($n) if defined $bignum::VERSION ||
ref($n) eq 'Math::BigFloat';
- return _XS_RiemannZeta($n) if $n <= $_XS_MAXVAL;
+ return _XS_RiemannZeta($n)
+ if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $n <=
$_XS_MAXVAL;
return Math::Prime::Util::PP::RiemannZeta($n);
}
sub RiemannR {
- my($n, $tol) = @_;
+ my($n) = @_;
croak("Invalid input to ReimannR: x must be > 0") if $n <= 0;
- return Math::Prime::Util::PP::RiemannR($n, $tol) if defined $bignum::VERSION
|| ref($n) eq 'Math::BigFloat';
- return _XS_RiemannR($n) if $n <= $_XS_MAXVAL;
- return Math::Prime::Util::PP::RiemannR($n, $tol);
+ return _XS_RiemannR($n)
+ if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' && $n <=
$_XS_MAXVAL;
+ return Math::Prime::Util::PP::RiemannR($n);
}
sub ExponentialIntegral {
@@ -1438,9 +1440,10 @@ sub ExponentialIntegral {
return 0 if $n == $_Neg_Infinity;
return $_Infinity if $n == $_Infinity;
- return Math::Prime::Util::PP::ExponentialIntegral($n) if defined
$bignum::VERSION || ref($n) eq 'Math::BigFloat';
- return Math::Prime::Util::PP::ExponentialIntegral($n) if !$_Config{'xs'};
- return _XS_ExponentialIntegral($n);
+ return _XS_ExponentialIntegral($n)
+ if !defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' &&
$_Config{'xs'};
+
+ return Math::Prime::Util::PP::ExponentialIntegral($n);
}
sub LogarithmicIntegral {
@@ -1448,17 +1451,15 @@ sub LogarithmicIntegral {
return 0 if $n == 0;
return $_Neg_Infinity if $n == 1;
return $_Infinity if $n == $_Infinity;
- if ($n == 2) {
- return (defined $bignum::VERSION || ref($n) eq 'Math::BigFloat')
- ?
Math::BigFloat->new('1.045163780117492784844588889194613136522615578151201575832909144075013205210359530172717405626383356306')
- : 1.045163780117492784844588889194613136522615578151;
- }
croak("Invalid input to LogarithmicIntegral: x must be >= 0") if $n <= 0;
- return Math::Prime::Util::PP::LogarithmicIntegral($n)
- if defined $bignum::VERSION || ref($n) eq 'Math::BigFloat' ||
!$_Config{'xs'};
- return _XS_LogarithmicIntegral($n);
+ if (!defined $bignum::VERSION && ref($n) ne 'Math::BigFloat' &&
$_Config{'xs'}) {
+ return 1.045163780117492784844588889194613136522615578151 if $n == 2;
+ return _XS_LogarithmicIntegral($n);
+ }
+
+ return Math::Prime::Util::PP::LogarithmicIntegral($n);
}
#############################################################################
@@ -2357,12 +2358,25 @@ find a factor C<p> of C<n> where C<p-1> is smooth (it
has no large factors).
Given a non-zero floating point input C<x>, this returns the real-valued
exponential integral of C<x>, defined as the integral of C<e^t/t dt>
from C<-infinity> to C<x>.
-Depending on the input, the integral is calculated using
+
+If the bignum module has been loaded, all inputs will be treated as if they
+were Math::BigFloat objects.
+
+For non-BigInt/BigFloat objects, the result should be accurate to at least 14
+digits.
+
+For BigInt / BigFloat objects, we first check to see if the Math::MPFR module
+is installed. If so, then it is used, as it will return results much faster
+and can be more accurate. Accuracy when using MPFR will be equal to the
+C<accuracy()> value of the input (or the default BigFloat accuracy, which
+is 40 by default).
+
+MPFR is used for positive inputs only. If Math::MPFR is not installed or the
+input is negative, then other methods are used:
continued fractions (C<x E<lt> -1>),
rational Chebyshev approximation (C< -1 E<lt> x E<lt> 0>),
a convergent series (small positive C<x>),
or an asymptotic divergent series (large positive C<x>).
-
Accuracy should be at least 14 digits.
@@ -2382,11 +2396,20 @@ term C<li0> for this function, and define C<li(x) =
Li0(x) - li0(2)>. Due to
this terminilogy confusion, it is important to check which exact definition is
being used.
+If the bignum module has been loaded, all inputs will be treated as if they
+were Math::BigFloat objects.
-This function is implemented as C<li(x) = Ei(ln x)> after handling special
-values.
+For non-BigInt/BigFloat objects, the result should be accurate to at least 14
+digits.
-Accuracy should be at least 14 digits.
+For BigInt / BigFloat objects, we first check to see if the Math::MPFR module
+is installed. If so, then it is used, as it will return results much faster
+and can be more accurate. Accuracy when using MPFR will be equal to the
+C<accuracy()> value of the input (or the default BigFloat accuracy, which
+is 40 by default).
+
+MPFR is used for inputs greater than 1 only. If Math::MPFR is not installed or
+the input is less than 1, results will be calculated as C<Ei(ln x)>.
=head2 RiemannZeta
@@ -2399,11 +2422,24 @@ subtracted to ensure maximum precision for large values
of C<s>. The zeta
function is the sum from k=1 to infinity of C<1 / k^s>. This function only
uses real arguments, so is basically the Euler Zeta function.
-Accuracy should be at least 14 digits with native numbers and 35 digits with
-bignum or a BigInt/BigFloat argument. Small integer values are returned from
-a table. For native-precision numbers, a rational Chebyshev approximation is
-used between 0.5 and 5, while larger values use a series. Multiple precision
-numbers use Borwein (1991) algorithm 2 or the basic series.
+If the bignum module has been loaded, all inputs will be treated as if they
+were Math::BigFloat objects.
+
+For non-BigInt/BigFloat objects, the result should be accurate to at least 14
+digits. The XS code uses a rational Chebyshev approximation between 0.5 and 5,
+and a series for larger values.
+
+For BigInt / BigFloat objects, we first check to see if the Math::MPFR module
+is installed. If so, then it is used, as it will return results much faster
+and can be more accurate. Accuracy when using MPFR will be equal to the
+C<accuracy()> value of the input (or the default BigFloat accuracy, which
+is 40 by default).
+
+If Math::MPFR is not installed, then results are calculated using either
+Borwein (1991) algorithm 2, or the basic series. Full input accuracy is
+attempted, but there are defects in Math::BigFloat with high accuracy
+computations that make this difficult. It is also very slow. I highly
+recommend installing Math::MPFR for BigFloat computations.
=head2 RiemannR
@@ -2414,8 +2450,18 @@ Given a positive non-zero floating point input, returns
the floating
point value of Riemann's R function. Riemann's R function gives a very close
approximation to the prime counting function.
-Accuracy should be at least 14 digits for native numbers and 35 digits with
-bignum or a BigInt/BigFloat argument.
+If the bignum module has been loaded, all inputs will be treated as if they
+were Math::BigFloat objects.
+
+For non-BigInt/BigFloat objects, the result should be accurate to at least 14
+digits.
+
+For BigInt / BigFloat objects, we first check to see if the Math::MPFR module
+is installed. If so, then it is used, as it will return results much faster
+and can be more accurate. Accuracy when using MPFR will be equal to the
+C<accuracy()> value of the input (or the default BigFloat accuracy, which
+is 40 by default). Accuracy without MPFR should be 35 digits.
+
=head1 EXAMPLES
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 253d15b..5de0e92 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -43,6 +43,8 @@ our $_Infinity = 0+'inf';
$_Infinity = 20**20**20 if 65535 > $_Infinity; # E.g. Windows
our $_Neg_Infinity = -$_Infinity;
+my $_have_MPFR = -1;
+
my $_precalc_size = 0;
sub prime_precalc {
my($n) = @_;
@@ -1436,17 +1438,47 @@ my $_const_euler =
0.57721566490153286060651209008240243104215933593992;
my $_const_li2 = 1.045163780117492784844588889194613136522615578151;
sub ExponentialIntegral {
- my($x, $tol) = @_;
+ my($x) = @_;
return $_Neg_Infinity if $x == 0;
return 0 if $x == $_Neg_Infinity;
return $_Infinity if $x == $_Infinity;
- $tol = 1e-16 unless defined $tol;
- my $sum = 0.0;
- my($y, $t);
- my $c = 0.0;
+
+ # Use MPFR if possible.
+ if ($_have_MPFR < 0) {
+ $_have_MPFR = 1;
+ eval { require Math::MPFR; 1; } or do { $_have_MPFR = 0; };
+ }
+ # Gotcha -- MPFR decided to make negative inputs return NaN. Grrr.
+ if ($_have_MPFR && $x > 0) {
+ my $wantbf = 0;
+ my $xdigits = 17;
+ if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
+ if (!defined $MATH::BigFloat::VERSION) {
+ eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
+ or do { croak "Cannot load Math::BigFloat "; }
+ }
+ $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
+ $wantbf = 1;
+ $xdigits = $x->accuracy || Math::BigFloat->accuracy() ||
Math::BigFloat->div_scale();
+ }
+ 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) : 0.0 + $strval;
+ }
$x = new Math::BigFloat "$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) {
@@ -1518,13 +1550,52 @@ sub LogarithmicIntegral {
return 0 if $x == 0;
return $_Neg_Infinity if $x == 1;
return $_Infinity if $x == $_Infinity;
- return $_const_li2 if $x == 2;
croak "Invalid input to LogarithmicIntegral: x must be > 0" if $x <= 0;
+ # Use MPFR if possible.
+ if ($_have_MPFR < 0) {
+ $_have_MPFR = 1;
+ eval { require Math::MPFR; 1; } or do { $_have_MPFR = 0; };
+ }
+ # Remember MPFR eint doesn't handle negative inputs
+ if ($_have_MPFR && $x >= 1) {
+ my $wantbf = 0;
+ my $xdigits = 17;
+ if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
+ if (!defined $MATH::BigFloat::VERSION) {
+ eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
+ or do { croak "Cannot load Math::BigFloat "; }
+ }
+ $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
+ $wantbf = 1;
+ $xdigits = $x->accuracy || Math::BigFloat->accuracy() ||
Math::BigFloat->div_scale();
+ }
+ $x = log($x); # Use BigFloat to do the log to simplify precision tracking.
+ 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_eint($lix, $rx, $rnd);
+ my $strval = Math::MPFR::Rmpfr_get_str($lix, 10, 0, $rnd);
+ return ($wantbf) ? Math::BigFloat->new($strval) : 0.0 + $strval;
+ }
+
+ if ($x == 2) {
+ my $li2const = (ref($x) eq 'Math::BigFloat') ?
Math::BigFloat->new('1.04516378011749278484458888919461313652261557815120157583290914407501320521')
: $_const_li2;
+ return $li2const;
+ }
+
$x = new Math::BigFloat "$x" if defined $bignum::VERSION && ref($x) ne
'Math::BigFloat';
my $logx = log($x);
# 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.
if ($x > 1e16) {
my $tol = 1e-20;
my $invx = 1.0 / $logx;
@@ -1611,16 +1682,50 @@ my @_Riemann_Zeta_Table = (
sub RiemannZeta {
- my($x, $tol) = @_;
+ my($x) = @_;
+
+ # Use MPFR if possible.
+ if ($_have_MPFR < 0) {
+ $_have_MPFR = 1;
+ eval { require Math::MPFR; 1; } or do { $_have_MPFR = 0; };
+ }
+ if ($_have_MPFR) {
+ my $wantbf = 0;
+ my $xdigits = 17;
+ if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
+ if (!defined $MATH::BigFloat::VERSION) {
+ eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
+ or do { croak "Cannot load Math::BigFloat "; }
+ }
+ $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
+ $wantbf = 1;
+ $xdigits = $x->accuracy || Math::BigFloat->accuracy() ||
Math::BigFloat->div_scale();
+ }
+ 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 $zetax = Math::MPFR->new();
+ # Add more bits to account for the leading zeros.
+ my $extra_bits = int(abs($x));
+ Math::MPFR::Rmpfr_set_prec($zetax, $bit_precision + $extra_bits);
+ 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) : 0.0 + $strval;
+ }
if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
+ # No MPFR, BigFloat
require Math::Prime::Util::ZetaBigFloat;
- return Math::Prime::Util::ZetaBigFloat::RiemannZeta($x, $tol);
+ 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];
- $tol = 1e-16 unless defined $tol;
+ my $tol = 1e-16;
my($y, $t);
my $sum = 0.0;
my $c = 0.0;
@@ -1640,17 +1745,76 @@ sub RiemannZeta {
# Riemann R function
sub RiemannR {
- my($x, $tol) = @_;
+ my($x) = @_;
croak "Invalid input to ReimannR: x must be > 0" if $x <= 0;
+ # Use MPFR if possible.
+ if ($_have_MPFR < 0) {
+ $_have_MPFR = 1;
+ eval { require Math::MPFR; 1; } or do { $_have_MPFR = 0; };
+ }
+ if ($_have_MPFR) {
+ my $wantbf = 0;
+ my $xdigits = 17;
+ if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
+ if (!defined $MATH::BigFloat::VERSION) {
+ eval { require Math::BigFloat; Math::BigFloat->import(); 1; }
+ or do { croak "Cannot load Math::BigFloat "; }
+ }
+ $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
+ $wantbf = 1;
+ $xdigits = $x->accuracy || Math::BigFloat->accuracy() ||
Math::BigFloat->div_scale();
+ }
+ 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 .. 10000) {
+ 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) : 0.0 + $strval;
+ }
+
if (defined $bignum::VERSION || ref($x) =~ /^Math::Big/) {
require Math::Prime::Util::ZetaBigFloat;
- return Math::Prime::Util::ZetaBigFloat::RiemannR($x, $tol);
+ return Math::Prime::Util::ZetaBigFloat::RiemannR($x);
}
- $tol = 1e-16 unless defined $tol;
+ my $tol = 1e-16;
my $sum = 0.0;
my($y, $t);
my $c = 0.0;
@@ -1850,6 +2014,8 @@ math packages. When given two arguments, it returns the
inclusive
count of primes between the ranges (e.g. C<(13,17)> returns 2, C<14,17>
and C<13,16> return 1, and C<14,16> returns 0).
+The Lehmer method is used for large values, which speeds up results greatly.
+
=head2 nth_prime
@@ -1858,8 +2024,9 @@ and C<13,16> return 1, and C<14,16> returns 0).
Returns the prime that lies in index C<n> in the array of prime numbers. Put
another way, this returns the smallest C<p> such that C<Pi(p) E<gt>= n>.
-This relies on generating primes, so can require a lot of time and space for
-large inputs.
+The Lehmer prime count is used to speed up results for large inputs, but both
+methods take quite a bit of time and space. Think abut whether a bound or
+approximation would be acceptable instead.
=head2 is_strong_pseudoprime
@@ -1999,12 +2166,22 @@ others they succeed in a remarkably short time.
Given a non-zero floating point input C<x>, this returns the real-valued
exponential integral of C<x>, defined as the integral of C<e^t/t dt>
from C<-infinity> to C<x>.
-Depending on the input, the integral is calculated using
+
+If the bignum module has been loaded, all inputs will be treated as if they
+were Math::BigFloat objects.
+
+We first check to see if the Math::MPFR module is installed. If so, then
+it is used, as it will return results much faster and can be more accurate.
+Accuracy when using MPFR will be 17 digits for non-BigInt/BigFloats, and
+for BigInt/BigFloat inputs will be equal to the C<accuracy()> value of the
+input (or the default BigFloat accuracy, which is 40 by default).
+
+MPFR is used for positive inputs only. If Math::MPFR is not installed or the
+input is negative, then other methods are used:
continued fractions (C<x E<lt> -1>),
rational Chebyshev approximation (C< -1 E<lt> x E<lt> 0>),
a convergent series (small positive C<x>),
or an asymptotic divergent series (large positive C<x>).
-
Accuracy should be at least 14 digits.
@@ -2021,10 +2198,14 @@ This is often known as C<li(x)>. A related function is
the offset logarithmic
integral, sometimes known as C<Li(x)> which avoids the singularity at 1. It
may be defined as C<Li(x) = li(x) - li(2)>.
-This function is implemented as C<li(x) = Ei(ln x)> after handling special
-values.
+We first check to see if the Math::MPFR module is installed. If so, then
+it is used, as it will return results much faster and can be more accurate.
+Accuracy when using MPFR will be 17 digits for non-BigInt/BigFloats, and
+for BigInt/BigFloat inputs will be equal to the C<accuracy()> value of the
+input (or the default BigFloat accuracy, which is 40 by default).
-Accuracy should be at least 14 digits.
+MPFR is used for inputs greater than 1 only. If Math::MPFR is not installed or
+the input is less than 1, results will be calculated as C<Ei(ln x)>.
=head2 RiemannZeta
@@ -2035,10 +2216,17 @@ point value of ζ(s)-1, where ζ(s) is the Riemann zeta
function. One is
subtracted to ensure maximum precision for large values of C<s>. The zeta
function is the sum from k=1 to infinity of C<1 / k^s>
-Accuracy should be at least 14 digits, but currently does not increase
-accuracy with big floats. Small integer values are returned from a table,
-values between 0.5 and 5 use rational Chebyshev approximation, and larger
-values use a series.
+If the bignum module has been loaded, all inputs will be treated as if they
+were Math::BigFloat objects.
+
+We first check to see if the Math::MPFR module is installed. If so, then
+it is used, as it will return results much faster and can be more accurate.
+Accuracy when using MPFR will be 17 digits for non-BigInt/BigFloats, and
+for BigInt/BigFloat inputs will be equal to the C<accuracy()> value of the
+input (or the default BigFloat accuracy, which is 40 by default).
+
+If Math::MPFR is not installed, then results are calculated either from a
+table, rational Chebyshev approximation, or via a series.
=head2 RiemannR
@@ -2049,7 +2237,14 @@ Given a positive non-zero floating point input, returns
the floating
point value of Riemann's R function. Riemann's R function gives a very close
approximation to the prime counting function.
-Accuracy should be at least 14 digits.
+If the bignum module has been loaded, all inputs will be treated as if they
+were Math::BigFloat objects.
+
+We first check to see if the Math::MPFR module is installed. If so, then
+it is used, as it will return results much faster and can be more accurate.
+Accuracy when using MPFR will be 17 digits for non-BigInt/BigFloats, and
+for BigInt/BigFloat inputs will be equal to the C<accuracy()> value of the
+input (or the default BigFloat accuracy, which is 40 by default).
=head1 LIMITATIONS
diff --git a/lib/Math/Prime/Util/ZetaBigFloat.pm
b/lib/Math/Prime/Util/ZetaBigFloat.pm
index 05aa3d2..54f036f 100644
--- a/lib/Math/Prime/Util/ZetaBigFloat.pm
+++ b/lib/Math/Prime/Util/ZetaBigFloat.pm
@@ -284,12 +284,17 @@ sub _Recompute_Dk {
}
sub RiemannZeta {
- my($x, $tol) = @_;
+ my($x) = @_;
+
+ $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
+ my $xdigits = $x->accuracy || Math::BigFloat->accuracy() ||
Math::BigFloat->div_scale();
return $_Riemann_Zeta_Table[int($x)-2]
- if $x == int($x) && defined $_Riemann_Zeta_Table[int($x)-2];
+ if $x == int($x)
+ && defined $_Riemann_Zeta_Table[int($x)-2]
+ && $xdigits <= 44;
- $tol = 1e-40 unless defined $tol;
+ my $tol = 0.0 + "1e-$xdigits";
# Trying to work around Math::BigFloat bugs RT 43692 and RT 77105 which make
# a right mess of things. Watch this:
@@ -299,7 +304,6 @@ sub RiemannZeta {
# into (6^-(40.5/4))^4 (assuming the base is positive). Without that hack,
# none of this would work at all.
- $x = Math::BigFloat->new("$x");
my $superx = 1;
my $subx = $x->copy;
while ($subx > 8) {
@@ -310,7 +314,8 @@ sub RiemannZeta {
# Go with the basic formula for large x, as it best works around the mess,
# though is unfortunately much slower.
if ($x > 30) {
- my $sum = 0.0;
+ my $sum = Math::BigFloat->bzero;
+ $sum->accuracy($xdigits);
for my $k (4 .. 1000) {
my $term = ( $k ** -$subx ) ** $superx;
$sum += $term;
@@ -335,9 +340,11 @@ sub RiemannZeta {
# return $sum;
#}
- # If we wanted to change the Borwein series being used:
- # _Recompute_Dk(55);
-
+ {
+ my $dig = int($_Borwein_n / 1.3)+1;
+ _Recompute_Dk( int($xdigits * 1.3) + 4 ) if $dig < $xdigits;
+ }
+
if (ref $_Borwein_dk[0] ne 'Math::BigInt') {
@_Borwein_dk = map { Math::BigInt->new("$_") } @_Borwein_dk;
}
@@ -371,23 +378,27 @@ sub RiemannZeta {
$sum += Math::BigFloat->new( $one - $_Borwein_dk[$n] ); # term k=0
$sum->bdiv( $divisor );
$sum->bsub(1);
+ #$sum->fround($xdigits);
return $sum;
}
# Riemann R function
sub RiemannR {
- my($x, $tol) = @_;
+ my($x) = @_;
+ $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
+ my $xdigits = $x->accuracy || Math::BigFloat->accuracy() ||
Math::BigFloat->div_scale();
+ my $tol = 0.0 + "1e-$xdigits";
+
+ # TODO: Changing input accuracy may mean we should recalculate this. Also,
+ # the default table is only 44 digits.
if (scalar @_Riemann_Zeta_Premult == 0) {
@_Riemann_Zeta_Premult = map { my $v = Math::BigFloat->bone;
- $v->accuracy($x->accuracy() || 45);
+ $v->accuracy($xdigits);
$v / ($_Riemann_Zeta_Table[$_-1] * $_ + $_)
}
(1 .. @_Riemann_Zeta_Table);
}
- $x = new Math::BigFloat "$x" if ref($x) ne 'Math::BigFloat';
-
- $tol = 1e-35 unless defined $tol;
my $sum = Math::BigFloat->bone;
my $flogx = log($x);
@@ -396,9 +407,19 @@ sub RiemannR {
for my $k (1 .. 10000) {
my $zeta_term = $_Riemann_Zeta_Premult[$k-1];
if (!defined $zeta_term) {
- my $zeta = (($k-1) <= $#_Riemann_Zeta_Table)
- ? $_Riemann_Zeta_Table[$k-1]
- : RiemannZeta( $k+1 );
+ my $zeta = $_Riemann_Zeta_Table[$k-1];
+ if (!defined $zeta) {
+ my $kz = $k+1;
+ if ($kz >= 100 && $xdigits <= 40) {
+ # For this accuracy level, two terms are more than enough. Also,
+ # we should be able to miss the Math::BigFloat accuracy bug. If we
+ # try to do this for higher accuracy, things will go very bad.
+ $zeta = Math::BigFloat->new(3)->bpow(-$kz)
+ + Math::BigFloat->new(2)->bpow(-$kz);
+ } else {
+ $zeta = Math::Prime::Util::ZetaBigFloat::RiemannZeta( $kz );
+ }
+ }
$zeta_term = Math::BigFloat->bone / ($zeta * $k + $k);
}
$part_term *= $flogx / $k;
@@ -436,12 +457,26 @@ Version 0.14
Math::BigFloat versions`of the Riemann Zeta and Riemann R functions. These
are kept in a separate module because they use a lot of big tables that we'd
-prefer not to have loaded all the time.
+prefer to only load if needed.
=head1 DESCRIPTION
Pure Perl implementations of Riemann Zeta and Riemann R using Math::BigFloat.
+These functions are used if:
+
+=over 4
+
+=item The input is a BigInt, a BigFloat, or the bignum module has been loaded.
+
+=item The Math::MPFR module is not available.
+
+=back
+
+If you use these functions a lot, I B<highly> recommend you install
+L<Math::MPFR>, which the main L<Math::Prime::Util> functions will find.
+These give B<much> better performance, and better accuracy. You can also
+use L<Math::Pari> for the Riemann Zeta function.
=head1 FUNCTIONS
@@ -455,10 +490,9 @@ point value of ζ(s)-1, where ζ(s) is the Riemann zeta
function. One is
subtracted to ensure maximum precision for large values of C<s>. The zeta
function is the sum from k=1 to infinity of C<1 / k^s>
-Accuracy should be at least 14 digits, but currently does not increase
-accuracy with big floats. Small integer values are returned from a table,
-values between 0.5 and 5 use rational Chebyshev approximation, and larger
-values use a series.
+Results are calculated using either Borwein (1991) algorithm 2, or the basic
+series. Full input accuracy is attempted, but there are defects in
+Math::BigFloat with high accuracy computations that make this difficult.
=head2 RiemannR
@@ -469,7 +503,7 @@ Given a positive non-zero floating point input, returns the
floating
point value of Riemann's R function. Riemann's R function gives a very close
approximation to the prime counting function.
-Accuracy should be at least 14 digits.
+Accuracy should be about 35 digits.
=head1 LIMITATIONS
@@ -478,20 +512,20 @@ Bugs in Math::BigFloat (RT 43692, RT 77105) cause many
problems with this code.
I've attempted to apply workarounds, but it is possible there are cases they
miss.
-The accuracy goals (35 digits) are sometimes missed by a digit or two, and
-extensive testing needs to be done to ensure we meet the goals.
+The accuracy goals (35 digits) are sometimes missed by a digit or two.
=head1 PERFORMANCE
-Performance is not good at all. A version using XS+GMP would be good to have.
-Pari can give better accuracy in a miniscule fraction of the time.
+Performance is quite bad.
=head1 SEE ALSO
L<Math::Prime::Util>
+L<Math::MPFR>
+
L<Math::Pari>
diff --git a/util.c b/util.c
index a67b506..9bedf96 100644
--- a/util.c
+++ b/util.c
@@ -846,22 +846,31 @@ long double ld_riemann_zeta(long double x) {
long double sumd =
C8q[0]+x*(C8q[1]+x*(C8q[2]+x*(C8q[3]+x*(C8q[4]+x*(C8q[5]+x*(C8q[6]+x*(C8q[7]+x*C8q[8])))))));
sum = (sumn - (x-1)*sumd) / ((x-1)*sumd);
+ } else if (x > 2000.0) {
+
+ /* 1) zeta(2000)-1 is about 8.7E-603, which is far less than a IEEE-754
+ * 64-bit double can represent. A 128-bit quad could go to ~16000.
+ * 2) pow / powl start getting obnoxiously slow with values like -7500.
+ */
+ return 0.0;
+
} else {
- /* This series needs about half the number of terms as the usual k^-x,
- * and gets slightly better numerical results.
+ /* I was using this series:
* functions.wolfram.com/ZetaFunctionsandPolylogarithms/Zeta/06/04/0003
+ * which for integer values is great -- it needs half the number of terms
+ * and gives slightly better numerical results. However, for non-integer
+ * values using double precision, it's awful.
+ * Back to the defining equation.
*/
- for (k = 2; k <= 1000000; k++) {
- term = powl(2*k+1, -x);
+ for (k = 5; k <= 1000000; k++) {
+ term = powl(k, -x);
KAHAN_SUM(sum, term);
if (term < tol*sum) break;
}
+ KAHAN_SUM(sum, powl(4, -x) );
KAHAN_SUM(sum, powl(3, -x) );
- term = 1.0L / (1.0L - powl(2, -x));
- sum *= term;
- sum += (term - 1.0L);
- /* printf("zeta(%lf) = %.15lf in %d iterations\n", x, sum, k); */
+ KAHAN_SUM(sum, powl(2, -x) );
}
--
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
