This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.37 in repository libmath-prime-util-perl.
commit f2d43c34330f108b3330e09a6688406bf7b35f0e Author: Dana Jacobsen <d...@acm.org> Date: Tue Jan 21 15:40:27 2014 -0800 Move some functions from Util.pm to util.c. Faster, and reduces startup bloat. --- Changes | 14 +- XS.xs | 33 +++- bin/factor.pl | 5 +- examples/test-factor-gnufactor.pl | 2 +- lib/Math/Prime/Util.pm | 337 +------------------------------------- lib/Math/Prime/Util/PP.pm | 317 ++++++++++++++++++++++++++++++++++- lmo.c | 6 +- util.c | 161 ++++++++++++++++-- util.h | 8 +- 9 files changed, 522 insertions(+), 361 deletions(-) diff --git a/Changes b/Changes index d81e0b2..f0b0955 100644 --- a/Changes +++ b/Changes @@ -10,11 +10,21 @@ Revision history for Perl module Math::Prime::Util - Dynamically loads the PP code and Math::BigInt only when needed. This removes a lot of bloat for the usual cases: 2.0 MB perl -E 'say 1' - 4.2 MB Math::Prime::XS - 4.6 MB MPU 0.37 + 4.4 MB MPU 0.37 + 4.5 MB Math::Prime::XS + Math::Factor::XS + 5.3 MB Math::Pari 9.6 MB MPU 0.36 10.0 MB MPU 0.35 + - Together with the previous, adjusted factor script. Over 2x faster + with native integer (non-expression) arguments. This is just not + loading thousands of lines of Perl code that aren't used, which + was more time-consuming than the actual factoring. + + - nth_prime_{lower,upper,approx} and prime_count_{lower,upper,approx} + moved to XS->PP. This helps us slim down and cut startup overhead. + + 0.36 2014-01-13 [API Changes] diff --git a/XS.xs b/XS.xs index cb234f6..d148121 100644 --- a/XS.xs +++ b/XS.xs @@ -589,27 +589,46 @@ next_prime(IN SV* svn) ALIAS: prev_prime = 1 nth_prime = 2 + nth_prime_upper = 3 + nth_prime_lower = 4 + nth_prime_approx = 5 + prime_count_upper = 6 + prime_count_lower = 7 + prime_count_approx = 8 PPCODE: if (_validate_int(aTHX_ svn, 0)) { UV n = my_svuv(svn); - if ( ((ix == 0) && (n >= MPU_MAX_PRIME)) || - ((ix == 2) && (n >= MPU_MAX_PRIME_IDX)) ) { + if ( (n >= MPU_MAX_PRIME && ix == 0) || + (n >= MPU_MAX_PRIME_IDX && (ix==2 || ix==3 || ix==4 || ix==5)) ) { /* Out of range. Fall through to Perl. */ } else { UV ret; switch (ix) { case 0: ret = next_prime(n); break; case 1: ret = (n < 3) ? 0 : prev_prime(n); break; - case 2: - default:ret = _XS_nth_prime(n); break; + case 2: ret = nth_prime(n); break; + case 3: ret = nth_prime_upper(n); break; + case 4: ret = nth_prime_lower(n); break; + case 5: ret = nth_prime_approx(n); break; + case 6: ret = prime_count_upper(n); break; + case 7: ret = prime_count_lower(n); break; + case 8: + default:ret = prime_count_approx(n); break; } XSRETURN_UV(ret); } } switch (ix) { - case 0: _vcallsub("_generic_next_prime"); break; - case 1: _vcallsub("_generic_prev_prime"); break; - default: _vcallsub_with_pp("nth_prime"); break; + case 0: _vcallsub("_generic_next_prime"); break; + case 1: _vcallsub("_generic_prev_prime"); break; + case 2: _vcallsub_with_pp("nth_prime"); break; + case 3: _vcallsub_with_pp("nth_prime_upper"); break; + case 4: _vcallsub_with_pp("nth_prime_lower"); break; + case 5: _vcallsub_with_pp("nth_prime_approx"); break; + case 6: _vcallsub_with_pp("prime_count_upper"); break; + case 7: _vcallsub_with_pp("prime_count_lower"); break; + case 8: + default: _vcallsub_with_pp("prime_count_approx"); break; } return; /* skip implicit PUTBACK */ diff --git a/bin/factor.pl b/bin/factor.pl index 6fca37b..779c2db 100755 --- a/bin/factor.pl +++ b/bin/factor.pl @@ -2,10 +2,8 @@ use strict; use warnings; use Getopt::Long; -use bigint try => 'GMP'; use Math::Prime::Util qw/factor nth_prime prime_set_config/; $| = 1; -no bigint; # Allow execution of any of these functions in the command line my @mpu_funcs = (qw/next_prime prev_prime prime_count nth_prime random_prime @@ -23,7 +21,7 @@ GetOptions(\%opts, ) || die_usage(); if (exists $opts{'version'}) { my $version_str = - "factor.pl version 1.1 using Math::Prime::Util $Math::Prime::Util::VERSION"; + "factor.pl version 1.2 using Math::Prime::Util $Math::Prime::Util::VERSION"; $version_str .= " and MPU::GMP $Math::Prime::Util::GMP::VERSION" if Math::Prime::Util::prime_get_config->{'gmp'}; $version_str .= "\nWritten by Dana Jacobsen.\n"; @@ -69,6 +67,7 @@ sub eval_expr { $expr =~ s/:$mpu_func_map{$func}\(/Math::Prime::Util::$func(/g; } $expr =~ s/(\d+)/ Math::BigInt->new("$1") /g; + $expr = 'use Math::BigInt try=>"GMP"; ' . $expr; my $res = eval $expr; ## no critic die "Cannot eval: $expr\n" if !defined $res; $res = int($res->bstr) if ref($res) eq 'Math::BigInt' && $res <= ~0; diff --git a/examples/test-factor-gnufactor.pl b/examples/test-factor-gnufactor.pl index 229461e..f5614ba 100755 --- a/examples/test-factor-gnufactor.pl +++ b/examples/test-factor-gnufactor.pl @@ -150,7 +150,7 @@ sub mpu_factors { my @ns = @_; my $numpercommand = int( (4000-30)/(length($ns[-1])+1) ); while (@ns) { - my $cs = join(" ", 'factor.pl', splice(@ns, 0, $numpercommand)); + my $cs = join(" ", 'perl -Iblib/lib -Iblib/arch bin/factor.pl', splice(@ns, 0, $numpercommand)); my $fout = qx{$cs}; my @flines = split(/\n/, $fout); foreach my $fline (@flines) { diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm index d2f5ec6..bcd2488 100644 --- a/lib/Math/Prime/Util.pm +++ b/lib/Math/Prime/Util.pm @@ -114,6 +114,12 @@ BEGIN { *mertens = \&Math::Prime::Util::_generic_mertens; *prime_count = \&Math::Prime::Util::_generic_prime_count; *nth_prime = \&Math::Prime::Util::PP::nth_prime; + *nth_prime_upper=\&Math::Prime::Util::PP::nth_prime_upper; + *nth_prime_lower=\&Math::Prime::Util::PP::nth_prime_lower; + *nth_prime_approx=\&Math::Prime::Util::PP::nth_prime_approx; + *prime_count_upper=\&Math::Prime::Util::PP::prime_count_upper; + *prime_count_lower=\&Math::Prime::Util::PP::prime_count_lower; + *prime_count_approx=\&Math::Prime::Util::PP::prime_count_approx; *carmichael_lambda = \&Math::Prime::Util::_generic_carmichael_lambda; *kronecker = \&Math::Prime::Util::_generic_kronecker; *divisor_sum = \&Math::Prime::Util::_generic_divisor_sum; @@ -286,29 +292,6 @@ sub _upgrade_to_float { return Math::BigFloat->new($_[0]); } -# TODO: Remove these, push more funcs into XS -my @_primes_small = ( - 0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97, - 101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191, - 193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283, - 293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401, - 409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509); -sub _tiny_prime_count { - my($n) = @_; - return if $n >= $_primes_small[-1]; - my $j = $#_primes_small; - my $i = 1 + ($n >> 4); - while ($i < $j) { - my $mid = ($i+$j)>>1; - if ($_primes_small[$mid] <= $n) { $i = $mid+1; } - else { $j = $mid; } - } - return $i-1; -} - - - - ############################################################################# @@ -822,36 +805,8 @@ sub _generic_divisors { my($n) = @_; _validate_num($n) || _validate_positive_integer($n); - # In scalar context, returns sigma_0(n). Very fast. - return divisor_sum($n,0) unless wantarray; - return ($n == 0) ? (0,1) : (1) if $n <= 1; - - my %all_factors; - my @factors = factor($n); - return (1,$n) if scalar @factors == 1; - - if (ref($n) eq 'Math::BigInt') { - foreach my $f1 (@factors) { - my $big_f1 = Math::BigInt->new("$f1"); - my @to_add = map { ($_ <= ''.~0) ? _bigint_to_int($_) : $_ } - grep { $_ < $n } - map { $big_f1 * $_ } - keys %all_factors; - undef @all_factors{ $f1, @to_add }; - } - } else { - foreach my $f1 (@factors) { - my @to_add = grep { $_ < $n } - map { $f1 * $_ } - keys %all_factors; - undef @all_factors{ $f1, @to_add }; - } - } - # Add 1 and n - undef $all_factors{1}; - undef $all_factors{$n}; - my @divisors = sort {$a<=>$b} keys %all_factors; - return @divisors; + require Math::Prime::Util::PP; + return Math::Prime::Util::PP::divisors($n); } @@ -1010,282 +965,6 @@ sub verify_prime { ############################################################################# -sub prime_count_approx { - my($x) = @_; - _validate_num($x) || _validate_positive_integer($x); - - # With XS using 30k tables, this is super fast. - return prime_count($x) if $x < $_XS_MAXVAL && $x < 3_000_000; - # Give an exact answer for what we have in our little table. - return _tiny_prime_count($x) if $x < $_primes_small[-1]; - - # Below 2^58th or so, all differences between the high precision and C double - # precision result are less than 0.5. - if ($x <= $_XS_MAXVAL && $x <= 144115188075855872) { - return int(_XS_RiemannR($x) + 0.5); - } - - # Turn on high precision FP if they gave us a big number. - $x = _upgrade_to_float($x) if ref($_[0]) eq 'Math::BigInt'; - - # Method 10^10 %error 10^19 %error - # ----------------- ------------ ------------ - # n/(log(n)-1) .22% .06% - # average bounds .01% .0002% - # li(n) .0007% .00000004% - # li(n)-li(n^.5)/2 .0004% .00000001% - # R(n) .0004% .00000001% - # - # Also consider: http://trac.sagemath.org/sage_trac/ticket/8135 - - # my $result = int( (prime_count_upper($x) + prime_count_lower($x)) / 2); - # my $result = int( LogarithmicIntegral($x) ); - # my $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2); - # my $result = RiemannR($x) + 0.5; - - # Sadly my Perl RiemannR function is really slow for big values. If MPFR - # is available, then use it -- it rocks. Otherwise, switch to LiCorr for - # very big values. This is hacky and shouldn't be necessary. - my $result; - if ( $x < 1e36 || Math::Prime::Util::PP::_MPFR_available() ) { - 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); - } - $result = RiemannR($x) + 0.5; - } else { - $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2); - } - - return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat'; - return int($result); -} - -sub prime_count_lower { - my($x) = @_; - _validate_num($x) || _validate_positive_integer($x); - - # With XS using 30k tables, this is super fast. - return prime_count($x) if $x < $_XS_MAXVAL && $x < 3_000_000; - # Give an exact answer for what we have in our little table. - return _tiny_prime_count($x) if $x < $_primes_small[-1]; - - $x = _upgrade_to_float($x) - if ref($x) eq 'Math::BigInt' || ref($_[0]) eq 'Math::BigInt'; - - my $flogx = log($x); - - # Chebyshev: 1*x/logx x >= 17 - # 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. - - 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); -} - -sub prime_count_upper { - my($x) = @_; - _validate_num($x) || _validate_positive_integer($x); - - # With XS using 30k tables, this is super fast. - return prime_count($x) if $x < $_XS_MAXVAL && $x < 3_000_000; - # Give an exact answer for what we have in our little table. - return _tiny_prime_count($x) if $x < $_primes_small[-1]; - - $x = _upgrade_to_float($x) - if ref($x) eq 'Math::BigInt' || ref($_[0]) eq 'Math::BigInt'; - - # Chebyshev: 1.25506*x/logx x >= 17 - # Rosser & Schoenfeld: x/(logx-3/2) x >= 67 - # Dusart 1999: x/logx*(1+1/logx+2.51/logxlogx) x >= 355991 - # Dusart 2010: x/logx*(1+1/logx+2.334/logxlogx) x >= 2_953_652_287 - - # As with the lower bounds, Dusart bounds are best by far. - - # Another possibility here for numbers under 3000M is to use Li(x) - # minus a correction. - - my $flogx = log($x); - - 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. - $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); -} - -############################################################################# - -sub nth_prime_approx { - my($n) = @_; - _validate_num($n) || _validate_positive_integer($n); - - return $_primes_small[$n] if $n <= $#_primes_small; - - $n = _upgrade_to_float($n) - if ref($n) eq 'Math::BigInt' || $n >= MPU_MAXPRIMEIDX; - - my $flogn = log($n); - my $flog2n = log($flogn); - - # Cipolla 1902: - # m=0 fn * ( flogn + flog2n - 1 ); - # m=1 + ((flog2n - 2)/flogn) ); - # m=2 - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn)) - # + O((flog2n/flogn)^3) - # - # Shown in Dusart 1999 page 12, as well as other sources such as: - # http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf - # where the main issue you run into is that you're doing polynomial - # interpolation, so it oscillates like crazy with many high-order terms. - # Hence I'm leaving it at m=2. - # - - my $approx = $n * ( $flogn + $flog2n - 1 - + (($flog2n - 2)/$flogn) - - ((($flog2n*$flog2n) - 6*$flog2n + 11) / (2*$flogn*$flogn)) - ); - - # Apply a correction to help keep values close. - my $order = $flog2n/$flogn; - $order = $order*$order*$order * $n; - - if ($n < 259) { $approx += 10.4 * $order; } - elsif ($n < 775) { $approx += 7.52* $order; } - elsif ($n < 1271) { $approx += 5.6 * $order; } - elsif ($n < 2000) { $approx += 5.2 * $order; } - elsif ($n < 4000) { $approx += 4.3 * $order; } - elsif ($n < 12000) { $approx += 3.0 * $order; } - elsif ($n < 150000) { $approx += 2.1 * $order; } - elsif ($n < 200000000) { $approx += 0.0 * $order; } - else { $approx += -0.010 * $order; } - # $approx = -0.025 is better for the last, but it gives problems with some - # other code that always wants the asymptotic approximation to be >= actual. - - return int($approx + 0.5); -} - -# The nth prime will be greater than or equal to this number -sub nth_prime_lower { - my($n) = @_; - _validate_num($n) || _validate_positive_integer($n); - - return $_primes_small[$n] if $n <= $#_primes_small; - - $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45; - - my $flogn = log($n); - my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n) - - # Dusart 1999 page 14, for all n >= 2 - #my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.25)/$flogn)); - # Dusart 2010 page 2, for all n >= 3 - my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.10)/$flogn)); - - return int($lower); -} - -# The nth prime will be less or equal to this number -sub nth_prime_upper { - my($n) = @_; - _validate_num($n) || _validate_positive_integer($n); - - return $_primes_small[$n] if $n <= $#_primes_small; - - $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45; - - my $flogn = log($n); - my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n) - - my $upper; - if ($n >= 688383) { # Dusart 2010 page 2 - $upper = $n * ( $flogn + $flog2n - 1.0 + (($flog2n-2.00)/$flogn) ); - } elsif ($n >= 178974) { # Dusart 2010 page 7 - $upper = $n * ( $flogn + $flog2n - 1.0 + (($flog2n-1.95)/$flogn) ); - } elsif ($n >= 39017) { # Dusart 1999 page 14 - $upper = $n * ( $flogn + $flog2n - 0.9484 ); - } elsif ($n >= 6) { # Modified Robin 1983, for 6-39016 only - $upper = $n * ( $flogn + 0.6000 * $flog2n ); - } else { - $upper = $n * ( $flogn + $flog2n ); - } - - return int($upper + 1.0); -} - - -############################################################################# - - ############################################################################# sub RiemannZeta { diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm index e72ae6d..c54d71e 100644 --- a/lib/Math/Prime/Util/PP.pm +++ b/lib/Math/Prime/Util/PP.pm @@ -76,6 +76,12 @@ sub _bigint_to_int { : int("$_[0]"); } +sub _upgrade_to_float { + do { require Math::BigFloat; Math::BigFloat->import(); } + if !defined $Math::BigFloat::VERSION; + return Math::BigFloat->new($_[0]); +} + sub _validate_num { my($n, $min, $max) = @_; croak "Parameter must be defined" if !defined $n; @@ -124,11 +130,7 @@ my @_primes_small = ( 101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191, 193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283, 293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401, - 409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499); -my @_prime_count_small = ( - 0,0,1,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,8,8,8,8,9,9,9,9,9,9,10,10, - 11,11,11,11,11,11,12,12,12,12,13,13,14,14,14,14,15,15,15,15,15,15, - 16,16,16,16,16,16,17,17,18,18,18,18,18,18,19); + 409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509); my @_prime_next_small = ( 2,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23, 29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41,41,43,43,47, @@ -139,6 +141,19 @@ my @_prime_indices = (1, 7, 11, 13, 17, 19, 23, 29); my @_nextwheel30 = (1,7,7,7,7,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23,29,29,29,29,29,29,1); my @_prevwheel30 = (29,29,1,1,1,1,1,1,7,7,7,7,11,11,13,13,13,13,17,17,19,19,19,19,23,23,23,23,23,23); +sub _tiny_prime_count { + my($n) = @_; + return if $n >= $_primes_small[-1]; + my $j = $#_primes_small; + my $i = 1 + ($n >> 4); + while ($i < $j) { + my $mid = ($i+$j)>>1; + if ($_primes_small[$mid] <= $n) { $i = $mid+1; } + else { $j = $mid; } + } + return $i-1; +} + sub _is_prime7 { # n must not be divisible by 2, 3, or 5 my($n) = @_; @@ -966,6 +981,263 @@ sub nth_prime { $prime; } +# The nth prime will be less or equal to this number +sub nth_prime_upper { + my($n) = @_; + _validate_positive_integer($n); + + return $_primes_small[$n] if $n <= $#_primes_small; + + $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45; + + my $flogn = log($n); + my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n) + + my $upper; + if ($n >= 688383) { # Dusart 2010 page 2 + $upper = $n * ( $flogn + $flog2n - 1.0 + (($flog2n-2.00)/$flogn) ); + } elsif ($n >= 178974) { # Dusart 2010 page 7 + $upper = $n * ( $flogn + $flog2n - 1.0 + (($flog2n-1.95)/$flogn) ); + } elsif ($n >= 39017) { # Dusart 1999 page 14 + $upper = $n * ( $flogn + $flog2n - 0.9484 ); + } elsif ($n >= 6) { # Modified Robin 1983, for 6-39016 only + $upper = $n * ( $flogn + 0.6000 * $flog2n ); + } else { + $upper = $n * ( $flogn + $flog2n ); + } + + return int($upper + 1.0); +} + +# The nth prime will be greater than or equal to this number +sub nth_prime_lower { + my($n) = @_; + _validate_num($n) || _validate_positive_integer($n); + + return $_primes_small[$n] if $n <= $#_primes_small; + + $n = _upgrade_to_float($n) if $n > MPU_MAXPRIMEIDX || $n > 2**45; + + my $flogn = log($n); + my $flog2n = log($flogn); # Note distinction between log_2(n) and log^2(n) + + # Dusart 1999 page 14, for all n >= 2 + #my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.25)/$flogn)); + # Dusart 2010 page 2, for all n >= 3 + my $lower = $n * ($flogn + $flog2n - 1.0 + (($flog2n-2.10)/$flogn)); + + return int($lower); +} + +sub nth_prime_approx { + my($n) = @_; + _validate_num($n) || _validate_positive_integer($n); + + return $_primes_small[$n] if $n <= $#_primes_small; + + $n = _upgrade_to_float($n) + if ref($n) eq 'Math::BigInt' || $n >= MPU_MAXPRIMEIDX; + + my $flogn = log($n); + my $flog2n = log($flogn); + + # Cipolla 1902: + # m=0 fn * ( flogn + flog2n - 1 ); + # m=1 + ((flog2n - 2)/flogn) ); + # m=2 - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn)) + # + O((flog2n/flogn)^3) + # + # Shown in Dusart 1999 page 12, as well as other sources such as: + # http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf + # where the main issue you run into is that you're doing polynomial + # interpolation, so it oscillates like crazy with many high-order terms. + # Hence I'm leaving it at m=2. + + my $approx = $n * ( $flogn + $flog2n - 1 + + (($flog2n - 2)/$flogn) + - ((($flog2n*$flog2n) - 6*$flog2n + 11) / (2*$flogn*$flogn)) + ); + + # Apply a correction to help keep values close. + my $order = $flog2n/$flogn; + $order = $order*$order*$order * $n; + + if ($n < 259) { $approx += 10.4 * $order; } + elsif ($n < 775) { $approx += 7.52* $order; } + elsif ($n < 1271) { $approx += 5.6 * $order; } + elsif ($n < 2000) { $approx += 5.2 * $order; } + elsif ($n < 4000) { $approx += 4.3 * $order; } + elsif ($n < 12000) { $approx += 3.0 * $order; } + elsif ($n < 150000) { $approx += 2.1 * $order; } + elsif ($n < 200000000) { $approx += 0.0 * $order; } + else { $approx += -0.010 * $order; } + # $approx = -0.025 is better for the last, but it gives problems with some + # other code that always wants the asymptotic approximation to be >= actual. + + return int($approx + 0.5); +} + +############################################################################# + +sub prime_count_approx { + my($x) = @_; + _validate_num($x) || _validate_positive_integer($x); + + # Turn on high precision FP if they gave us a big number. + $x = _upgrade_to_float($x) if ref($_[0]) eq 'Math::BigInt'; + # Method 10^10 %error 10^19 %error + # ----------------- ------------ ------------ + # n/(log(n)-1) .22% .06% + # average bounds .01% .0002% + # li(n) .0007% .00000004% + # li(n)-li(n^.5)/2 .0004% .00000001% + # R(n) .0004% .00000001% + # + # Also consider: http://trac.sagemath.org/sage_trac/ticket/8135 + + # my $result = int( (prime_count_upper($x) + prime_count_lower($x)) / 2); + # my $result = int( LogarithmicIntegral($x) ); + # my $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2); + # my $result = RiemannR($x) + 0.5; + + # Sadly my Perl RiemannR function is really slow for big values. If MPFR + # is available, then use it -- it rocks. Otherwise, switch to LiCorr for + # very big values. This is hacky and shouldn't be necessary. + my $result; + if ( $x < 1e36 || _MPFR_available() ) { + 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); + } + $result = RiemannR($x) + 0.5; + } else { + $result = int(LogarithmicIntegral($x) - LogarithmicIntegral(sqrt($x))/2); + } + + return Math::BigInt->new($result->bfloor->bstr()) if ref($result) eq 'Math::BigFloat'; + return int($result); +} + +sub prime_count_lower { + my($x) = @_; + _validate_num($x) || _validate_positive_integer($x); + + return _tiny_prime_count($x) if $x < $_primes_small[-1]; + + $x = _upgrade_to_float($x) + if ref($x) eq 'Math::BigInt' || ref($_[0]) eq 'Math::BigInt'; + + my $flogx = log($x); + + # Chebyshev: 1*x/logx x >= 17 + # 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. + + my $result; + if ($x > 1000_000_000_000 && Math::Prime::Util::prime_get_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); +} + +sub prime_count_upper { + my($x) = @_; + _validate_num($x) || _validate_positive_integer($x); + + # Give an exact answer for what we have in our little table. + return _tiny_prime_count($x) if $x < $_primes_small[-1]; + + $x = _upgrade_to_float($x) + if ref($x) eq 'Math::BigInt' || ref($_[0]) eq 'Math::BigInt'; + + # Chebyshev: 1.25506*x/logx x >= 17 + # Rosser & Schoenfeld: x/(logx-3/2) x >= 67 + # Dusart 1999: x/logx*(1+1/logx+2.51/logxlogx) x >= 355991 + # Dusart 2010: x/logx*(1+1/logx+2.334/logxlogx) x >= 2_953_652_287 + + # As with the lower bounds, Dusart bounds are best by far. + + # Another possibility here for numbers under 3000M is to use Li(x) + # minus a correction. + + my $flogx = log($x); + + my $result; + if ($x > 10000_000_000_000 && Math::Prime::Util::prime_get_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. + $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); +} + + +############################################################################# + sub _mulmod { my($x, $y, $n) = @_; return (($x * $y) % $n) if ($x|$y) < MPU_HALFWORD; @@ -2553,6 +2825,41 @@ sub ecm_factor { @factors; } +sub divisors { + my($n) = @_; + _validate_num($n) || _validate_positive_integer($n); + + # In scalar context, returns sigma_0(n). Very fast. + return Math::Prime::Util::divisor_sum($n,0) unless wantarray; + return ($n == 0) ? (0,1) : (1) if $n <= 1; + + my %all_factors; + my @factors = Math::Prime::Util::factor($n); + return (1,$n) if scalar @factors == 1; + + if (ref($n) eq 'Math::BigInt') { + foreach my $f1 (@factors) { + my $big_f1 = Math::BigInt->new("$f1"); + my @to_add = map { ($_ <= ''.~0) ? _bigint_to_int($_) : $_ } + grep { $_ < $n } + map { $big_f1 * $_ } + keys %all_factors; + undef @all_factors{ $f1, @to_add }; + } + } else { + foreach my $f1 (@factors) { + my @to_add = grep { $_ < $n } + map { $f1 * $_ } + keys %all_factors; + undef @all_factors{ $f1, @to_add }; + } + } + # Add 1 and n + undef $all_factors{1}; + undef $all_factors{$n}; + my @divisors = sort {$a<=>$b} keys %all_factors; + return @divisors; +} sub ExponentialIntegral { diff --git a/lmo.c b/lmo.c index 3d40af5..3db1ad1 100644 --- a/lmo.c +++ b/lmo.c @@ -293,8 +293,8 @@ static UV _phi_recurse(UV x, UV a) { UV i, c = (a > PHIC) ? PHIC : a; UV sum = tablephi(x, c); if (a > c) { - UV p = _XS_nth_prime(c); - UV pa = _XS_nth_prime(a); + UV p = nth_prime(c); + UV pa = nth_prime(a); for (i = c+1; i <= a; i++) { UV xp; p = next_prime(p); @@ -344,7 +344,7 @@ UV legendre_phi(UV x, UV a) { uint32_t lastidx; UV res, max_cache_a = (a >= PHICACHEA) ? PHICACHEA : a+1; Newz(0, cache, PHICACHEX * max_cache_a, uint16_t); - primes = make_primelist(_XS_nth_prime(a+1), &lastidx); + primes = make_primelist(nth_prime(a+1), &lastidx); res = (UV) _phi(x, a, 1, primes, lastidx, cache); Safefree(primes); Safefree(cache); diff --git a/util.c b/util.c index 3a6dca8..e8e8e31 100644 --- a/util.c +++ b/util.c @@ -29,12 +29,14 @@ extern long double logl(long double); extern long double fabsl(long double); extern long double floorl(long double); + extern long double ceill(long double); #else #define powl(x, y) (long double) pow( (double) (x), (double) (y) ) #define expl(x) (long double) exp( (double) (x) ) #define logl(x) (long double) log( (double) (x) ) #define fabsl(x) (long double) fabs( (double) (x) ) #define floorl(x) (long double) floor( (double) (x) ) + #define ceill(x) (long double) ceil( (double) (x) ) #endif #ifdef LDBL_INFINITY @@ -570,7 +572,90 @@ UV _XS_prime_count(UV low, UV high) return count; } +UV prime_count_approx(UV n) +{ + if (n < 3000000) return _XS_prime_count(2, n); + return (UV) (_XS_RiemannR( (long double) n ) + 0.5 ); +} + +UV prime_count_lower(UV n) +{ + long double fn, flogn, lower, a; + + if (n < 33000) return _XS_prime_count(2, n); + + fn = (long double) n; + flogn = logl(n); + + if (n < 176000) a = 1.80; + else if (n < 315000) a = 2.10; + else if (n < 1100000) a = 2.20; + else if (n < 4500000) a = 2.31; + else if (n <233000000) a = 2.36; +#if BITS_PER_WORD == 32 + else a = 2.32; +#else + else if (n < UVCONST( 5433800000)) a = 2.32; + else if (n < UVCONST(60000000000)) a = 2.15; + else a = 2.00; +#endif + + lower = fn/flogn * (1.0 + 1.0/flogn + a/(flogn*flogn)); + return (UV) floorl(lower); +} + +typedef struct { + UV thresh; + float aval; +} thresh_t; + +static const thresh_t _upper_thresh[] = { + { 59000, 2.48 }, + { 350000, 2.52 }, + { 355991, 2.54 }, + { 356000, 2.51 }, + { 3550000, 2.50 }, + { 3560000, 2.49 }, + { 5000000, 2.48 }, + { 8000000, 2.47 }, + { 13000000, 2.46 }, + { 18000000, 2.45 }, + { 31000000, 2.44 }, + { 41000000, 2.43 }, + { 48000000, 2.42 }, + { 119000000, 2.41 }, + { 182000000, 2.40 }, + { 192000000, 2.395 }, + { 213000000, 2.390 }, + { 271000000, 2.385 }, + { 322000000, 2.380 }, + { 400000000, 2.375 }, + { 510000000, 2.370 }, + { 682000000, 2.367 }, + { UVCONST(2953652287), 2.362 } +}; +#define NUPPER_THRESH (sizeof(_upper_thresh)/sizeof(_upper_thresh[0])) + +UV prime_count_upper(UV n) +{ + int i; + long double fn, flogn, upper, a; + + if (n < 33000) return _XS_prime_count(2, n); + + fn = (long double) n; + flogn = logl(n); + + for (i = 0; i < NUPPER_THRESH; i++) + if (n < _upper_thresh[i].thresh) + break; + + if (i < NUPPER_THRESH) a = _upper_thresh[i].aval; + else a = 2.334; /* Dusart 2010, page 2 */ + upper = fn/flogn * (1.0 + 1.0/flogn + a/(flogn*flogn)); + return (UV) ceill(upper); +} static const unsigned short primes_small[] = {0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97, @@ -580,18 +665,17 @@ static const unsigned short primes_small[] = 409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499}; #define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0])) -/* Note: We're keeping this here because we use it for nth_prime */ /* The nth prime will be less or equal to this number */ -static UV _XS_nth_prime_upper(UV n) +UV nth_prime_upper(UV n) { - double fn, flogn, flog2n, upper; + long double fn, flogn, flog2n, upper; if (n < NPRIMES_SMALL) return primes_small[n]; - fn = (double) n; - flogn = log(n); - flog2n = log(flogn); /* Note distinction between log_2(n) and log^2(n) */ + fn = (long double) n; + flogn = logl(n); + flog2n = logl(flogn); /* Note distinction between log_2(n) and log^2(n) */ if (n >= 688383) /* Dusart 2010 page 2 */ upper = fn * (flogn + flog2n - 1.0 + ((flog2n-2.00)/flogn)); @@ -615,16 +699,73 @@ static UV _XS_nth_prime_upper(UV n) * nth_prime_lower(n) <= nth_prime(n) <= nth_prime_upper(n) */ /* Watch out for overflow */ - if (upper >= (double)UV_MAX) { + if (upper >= (long double)UV_MAX) { if (n <= MPU_MAX_PRIME_IDX) return MPU_MAX_PRIME; croak("nth_prime_upper(%"UVuf") overflow", n); } - return (UV) ceil(upper); + return (UV) ceill(upper); +} + +/* The nth prime will be greater than or equal to this number */ +UV nth_prime_lower(UV n) +{ + long double fn, flogn, flog2n, lower; + + if (n < NPRIMES_SMALL) + return primes_small[n]; + + fn = (long double) n; + flogn = logl(n); + flog2n = logl(flogn); /* Note distinction between log_2(n) and log^2(n) */ + + /* Dusart 2010 page 2, for all n >= 3 */ + lower = fn * (flogn + flog2n - 1.0 + ((flog2n-2.10)/flogn)); + + return (UV) floorl(lower); +} + +UV nth_prime_approx(UV n) +{ + long double fn, flogn, flog2n, approx, order; + + if (n < NPRIMES_SMALL) + return primes_small[n]; + + fn = (long double) n; + flogn = logl(n); + flog2n = logl(flogn); /* Note distinction between log_2(n) and log^2(n) */ + + /* Cipolla 1902: + * m=0 fn * ( flogn + flog2n - 1 ); + * m=1 + ((flog2n - 2)/flogn) ); + * m=2 - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn)) + * + O((flog2n/flogn)^3) + */ + + approx = fn * ( flogn + flog2n - 1.0 + + ((flog2n - 2.0) / flogn) + - (((flog2n*flog2n) - 6.0*flog2n + 11.0) / (2*flogn*flogn)) + ); + + /* Apply a correction */ + order = flog2n / flogn; + order = order * order * order * fn; + if (n < 259) { approx += 10.4 * order; } + else if (n < 775) { approx += 7.52 * order; } + else if (n < 1271) { approx += 5.6 * order; } + else if (n < 2000) { approx += 5.2 * order; } + else if (n < 4000) { approx += 4.3 * order; } + else if (n < 12000) { approx += 3.0 * order; } + else if (n < 150000) { approx += 2.1 * order; } + else if (n <200000000) { } + else { approx += -0.01 * order; } /* -0.25 is closer */ + + return (UV) floorl(approx + 0.5); } -UV _XS_nth_prime(UV n) +UV nth_prime(UV n) { const unsigned char* cache_sieve; unsigned char* segment; @@ -638,7 +779,7 @@ UV _XS_nth_prime(UV n) return primes_small[n]; /* Determine a bound on the nth prime. We know it comes before this. */ - upper_limit = _XS_nth_prime_upper(n); + upper_limit = nth_prime_upper(n); MPUassert(upper_limit > 0, "nth_prime got an upper limit of 0"); /* For relatively small values, generate a sieve and count the results. diff --git a/util.h b/util.h index 715e4b7..76b9fb5 100644 --- a/util.h +++ b/util.h @@ -13,7 +13,13 @@ extern UV next_prime(UV x); extern UV prev_prime(UV x); extern UV _XS_prime_count(UV low, UV high); -extern UV _XS_nth_prime(UV x); +extern UV nth_prime(UV x); +extern UV nth_prime_upper(UV x); +extern UV nth_prime_lower(UV x); +extern UV nth_prime_approx(UV x); +extern UV prime_count_lower(UV x); +extern UV prime_count_upper(UV x); +extern UV prime_count_approx(UV x); extern signed char* _moebius_range(UV low, UV high); extern UV* _totient_range(UV low, UV high); -- 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