This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.13 in repository libmath-prime-util-perl.
commit 9c52e2db3823fe073cb148b59c8e89608e288cbd Author: Dana Jacobsen <d...@acm.org> Date: Tue Jul 24 18:47:54 2012 -0600 Add prime_set_config, add assume_rh, use Schoenfeld bounds --- Changes | 3 + lib/Math/Prime/Util.pm | 174 +++++++++++++++++++++++++++++-------------------- 2 files changed, 107 insertions(+), 70 deletions(-) diff --git a/Changes b/Changes index 84002fa..8a38333 100644 --- a/Changes +++ b/Changes @@ -4,6 +4,9 @@ Revision history for Perl extension Math::Prime::Util. - Add RiemannZeta function. We used this before for Riemann R, but now it is exported and should be accurate for small numbers (it was only used for integers > 40 before). + - add prime_set_config(), including ability to set assume_rh to let all + functions assume the Riemann Hypothesis. + - Use the Schoenfeld bound for Pi(x) (x large) if assume_rh is true. 0.11 23 July 2012 - Turn off threading tests on Cygwin, as threads on some Cygwin platforms diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index b1eb9c1..e1b90a4 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -12,7 +12,7 @@ BEGIN { # use parent qw( Exporter ); use base qw( Exporter ); our @EXPORT_OK = qw( - prime_get_config + prime_get_config prime_set_config prime_precalc prime_memfree is_prime is_prob_prime is_strong_pseudoprime is_strong_lucas_pseudoprime @@ -101,6 +101,7 @@ if ($_Config{'maxbits'} == 32) { $_Config{'maxprime'} = 18446744073709551557; $_Config{'maxprimeidx'} = 425656284035217743; } +$_Config{'assume_rh'} = 0; # used for code like: # return _XS_foo($n) if $n <= $_XS_MAXVAL @@ -140,7 +141,28 @@ sub prime_get_config { : Math::Prime::Util::PP::_get_prime_cache_size; return \%config; +} +# Note: You can cause yourself pain if you turn on xs or gmp when they're not +# loaded. Your calls will probably die. +sub prime_set_config { + my %params = (@_); # no defaults + while (my($param, $value) = each %params) { + $param = lc $param; + # dispatch table should go here. + if ($param eq 'xs') { + $_Config{'xs'} = ($value) ? 1 : 0; + $_XS_MAXVAL = $_Config{'xs'} ? $_Config{'maxparam'} : -1; + } elsif ($param eq 'gmp') { + $_Config{'gmp'} = ($value) ? 1 : 0; + $_HAVE_GMP = $_Config{'gmp'}; + } elsif ($param =~ /^(assume[_ ]?)?rh$/ || $param =~ /riemann\s*h/) { + $_Config{'assume_rh'} = ($value) ? 1 : 0; + } else { + croak "Unknown or invalid configuration setting: $param\n"; + } + } + 1; } sub _validate_positive_integer { @@ -940,36 +962,35 @@ sub prime_count_lower { # Rosser & Schoenfeld: x/(logx-1/2) x >= 67 # Dusart 1999: x/logx*(1+1/logx+1.8/logxlogx) x >= 32299 # Dusart 2010: x/logx*(1+1/logx+2.0/logxlogx) x >= 88783 - # The Dusart (1999 or 2010) bounds are far, far better than the others. - # TODO: - # We need a assume_riemann_hypothesis(bool) function, which would let - # these bounds return the Schoenfeld or Stoll limits. The former are - # better for n > ~10^12, the latter for n > ~10^8. - # Given the ability to hand test to ~100_000M, if the Stoll limits are - # better then we can always use them up to the verification point. - - # For smaller numbers this works out well. - return int( $x / ($flogx - 0.7) ) if $x < 599; - - my $a; - # Hand tuned for small numbers (< 60_000M) - if ($x < 2700) { $a = 0.30; } - elsif ($x < 5500) { $a = 0.90; } - elsif ($x < 19400) { $a = 1.30; } - elsif ($x < 32299) { $a = 1.60; } - elsif ($x < 176000) { $a = 1.80; } - elsif ($x < 315000) { $a = 2.10; } - elsif ($x < 1100000) { $a = 2.20; } - elsif ($x < 4500000) { $a = 2.31; } - elsif ($x < 233000000) { $a = 2.36; } - elsif ($x < 5433800000) { $a = 2.32; } - elsif ($x <60000000000) { $a = 2.15; } - else { $a = 2.00; } # Dusart 2010, page 2 - - my $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)); - $result = Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat'; + my $result; + if ($x > 1000_000_000_000 && $_Config{'assume_rh'}) { + my $lix = LogarithmicIntegral($x); + my $sqx = sqrt($x); + # Schoenfeld bound: (constant is 8 * Pi) + $result = $lix - (($sqx*$flogx) / 25.13274122871834590770114707); + } elsif ($x < 599) { + $result = $x / ($flogx - 0.7); # For smaller numbers this works out well. + } else { + my $a; + # Hand tuned for small numbers (< 60_000M) + if ($x < 2700) { $a = 0.30; } + elsif ($x < 5500) { $a = 0.90; } + elsif ($x < 19400) { $a = 1.30; } + elsif ($x < 32299) { $a = 1.60; } + elsif ($x < 176000) { $a = 1.80; } + elsif ($x < 315000) { $a = 2.10; } + elsif ($x < 1100000) { $a = 2.20; } + elsif ($x < 4500000) { $a = 2.31; } + elsif ($x < 233000000) { $a = 2.36; } + elsif ($x < 5433800000) { $a = 2.32; } + elsif ($x <60000000000) { $a = 2.15; } + else { $a = 2.00; } # Dusart 2010, page 2 + $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)); + } + + return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat'; return int($result); } @@ -993,47 +1014,53 @@ sub prime_count_upper { my $flogx = log($x); - # These work out well for small values - return int( ($x / ($flogx - 1.048)) + 1.0 ) if $x < 1621; - return int( ($x / ($flogx - 1.071)) + 1.0 ) if $x < 5000; - return int( ($x / ($flogx - 1.098)) + 1.0 ) if $x < 15900; - - my $a; - # Hand tuned for small numbers (< 60_000M) - if ($x < 24000) { $a = 2.30; } - elsif ($x < 59000) { $a = 2.48; } - elsif ($x < 350000) { $a = 2.52; } - elsif ($x < 355991) { $a = 2.54; } - elsif ($x < 356000) { $a = 2.51; } - elsif ($x < 3550000) { $a = 2.50; } - elsif ($x < 3560000) { $a = 2.49; } - elsif ($x < 5000000) { $a = 2.48; } - elsif ($x < 8000000) { $a = 2.47; } - elsif ($x < 13000000) { $a = 2.46; } - elsif ($x < 18000000) { $a = 2.45; } - elsif ($x < 31000000) { $a = 2.44; } - elsif ($x < 41000000) { $a = 2.43; } - elsif ($x < 48000000) { $a = 2.42; } - elsif ($x < 119000000) { $a = 2.41; } - elsif ($x < 182000000) { $a = 2.40; } - elsif ($x < 192000000) { $a = 2.395; } - elsif ($x < 213000000) { $a = 2.390; } - elsif ($x < 271000000) { $a = 2.385; } - elsif ($x < 322000000) { $a = 2.380; } - elsif ($x < 400000000) { $a = 2.375; } - elsif ($x < 510000000) { $a = 2.370; } - elsif ($x < 682000000) { $a = 2.367; } - elsif ($x < 2953652287) { $a = 2.362; } - else { $a = 2.334; } # Dusart 2010, page 2 - #elsif ($x <60000000000) { $a = 2.362; } - #else { $a = 2.51; } # Dusart 1999, page 14 - - # Old versions of Math::BigFloat will do the Wrong Thing with this. - #return int( ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0 ); - my $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0; + my $result; + if ($x > 10000_000_000_000 && $_Config{'assume_rh'}) { + my $lix = LogarithmicIntegral($x); + my $sqx = sqrt($x); + # Schoenfeld bound: (constant is 8 * Pi) + $result = $lix + (($sqx*$flogx) / 25.13274122871834590770114707); + } elsif ($x < 1621) { $result = ($x / ($flogx - 1.048)) + 1.0; } + elsif ($x < 5000) { $result = ($x / ($flogx - 1.071)) + 1.0; } + elsif ($x < 15900) { $result = ($x / ($flogx - 1.098)) + 1.0; } + else { + my $a; + # Hand tuned for small numbers (< 60_000M) + if ($x < 24000) { $a = 2.30; } + elsif ($x < 59000) { $a = 2.48; } + elsif ($x < 350000) { $a = 2.52; } + elsif ($x < 355991) { $a = 2.54; } + elsif ($x < 356000) { $a = 2.51; } + elsif ($x < 3550000) { $a = 2.50; } + elsif ($x < 3560000) { $a = 2.49; } + elsif ($x < 5000000) { $a = 2.48; } + elsif ($x < 8000000) { $a = 2.47; } + elsif ($x < 13000000) { $a = 2.46; } + elsif ($x < 18000000) { $a = 2.45; } + elsif ($x < 31000000) { $a = 2.44; } + elsif ($x < 41000000) { $a = 2.43; } + elsif ($x < 48000000) { $a = 2.42; } + elsif ($x < 119000000) { $a = 2.41; } + elsif ($x < 182000000) { $a = 2.40; } + elsif ($x < 192000000) { $a = 2.395; } + elsif ($x < 213000000) { $a = 2.390; } + elsif ($x < 271000000) { $a = 2.385; } + elsif ($x < 322000000) { $a = 2.380; } + elsif ($x < 400000000) { $a = 2.375; } + elsif ($x < 510000000) { $a = 2.370; } + elsif ($x < 682000000) { $a = 2.367; } + elsif ($x < 2953652287) { $a = 2.362; } + else { $a = 2.334; } # Dusart 2010, page 2 + #elsif ($x <60000000000) { $a = 2.362; } + #else { $a = 2.51; } # Dusart 1999, page 14 + + # Old versions of Math::BigFloat will do the Wrong Thing with this. + #return int( ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0 ); + $result = ($x/$flogx) * (1.0 + 1.0/$flogx + $a/($flogx*$flogx)) + 1.0; + } + return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat'; return int($result); - } ############################################################################# @@ -1508,7 +1535,9 @@ use the Dusart (2010) bounds of x/logx * (1 + 1/logx + 2.334/log^2x) >= Pi(x) -above that range. These bounds do not assume the Riemann Hypothesis. +above that range. These bounds do not assume the Riemann Hypothesis. If the +configuration option C<assume_rh> has been set (it is off by default), then +the Schoenfeld (1976) bounds are used for large values. =head2 prime_count_approx @@ -1958,7 +1987,10 @@ Accuracy should be at least 14 digits. Given a floating point input C<s> where C<s E<gt>= 0.5>, returns the floating 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> +function is the sum from k=1 to infinity of C<1 / k^s>. + +Since the argument and the result are real, this is technically the Euler Zeta +function. Accuracy should be at least 14 digits, but currently does not increase accuracy with big floats. Small integer values are returned from a table, @@ -2181,6 +2213,8 @@ excellent versions in the public domain). =item Pierre-Alain Fouque and Mehdi Tibouchi, "Close to Uniform Prime Number Generation With Fewer Random Bits", 2011. L<http://eprint.iacr.org/2011/481> +=item Douglas A. Stoll and Patrick Demichel , "The impact of ζ(s) complex zeros on π(x) for x E<lt> 10^{10^{13}}", Mathematics of Computation, v80, n276, pp 2381-2394, October 2011. L<http://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4/home.html> + =back -- 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