[libmath-prime-util-perl] 02/07: Pure Perl working on all tests
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:41 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.09
in repository libmath-prime-util-perl.
commit cc71cd27829abb53b6523bd3b7811858ebe9f796
Author: Dana Jacobsen <dana at acm.org>
Date: Sat Jun 23 11:47:20 2012 -0600
Pure Perl working on all tests
---
Changes | 2 +-
lib/Math/Prime/Util.pm | 27 +++--
lib/Math/Prime/Util/PP.pm | 267 +++++++++++++++++++++++++++++++++++++++++++---
util.c | 2 +-
4 files changed, 273 insertions(+), 25 deletions(-)
diff --git a/Changes b/Changes
index 99cde38..0a6c18f 100644
--- a/Changes
+++ b/Changes
@@ -1,7 +1,7 @@
Revision history for Perl extension Math::Prime::Util.
0.09 23 June 2012
- - Started work on pure perl code.
+ - Pure Perl code. Passes all tests, but 1 to 13000x slower.
0.08 22 June 2012
- Added thread safety and tested good concurrency.
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index d2d5e83..565396b 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -69,7 +69,9 @@ BEGIN {
*prho_factor = \&Math::Prime::Util::PP::prho_factor;
*pminus1_factor = \&Math::Prime::Util::PP::pminus1_factor;
- # TODO: Ei(x), li(x), R(x)
+ *RiemannR = \&Math::Prime::Util::PP::RiemannR;
+ *LogarithmicIntegral = \&Math::Prime::Util::PP::LogarithmicIntegral;
+ *ExponentialIntegral = \&Math::Prime::Util::PP::ExponentialIntegral;
# TODO: prime_count is horribly, horribly slow
# TODO: We should have some tests to verify XS vs. PP.
}
@@ -85,7 +87,7 @@ my $_maxprimeidx=(_maxbits == 32) ? 203280221 : 425656284035217743;
sub primes {
- my $optref = {}; $optref = shift if ref $_[0] eq 'HASH';
+ my $optref = (ref $_[0] eq 'HASH') ? shift : {};
croak "no parameters to primes" unless scalar @_ > 0;
croak "too many parameters to primes" unless scalar @_ <= 2;
my $low = (@_ == 2) ? shift : 2;
@@ -96,8 +98,8 @@ sub primes {
if ( (!defined $low) || (!defined $high) ||
($low =~ tr/0123456789//c) || ($high =~ tr/0123456789//c)
) {
- $low = 'undef' unless defined $low;
- $high = 'undef' unless defined $high;
+ $low = '<undef>' unless defined $low;
+ $high = '<undef>' unless defined $high;
croak "Parameters [ $low $high ] must be positive integers";
}
@@ -155,6 +157,7 @@ sub primes {
return $sref;
}
+
sub random_prime {
my $low = (@_ == 2) ? shift : 2;
my $high = shift;
@@ -170,8 +173,8 @@ sub random_prime {
# Make sure we have a valid range.
# TODO: this is is killing performance with large numbers
- $low = next_prime($low-1);
- $high = ($high < ~0) ? prev_prime($high+1) : prev_prime($high);
+ $low = next_prime($low - 1);
+ $high = ($high < ~0) ? prev_prime($high + 1) : prev_prime($high);
return $low if ($low == $high) && is_prime($low);
return if $low >= $high;
@@ -198,18 +201,19 @@ sub random_prime {
do {
$prime = $low + $irandf->($range);
croak "Random function broken?" if $loop_limit-- < 0;
- } while ( !($prime%2) || !($prime%3) || !is_prime($prime) );
+ } while ( !($prime % 2 ) || !($prime % 3) || !is_prime($prime) );
} else {
do {
my $rand = ( ($irandf->(4294967295) << 32) + $irandf->(4294967295) ) % $range;
$prime = $low + $rand;
croak "Random function broken?" if $loop_limit-- < 0;
- } while ( !($prime%2) || !($prime%3) || !is_prime($prime) );
+ } while ( !($prime % 2) || !($prime % 3) || !is_prime($prime) );
}
}
return $prime;
}
+
sub random_ndigit_prime {
my $digits = shift;
if ((!defined $digits) || ($digits > $_maxdigits) || ($digits < 1)) {
@@ -223,12 +227,14 @@ sub random_ndigit_prime {
# Perhaps a random_nbit_prime ? Definition?
+
sub all_factors {
my $n = shift;
my @factors = factor($n);
my %all_factors;
foreach my $f1 (@factors) {
next if $f1 >= $n;
+ # We're adding to %all_factors in the loop, so grab the keys now.
my @all = keys %all_factors;;
foreach my $f2 (@all) {
$all_factors{$f1*$f2} = 1 if ($f1*$f2) < $n;
@@ -835,21 +841,26 @@ C<is_prime>: my impressions:
Math::Primality Very slow Very slow
The differences are in the implementations:
+
- L<Math::Prime::FastSieve> only works in a sieved range, which is really
fast if you can do it (M::P::U will do the same if you call
C<prime_precalc>). Larger inputs just need too much time and memory
for the sieve.
+
- L<Math::Primality> uses a GMP BPSW test which is overkill for our 64-bit
range. It's generally an order of magnitude or two slower than any
of the others.
+
- L<Math::Pari> has some very effective code, but it has some overhead to get
to it from Perl. That means for small numbers it is relatively slow: an
order of magnitude slower than M::P::XS and M::P::Util (though arguably
this is only important for benchmarking since "slow" is ~2 microseconds).
Large numbers transition over to smarter tests so don't slow down much.
+
- L<Math::Prime::XS> does trial divisions, which is wonderful if the input
has a small factor (or is small itself). But it can take 1000x longer
if given a large prime.
+
- L<Math::Prime::Util> looks in the sieve for a fast bit lookup if that
exists (default up to 30,000 but it can be expanded, e.g.
C<prime_precalc>), uses trial division for numbers higher than this but
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 8e27c3f..93b604d 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -99,6 +99,44 @@ sub is_prime {
2*_is_prime7($n);
}
+sub _sieve_erat {
+ my($end) = @_;
+ my $last = int(($end+1)/2);
+
+ my $sieve = '';
+ my $n = 3;
+ while ( ($n*$n) <= $end ) {
+ my $s = $n*$n;
+ while ($s <= $end) {
+ vec($sieve, $s >> 1, 1) = 1;
+ $s += 2*$n;
+ }
+ do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0;
+ }
+ # Mark 1 as composite
+ vec($sieve, 0, 1) = 1;
+ $sieve;
+}
+# Uses 8x more memory, but about 50% faster
+sub _sieve_erat_array {
+ my($end) = @_;
+ my $last = int(($end+1)/2);
+
+ my @sieve;
+ my $n = 3;
+ while ( ($n*$n) <= $end ) {
+ my $s = $n*$n;
+ while ($s <= $end) {
+ $sieve[$s>>1] = 1;
+ $s += 2*$n;
+ }
+ do { $n += 2 } while $sieve[$n>>1];
+ }
+ # Mark 1 as composite
+ $sieve[0] = 1;
+ \@sieve;
+}
+
sub primes {
my $optref = (ref $_[0] eq 'HASH') ? shift : {};
croak "no parameters to primes" unless scalar @_ > 0;
@@ -119,14 +157,23 @@ sub primes {
push @$sref, 5 if ($low <= 5) && ($high >= 5);
$low = 7 if $low < 7;
- my $n = $low;
- my $base = 30 * int($n/30);
- my $in = 0; $in++ while ($n - $base) > $_prime_indices[$in];
- $n = $base + $_prime_indices[$in];
- while ($n <= $high) {
- push @$sref, $n if _is_prime7($n);
- if (++$in == 8) { $base += 30; $in = 0; }
+ if ($low == 7) {
+ my $sieve = _sieve_erat_array($high);
+ my $n = 7;
+ while ($n <= $high) {
+ push @$sref, $n if !$sieve->[$n>>1];
+ $n += 2;
+ }
+ } else {
+ my $n = $low;
+ my $base = 30 * int($n/30);
+ my $in = 0; $in++ while ($n - $base) > $_prime_indices[$in];
$n = $base + $_prime_indices[$in];
+ while ($n <= $high) {
+ push @$sref, $n if _is_prime7($n);
+ if (++$in == 8) { $base += 30; $in = 0; }
+ $n = $base + $_prime_indices[$in];
+ }
}
$sref;
}
@@ -203,14 +250,27 @@ sub prime_count {
$count++ if ($low <= 5) && ($high >= 5);
$low = 7 if $low < 7;
- my $n = $low;
- my $base = 30 * int($n/30);
- my $in = 0; $in++ while ($n - $base) > $_prime_indices[$in];
- $n = $base + $_prime_indices[$in];
- while ($n <= $high) {
- $count++ if _is_prime7($n);
- if (++$in == 8) { $base += 30; $in = 0; }
+ if ($low == 7) {
+ my $sieve = _sieve_erat_array($high);
+ my $n = 7;
+ while ($n <= $high) {
+ $count++ if !$sieve->[$n>>1];
+ $n += 4;
+ if ($n <= $high) {
+ $count++ if !$sieve->[$n>>1];
+ $n += 2;
+ }
+ }
+ } else {
+ my $n = $low;
+ my $base = 30 * int($n/30);
+ my $in = 0; $in++ while ($n - $base) > $_prime_indices[$in];
$n = $base + $_prime_indices[$in];
+ while ($n <= $high) {
+ $count++ if _is_prime7($n);
+ if (++$in == 8) { $base += 30; $in = 0; }
+ $n = $base + $_prime_indices[$in];
+ }
}
$count;
}
@@ -290,7 +350,8 @@ sub prime_count_approx {
return $_prime_count_small[$x] if $x <= $#_prime_count_small;
- return int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
+ #return int( (prime_count_upper($x) + prime_count_lower($x)) / 2);
+ return int(RiemannR($x)+0.5);
}
sub nth_prime_lower {
@@ -542,6 +603,182 @@ sub pbrent_factor { trial_factor(@_) }
sub prho_factor { trial_factor(@_) }
sub pminus1_factor { trial_factor(@_) }
+
+my $_const_euler = 0.57721566490153286060651209008240243104215933593992;
+my $_const_li2 = 1.045163780117492784844588889194613136522615578151;
+
+sub ExponentialIntegral {
+ my($x) = @_;
+ my $tol = 1e-16;
+ my $sum = 0.0;
+ my($y, $t);
+ my $c = 0.0;
+
+ croak "Invalid input to ExponentialIntegral: x must be != 0" if $x == 0;
+
+ my $val; # The result from one of the four methods
+
+ if ($x < -1) {
+ # Continued fraction
+ my $lc = 0;
+ my $ld = 1 / (1 - $x);
+ $val = $ld * (-exp($x));
+ for my $n (1 .. 100000) {
+ $lc = 1 / (2*$n + 1 - $x - $n*$n*$lc);
+ $ld = 1 / (2*$n + 1 - $x - $n*$n*$ld);
+ my $old = $val;
+ $val *= $ld/$lc;
+ last if abs($val - $old) <= ($tol * abs($val));
+ }
+ } elsif ($x < 0) {
+ # Rational Chebyshev approximation
+ my @C6p = ( -148151.02102575750838086,
+ 150260.59476436982420737,
+ 89904.972007457256553251,
+ 15924.175980637303639884,
+ 2150.0672908092918123209,
+ 116.69552669734461083368,
+ 5.0196785185439843791020);
+ my @C6q = ( 256664.93484897117319268,
+ 184340.70063353677359298,
+ 52440.529172056355429883,
+ 8125.8035174768735759866,
+ 750.43163907103936624165,
+ 40.205465640027706061433,
+ 1.0000000000000000000000);
+ my $sumn = $C6p[0]-$x*($C6p[1]-$x*($C6p[2]-$x*($C6p[3]-$x*($C6p[4]-$x*($C6p[5]-$x*$C6p[6])))));
+ my $sumd = $C6q[0]-$x*($C6q[1]-$x*($C6q[2]-$x*($C6q[3]-$x*($C6q[4]-$x*($C6q[5]-$x*$C6q[6])))));
+ $val = log(-$x) - ($sumn / $sumd);
+ } elsif ($x < -log($tol)) {
+ # Convergent series
+ my $fact_n = 1;
+ $y = $_const_euler-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ $y = log($x)-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ for my $n (1 .. 200) {
+ $fact_n *= $x/$n;
+ my $term = $fact_n / $n;
+ $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ last if $term < $tol;
+ }
+ $val = $sum;
+ } else {
+ # Asymptotic divergent series
+ $val = exp($x) / $x;
+ $y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ my $term = 1.0;
+ for my $n (1 .. 200) {
+ my $last_term = $term;
+ $term *= $n/$x;
+ last if $term < $tol;
+ if ($term < $last_term) {
+ $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ } else {
+ $y = (-$last_term/3)-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ last;
+ }
+ }
+ $val *= $sum;
+ }
+ $val;
+}
+
+sub LogarithmicIntegral {
+ my($x) = @_;
+ return 0 if $x == 0;
+ return 0+(-Infinity) if $x == 1;
+ return $_const_li2 if $x == 2;
+ croak "Invalid input to LogarithmicIntegral: x must be > 0" if $x <= 0;
+ ExponentialIntegral(log($x));
+}
+
+# Riemann Zeta function for integers, used for computing Riemann R
+my @_Riemann_Zeta_Table = (
+ 0.6449340668482264364724151666460251892, # zeta(2)
+ 0.2020569031595942853997381615114499908,
+ 0.0823232337111381915160036965411679028,
+ 0.0369277551433699263313654864570341681,
+ 0.0173430619844491397145179297909205279,
+ 0.0083492773819228268397975498497967596,
+ 0.0040773561979443393786852385086524653,
+ 0.0020083928260822144178527692324120605,
+ 0.0009945751278180853371459589003190170,
+ 0.0004941886041194645587022825264699365,
+ 0.0002460865533080482986379980477396710,
+ 0.0001227133475784891467518365263573957,
+ 0.0000612481350587048292585451051353337,
+ 0.0000305882363070204935517285106450626,
+ 0.0000152822594086518717325714876367220,
+ 0.0000076371976378997622736002935630292,
+ 0.0000038172932649998398564616446219397,
+ 0.0000019082127165539389256569577951013,
+ 0.0000009539620338727961131520386834493,
+ 0.0000004769329867878064631167196043730,
+ 0.0000002384505027277329900036481867530,
+ 0.0000001192199259653110730677887188823,
+ 0.0000000596081890512594796124402079358,
+ 0.0000000298035035146522801860637050694,
+ 0.0000000149015548283650412346585066307,
+ 0.0000000074507117898354294919810041706,
+ 0.0000000037253340247884570548192040184,
+ 0.0000000018626597235130490064039099454,
+ 0.0000000009313274324196681828717647350,
+ 0.0000000004656629065033784072989233251,
+ 0.0000000002328311833676505492001455976,
+ 0.0000000001164155017270051977592973835,
+ 0.0000000000582077208790270088924368599,
+ 0.0000000000291038504449709968692942523,
+ 0.0000000000145519218910419842359296322,
+ 0.0000000000072759598350574810145208690,
+ 0.0000000000036379795473786511902372363,
+ 0.0000000000018189896503070659475848321,
+ 0.0000000000009094947840263889282533118,
+);
+
+# Compute Riemann Zeta function. Slow and inaccurate near x = 1, but improves
+# very rapidly (x = 5 is quite reasonable).
+sub _evaluate_zeta {
+ my($x) = @_;
+ my $tol = 1e-16;
+ my $sum = 0.0;
+ my($y, $t);
+ my $c = 0.0;
+
+ $y = (1.0/(2.0**$x))-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+
+ for my $k (3 .. 100000) {
+ my $term = 1.0 / ($k ** $x);
+ $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ last if abs($term) < $tol;
+ }
+ $sum;
+}
+
+# Riemann R function
+sub RiemannR {
+ my($x) = @_;
+ my $tol = 1e-16;
+ my $sum = 0.0;
+ my($y, $t);
+ my $c = 0.0;
+
+ croak "Invalid input to ReimannR: x must be > 0" if $x <= 0;
+
+ $y = 1.0-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ my $flogx = log($x);
+ my $part_term = 1.0;
+
+ for my $k (1 .. 10000) {
+ # Small k from table, larger k from function
+ my $zeta = ($k <= $#_Riemann_Zeta_Table) ? $_Riemann_Zeta_Table[$k+1-2]
+ : _evaluate_zeta($k+1);
+ $part_term *= $flogx / $k;
+ my $term = $part_term / ($k + $k * $zeta);
+ $y = $term-$c; $t = $sum+$y; $c = ($t-$sum)-$y; $sum = $t;
+ last if abs($term) < $tol;
+ }
+ $sum;
+}
+
1;
__END__
diff --git a/util.c b/util.c
index a318f53..d3d040a 100644
--- a/util.c
+++ b/util.c
@@ -919,7 +919,7 @@ double LogarithmicIntegral(double x) {
if (x == 0) return 0;
if (x == 1) return -INFINITY;
if (x == 2) return li2;
- if (x <= 0) croak("Invalid input to ExponentialIntegral: x must be > 0");
+ if (x <= 0) croak("Invalid input to LogarithmicIntegral: x must be > 0");
return ExponentialIntegral(log(x));
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list