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