This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.23
in repository libmath-prime-util-perl.
commit 1d1d03d520817e1787beef1469a53f37e6267bd0
Author: Dana Jacobsen <d...@acm.org>
Date: Wed Feb 27 16:52:00 2013 -0800
Update PP Zeta -- much better now
---
Changes | 7 ++++--
lib/Math/Prime/Util/PP.pm | 57 ++++++++++++++++++++++++++++++++++-------------
util.c | 6 ++---
3 files changed, 49 insertions(+), 21 deletions(-)
diff --git a/Changes b/Changes
index 120c57c..c2f25eb 100644
--- a/Changes
+++ b/Changes
@@ -5,8 +5,11 @@ Revision history for Perl extension Math::Prime::Util.
- exp_mangoldt in XS, and speed up the PP version.
- Replace XS Zeta for x > 5 with series from Cephes. It is 1 eps more
- accurate for a small fracion of inputs. More importantly, it is much
- faster in range 5 < x < 10.
+ accurate for a small fraction of inputs. More importantly, it is much
+ faster in range 5 < x < 10. This only affects non-integer inputs.
+
+ - PP Zeta code replaced (for no-MPFR, non-bignums) with new series. The
+ new code is much more accurate for small values, and *much* faster.
0.22 26 February 2013
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 4f036d8..feeb506 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1744,22 +1744,47 @@ sub RiemannZeta {
# No MPFR, no BigFloat.
return 0.0 + $_Riemann_Zeta_Table[int($x)-2]
if $x == int($x) && defined $_Riemann_Zeta_Table[int($x)-2];
- my $tol = 1e-16;
- my($y, $t);
- my $sum = 0.0;
- my $c = 0.0;
-
- for my $k (2 .. 1000000) {
- my $term = (2*$k+1) ** -$x;
- $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
- last if $term < abs($tol*$sum);
- }
- my $term = 3 ** -$x;
- $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
- $t = 1.0 / (1.0 - (2 ** -$x));
- $sum *= $t;
- $sum += ($t - 1.0);
- return $sum;
+ my $tol = 1.11e-16;
+
+ # Series based on (2n)! / B_2n.
+ # This is a simplication of the Cephes zeta function.
+ my @A = (
+ 12.0,
+ -720.0,
+ 30240.0,
+ -1209600.0,
+ 47900160.0,
+ -1892437580.3183791606367583212735166426,
+ 74724249600.0,
+ -2950130727918.1642244954382084600497650,
+ 116467828143500.67248729113000661089202,
+ -4597978722407472.6105457273596737891657,
+ 181521054019435467.73425331153534235290,
+ -7166165256175667011.3346447367083352776,
+ 282908877253042996618.18640556532523927,
+ );
+ my $s = 0.0;
+ my $b = 0.0;
+ foreach my $i (2 .. 10) {
+ $b = $i ** -$x;
+ $s += $b;
+ return $s if abs($b/$s) < $tol;
+ }
+ my $w = 10.0;
+ $s = $s + $b*$w/($x-1.0) - 0.5*$b;
+ my $a = 1.0;
+ foreach my $i (0 .. 12) {
+ my $k = 2*$i;
+ $a *= $x + $k;
+ $b /= $w;
+ my $t = $a*$b/$A[$i];
+ $s += $t;
+ $t = abs($t/$s);
+ last if $t < $tol;
+ $a *= $x + $k + 1.0;
+ $b /= $w;
+ }
+ return $s;
}
# Riemann R function
diff --git a/util.c b/util.c
index b2056b2..fe2fc19 100644
--- a/util.c
+++ b/util.c
@@ -967,8 +967,9 @@ long double ld_riemann_zeta(long double x) {
-7166165256175667011.3346447367083352776L,
282908877253042996618.18640556532523927L,
};
- long double a, b, s, t, w;
- s = powl(1.0, -x) - 1.0;
+ long double a, b, s, t;
+ const long double w = 10.0;
+ s = 0.0;
b = 0.0;
for (i = 2; i < 11; i++) {
b = powl( i, -x );
@@ -976,7 +977,6 @@ long double ld_riemann_zeta(long double x) {
if (fabsl(b/s) < tol)
return s;
}
- w = i-1;
s = s + b*w/(x-1.0) - 0.5 * b;
a = 1.0;
for (i = 0; i < 13; i++) {
--
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
