[libmath-prime-util-perl] 39/72: Minor cleanup
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:49:39 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.32
in repository libmath-prime-util-perl.
commit 7be347d17aefbdffbdfe1b79673c0af7eef6fa8f
Author: Dana Jacobsen <dana at acm.org>
Date: Fri Sep 20 14:37:31 2013 -0700
Minor cleanup
---
TODO | 10 ++++++----
XS.xs | 14 +++++++-------
lib/Math/Prime/Util.pm | 44 ++++++++++++++++++++++----------------------
sieve.c | 7 +++++--
sieve.h | 4 ++--
5 files changed, 42 insertions(+), 37 deletions(-)
diff --git a/TODO b/TODO
index e9750a3..3ca6697 100644
--- a/TODO
+++ b/TODO
@@ -1,7 +1,4 @@
-- Examine behavior near 32-bit limit on 32-bit machines.
- (done for factoring)
-
- segment sieve should itself use a segment for its primes.
Today we'd need sqrt(2^64) max = 140MB. Segmenting would yield under 1MB.
@@ -43,4 +40,9 @@
- Use Montgomery routines in more places: Factoring.
-- Tests for PrimeIterator
+- More tests for PrimeIterator.
+
+- Put euler_phi and moebius directly in XS.
+ (1) optional second argument. Easily handled, and not hard to do in
+ generic sub call.
+ (2) generic sub returns an array. This is the sticking point.
diff --git a/XS.xs b/XS.xs
index ea74cb6..daf3125 100644
--- a/XS.xs
+++ b/XS.xs
@@ -665,21 +665,21 @@ _XS_moebius(IN UV lo, IN UV hi = 0)
Safefree(mu);
} else {
UV factors[MPU_MAX_FACTORS+1];
- UV nfactors, lastf;
+ UV nfactors;
UV n = lo;
if (n <= 1)
XSRETURN_IV(n);
- if ( (n >= 25) && ( !(n%4) || !(n%9) || !(n%25) ) )
+
+ if ( (!(n% 4) && n >= 4) || (!(n% 9) && n >= 9) ||
+ (!(n%25) && n >= 25) || (!(n%49) && n >= 49) )
XSRETURN_IV(0);
nfactors = factor(n, factors);
- lastf = 0;
- for (i = 0; i < nfactors; i++) {
- if (factors[i] == lastf)
+ for (i = 1; i < nfactors; i++)
+ if (factors[i] == factors[i-1])
XSRETURN_IV(0);
- lastf = factors[i];
- }
+
XSRETURN_IV( (nfactors % 2) ? -1 : 1 );
}
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 54076cf..6a0b325 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1254,7 +1254,8 @@ sub moebius {
my @mu = map { 1 } $lo .. $hi;
$mu[0] = 0 if $lo == 0;
my $sqrtn = int(sqrt($hi)+0.5);
- foreach my $p ( @{ primes($sqrtn) } ) {
+ forprimes {
+ my $p = $_;
my $i = $p * $p;
$i = $i * int($lo/$i) + (($lo % $i) ? $i : 0) if $i < $lo;
while ($i <= $hi) {
@@ -1267,7 +1268,7 @@ sub moebius {
$mu[$i-$lo] *= -$p;
$i += $p;
}
- }
+ } $sqrtn;
foreach my $i ($lo .. $hi) {
my $m = $mu[$i-$lo];
$m *= -1 if abs($m) != $i;
@@ -1280,10 +1281,8 @@ sub moebius {
# Quick check for small replicated factors
return 0 if ($n >= 25) && (!($n % 4) || !($n % 9) || !($n % 25));
my @factors = factor($n);
- my $lastf = 0;
- foreach my $factor (@factors) {
- return 0 if $factor == $lastf;
- $lastf = $factor;
+ foreach my $i (1 .. $#factors) {
+ return 0 if $factors[$i] == $factors[$i-1];
}
return (((scalar @factors) % 2) == 0) ? 1 : -1;
}
@@ -4245,9 +4244,6 @@ Here's the best way for PE187. Under 30 milliseconds from the command line:
=head1 LIMITATIONS
-I have not completed testing all the functions near the word size limit
-(e.g. C<2^32> for 32-bit machines). Please report any problems you find.
-
Perl versions earlier than 5.8.0 have a rather broken 64-bit implementation,
in that the values are accessed as doubles. Hence any value larger
than C<~ 2^49> will start losing bottom bits. This causes numerous functions
@@ -4283,7 +4279,8 @@ functionality with small integers. It's fast and simple, and has a good
set of features.
=item * L<Math::Primality> is the alternative module I use for primality
-testing on bigints.
+testing on bigints. The downside is that it can be slow, and the functions
+other than primality tests are I<very> slow.
=item * L<Math::Pari> if you want the kitchen sink and can install it and
handle using it. There are still some functions it doesn't do well
@@ -4292,15 +4289,17 @@ handle using it. There are still some functions it doesn't do well
=back
-L<Math::Prime::XS> has C<is_prime> and C<primes> functionality. There is no
-bigint support. The C<is_prime> function uses well-written trial division,
-meaning it is very fast for small numbers, but terribly slow for large
-64-bit numbers. Because MPU does input validation and bigint conversion,
-there is about 20 microseconds of additional overhead making MPXS a little
-faster for tiny inputs, but once over 700k, MPU is faster for all values.
-MPXS's prime sieve is an unoptimized non-segmented SoE which returns an
-array. It works well for 32-bit values, but speed and memory are problematic
-for larger values.
+L<Math::Prime::XS> has C<is_prime> and C<primes> functionality. There is
+no bigint support. The C<is_prime> function uses well-written trial
+division, meaning it is very fast for small numbers, but terribly slow for
+large 64-bit numbers. Because MPU does input validation and bigint
+conversion, there is about 20 microseconds of additional overhead making
+MPXS a little faster for tiny inputs, but once over 700k MPU is faster
+for all values. MPXS's prime sieve is an unoptimized non-segmented SoE
+which returns an array. Sieve bases larger than C<10^7> start taking
+inordinately long and using a lot of memory (gigabytes beyond C<10^10>).
+E.g. C<primes(10**9, 10**9+1000)> takes 36 seconds with MPXS, but only
+0.00015 seconds with MPU.
L<Math::Prime::FastSieve> supports C<primes>, C<is_prime>, C<next_prime>,
C<prev_prime>, C<prime_count>, and C<nth_prime>. The caveat is that all
@@ -4308,7 +4307,8 @@ functions only work within the sieved range, so are limited to about C<10^10>.
It uses a fast SoE to generate the main sieve. The sieve is 2-3x slower than
the base sieve for MPU, and is non-segmented so cannot be used for
larger values. Since the functions work with the sieve, they are very fast.
-All this functionality is present in MPU as well, though not required.
+The fast bit-vector-lookup functionality can be replicated in MPU using
+C<prime_precalc> but is not required.
L<Bit::Vector> supports the C<primes> and C<prime_count> functionality in a
somewhat similar way to L<Math::Prime::FastSieve>. It is the slowest of all
@@ -4335,7 +4335,7 @@ not support bigints. Both are implemented with trial division, meaning they
are very fast for really small values, but quickly become unusably slow
(factoring 19 digit semiprimes is over 700 times slower). It has additional
functions C<count_prime_factors> (possible in MPU using C<scalar factor($n)>)
-and C<matches> which has no direct equivalent.
+and C<matches> which has no equivalent.
L<Math::Big> version 1.12 includes C<primes> functionality. The current
code is only usable for very tiny inputs as it is incredibly slow and uses
@@ -4549,7 +4549,7 @@ C<is_prime>: my impressions for various sized inputs:
(2) Trial division only. Very fast if every factor is tiny.
(3) Too much memory to hold the sieve (11dig = 6GB, 12dig = ~50GB)
(4) By default L<Math::Pari> installs Pari 2.1.7, which uses 10 M-R tests
- for is_prime and is not fast. See notes below for 2.3.5.
+ for isprime and is not fast. See notes below for 2.3.5.
The differences are in the implementations:
diff --git a/sieve.c b/sieve.c
index 89caf16..95c90b9 100644
--- a/sieve.c
+++ b/sieve.c
@@ -233,8 +233,11 @@ int sieve_segment(unsigned char* mem, UV startd, UV endd)
START_DO_FOR_EACH_SIEVE_PRIME(sieve, 17, pcsize)
{
- /* p increments from 17 to at least sqrt(endp) */
- UV p2 = p*p; /* TODO: overflow */
+ /* p increments from 17 to at least sqrt(endp). Note on overflow:
+ * 32-bit: limit= 65535, max p = 65521, p*p = ~0-1965854
+ * 64-bit: limit=4294967295, max p = 4294967291, p*p = ~0-42949672934
+ * No overflow here, but possible after the incrementing below. */
+ UV p2 = p*p;
if (p2 > endp) break;
/* Find first multiple of p greater than p*p and larger than startp */
if (p2 < startp) {
diff --git a/sieve.h b/sieve.h
index ecf0093..7fd7610 100644
--- a/sieve.h
+++ b/sieve.h
@@ -132,8 +132,8 @@ static UV prev_prime_in_sieve(const unsigned char* sieve, UV p) {
} \
while ( p <= l_ ) {
-#define RETURN_FROM_EACH_PRIME(x) \
- do { release_prime_cache(sieve_); return(x); } while (0)
+#define RETURN_FROM_EACH_PRIME(retstmt) \
+ do { release_prime_cache(sieve_); retstmt; } while (0)
#define END_DO_FOR_EACH_PRIME \
if (p < 7) { \
--
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