[libmath-prime-util-perl] 05/25: Add commentary to measure_zeta_accuracy test

[Partha P. Mukherjee]

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

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

commit 622fa8f6a424b06793fad61b423f3adc4ad5835c
Author: Dana Jacobsen <dana at acm.org>
Date:   Thu Nov 21 17:38:13 2013 -0800

    Add commentary to measure_zeta_accuracy test
---
 TODO                                |  7 +++----
 lib/Math/Prime/Util/PP.pm           |  2 +-
 lib/Math/Prime/Util/ZetaBigFloat.pm | 11 +++++------
 xt/measure_zeta_accuracy.pl         | 20 +++++++++++++++++---
 xt/small-is-next-prev.pl            |  0
 5 files changed, 26 insertions(+), 14 deletions(-)

diff --git a/TODO b/TODO
index b336fdd..bc17e82 100644
--- a/TODO
+++ b/TODO
@@ -22,6 +22,8 @@
 
 - Test all routines for numbers on word-size boundary, or ranges that cross.
 
+- More testing on 32-bit machines.
+
 - An assembler version of mulmod for i386 would be _really_ helpful for
   all the non-x86-64 Intel machines.
 
@@ -57,7 +59,4 @@
   algorithm).  The PP code isn't doing that, which means we're doing lots of
   extra primality checks, which aren't cheap in PP.
 
-- use Math::BigInt instead of require, and return bigints as needed.
-  This is a fundamental change in how some functions operate, though likely
-  one that most people wouldn't see, and should be less surprise:
-  e.g. next_prime(~0) returns 18446744073709551629 rather than 0.
+- Consider using Test::Number::Delta for many tests
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index dea5f81..add2d00 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -2310,7 +2310,7 @@ sub RiemannZeta {
       $xdigits = $x->accuracy || Math::BigFloat->accuracy() || Math::BigFloat->div_scale();
     }
     my $rnd = 0;  # MPFR_RNDN;
-    my $bit_precision = int($xdigits * 3.322) + 5;
+    my $bit_precision = int($xdigits * 3.322) + 7;
     my $rx = Math::MPFR->new();
     Math::MPFR::Rmpfr_set_prec($rx, $bit_precision);
     Math::MPFR::Rmpfr_set_str($rx, "$x", 10, $rnd);
diff --git a/lib/Math/Prime/Util/ZetaBigFloat.pm b/lib/Math/Prime/Util/ZetaBigFloat.pm
index d837fcd..2997646 100644
--- a/lib/Math/Prime/Util/ZetaBigFloat.pm
+++ b/lib/Math/Prime/Util/ZetaBigFloat.pm
@@ -374,19 +374,18 @@ sub RiemannZeta {
   my $divisor = $one->copy->bsub($d1)->bmul(-$_Borwein_dk[$n]);
   $tol = $divisor->copy->bmul($tol)->babs();
 
-  my $sum = $zero->copy;
+  my ($sum, $bigk) = ($zero->copy, $one->copy);
   foreach my $k (1 .. $n-1) {
+    my $den = $bigk->binc()->copy->bpow($subx,$xdigits)->bpow($superx,$xdigits);
     my $term = ($k % 2)
              ? $zero->copy->badd($_Borwein_dk[$n])->bsub($_Borwein_dk[$k])
              : $zero->copy->badd($_Borwein_dk[$k])->bsub($_Borwein_dk[$n]);
-    my $den = $zero->copy->badd($k+1);
-    $den = ($den ** $subx) ** $superx;
-    $term /= $den;
-    $sum += $term;
+    $term->bdiv($den);
+    $sum->badd($term);
     last if $term->copy->babs() < $tol;
   }
   $sum->badd($one->copy->bsub($_Borwein_dk[$n]));  # term k=0
-  $sum->bdiv( $divisor )->bdec;
+  $sum->bdiv($divisor,$xdigits)->bdec;
   $sum->bround($xdigits-$extra_acc);
   return $sum;
 }
diff --git a/xt/measure_zeta_accuracy.pl b/xt/measure_zeta_accuracy.pl
old mode 100644
new mode 100755
index 0d50f96..ba5c1b3
--- a/xt/measure_zeta_accuracy.pl
+++ b/xt/measure_zeta_accuracy.pl
@@ -6,7 +6,8 @@ use Math::BigFloat lib=>'GMP';
 use Math::Prime::Util qw/:all/;
 use Term::ANSIColor;
 
-my $acc = 40;
+my $acc = shift || 40;
+die "Max accuracy for this test = 130 digits\n" if $acc > 130;
 
 # gp
 # \p 200
@@ -23,12 +24,25 @@ my %rvals = (
   '200' => '0.0000000000000000000000000000000000000000000000000000000000006223015277861141707144064053780124278238871664711431331935339387492776093057166188727575094880097645495454472391197851568776550275806071517',
 );
 
+my $acctext = ($acc == 40) ? "default 40-digit" : "$acc-digit";
+print <<EOT;
+Using $acctext accuracy.
+
+The first number is the precalculated correct value.
+The second number is the answer RiemannZeta is giving.
+Differences are highlighted in red.
+
+Either Math::MPFR or the Math::BigInt patch from RT43692
+are needed to prevent differences for accuracy > 38 digits.
+EOT
+
+
 foreach my $vstr (sort { $a <=> $b } keys %rvals) {
   my $zeta_str = $rvals{$vstr};
   my $lead = index($zeta_str, '.');
   my $v    = Math::BigFloat->new($vstr);
   my $zeta = Math::BigFloat->new($rvals{$vstr});
-  $v->accuracy($acc) if defined $acc && $acc != 40;
+  $v->accuracy($acc) if $acc != 40;
   #print "zeta($v) = $zeta\n";
   my $mpuzeta = RiemannZeta($v);
   my $mpuzeta_str = ref($mpuzeta) eq 'Math::BigFloat'
@@ -37,7 +51,7 @@ foreach my $vstr (sort { $a <=> $b } keys %rvals) {
   my $mzlen = length($mpuzeta_str);
   # Truncate zeta_str to length of mpuzeta_str, with rounding.
   {
-    $zeta_str = Math::BigFloat->new($zeta_str)->bmul(1,$acc||40)->bstr;
+    $zeta_str = Math::BigFloat->new($zeta_str)->bmul(1,$acc)->bstr;
   }
   if ($zeta_str ne $mpuzeta_str) {
     my $n = 0;
diff --git a/xt/small-is-next-prev.pl b/xt/small-is-next-prev.pl
old mode 100644
new mode 100755

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