Skip to content

Commit

Permalink
Don't compile Legendre/Meissel/Lehmer/LMOS prime counts by default
Browse files Browse the repository at this point in the history
  • Loading branch information
danaj committed Dec 31, 2013
1 parent df94359 commit 9d3bad1
Show file tree
Hide file tree
Showing 10 changed files with 59 additions and 52 deletions.
3 changes: 3 additions & 0 deletions Changes
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,9 @@ Revision history for Perl module Math::Prime::Util

- znorder uses Carmichael Lambda instead of Euler Phi. Faster.

- Legendre, Meissel, Lehmer, and LMOS methods are no longer compiled by
default. Add -DLEHMER if you really want access to them.


0.35 2013-12-08

Expand Down
6 changes: 3 additions & 3 deletions XS.xs
Original file line number Diff line number Diff line change
Expand Up @@ -231,9 +231,6 @@ _XS_nth_prime(IN UV n)
OUTPUT:
RETVAL

UV
_XS_legendre_phi(IN UV x, IN UV a)

SV*
sieve_primes(IN UV low, IN UV high)
ALIAS:
Expand Down Expand Up @@ -565,6 +562,9 @@ kronecker(IN SV* sva, IN SV* svb)
_vcallsubn(aTHX_ G_SCALAR, "_generic_kronecker", 2);
return; /* skip implicit PUTBACK */

UV
legendre_phi(IN UV n, IN UV a)

double
_XS_ExponentialIntegral(IN SV* x)
ALIAS:
Expand Down
2 changes: 1 addition & 1 deletion factor.c
Original file line number Diff line number Diff line change
Expand Up @@ -755,7 +755,7 @@ int pplus1_factor(UV n, UV *factors, UV B1)
{
UV X1, X2, f;
UV sqrtB1 = isqrt(B1);
MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in pminus1_factor");
MPUassert( (n >= 3) && ((n%2) != 0) , "bad n in pplus1_factor");

X1 = 7 % n;
X2 = 11 % n;
Expand Down
49 changes: 16 additions & 33 deletions lehmer.c
Original file line number Diff line number Diff line change
@@ -1,12 +1,10 @@
#if defined(LEHMER) || defined(PRIMESIEVE_STANDALONE)

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/* Below this size, just sieve (with table speedup). */
#define SIEVE_LIMIT 60000000
#define MAX_PHI_MEM (896*1024*1024)

/*****************************************************************************
*
* Lehmer prime counting utility. Calculates pi(x), count of primes <= x.
Expand Down Expand Up @@ -58,6 +56,10 @@
* Factorization", 2nd edition, 1994.
*/

/* Below this size, just sieve (with table speedup). */
#define SIEVE_LIMIT 60000000
#define MAX_PHI_MEM (896*1024*1024)

static int const verbose = 0;
#define STAGE_TIMING 0

Expand Down Expand Up @@ -844,35 +846,6 @@ UV _XS_LMOS_pi(UV n)
return sum;
}

static const unsigned char 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};
#define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))

UV _XS_legendre_phi(UV x, UV a) {
phitableinit();
/* For small values, calculate directly */
if (a <= PHIC) return tablephi(x, a);
/* For large values, do our non-recursive phi */
if (a >= NPRIMES_SMALL) return phi(x,a);
/* Otherwise, recurse */
{
UV i, sum = tablephi(x, PHIC);
for (i = PHIC+1; i <= a; i++) {
uint32_t p = primes_small[i];
UV xp = x/p;
if (xp < p) {
while (x < primes_small[a])
a--;
return (sum - a + i - 1);
}
sum -= _XS_legendre_phi(xp, i-1);
}
return sum;
}
}


#ifdef PRIMESIEVE_STANDALONE
int main(int argc, char *argv[])
{
Expand Down Expand Up @@ -907,3 +880,13 @@ int main(int argc, char *argv[])
return(0);
}
#endif

#else

#include "lehmer.h"
UV _XS_LMOS_pi(UV n) { croak("Not compiled with Lehmer support"); }
UV _XS_lehmer_pi(UV n) { croak("Not compiled with Lehmer support"); }
UV _XS_meissel_pi(UV n) { croak("Not compiled with Lehmer support"); }
UV _XS_legendre_pi(UV n) { croak("Not compiled with Lehmer support"); }

#endif
2 changes: 0 additions & 2 deletions lehmer.h
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,4 @@ extern UV _XS_meissel_pi(UV n);
extern UV _XS_lehmer_pi(UV n);
extern UV _XS_LMOS_pi(UV n);

extern UV _XS_legendre_phi(UV x, UV a);

#endif
24 changes: 24 additions & 0 deletions lmo.c
Original file line number Diff line number Diff line change
Expand Up @@ -259,6 +259,30 @@ static UV mapes(UV x, uint32 a)
return (UV) val;
}

/* TODO: This is pretty simplistic. */
UV legendre_phi(UV x, UV a) {
UV i, c = (a > 7) ? 7 : a;
UV sum = mapes(x, c);
if (a > c) {
UV p = _XS_nth_prime(c);
UV pa = _XS_nth_prime(a);
for (i = c+1; i <= a; i++) {
p = _XS_next_prime(p);
UV xp = x/p;
if (xp < p) {
while (x < pa) {
a--;
pa = _XS_prev_prime(pa);
}
return (sum - a + i - 1);
}
sum -= legendre_phi(xp, i-1);
}
}
return sum;
}


typedef struct {
sword_t *sieve; /* segment bit mask */
uint8 *word_count; /* bit count in each 64-bit word */
Expand Down
2 changes: 2 additions & 0 deletions lmo.h
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,6 @@

extern UV _XS_LMO_pi(UV n);

extern UV legendre_phi(UV n, UV a);

#endif
5 changes: 1 addition & 4 deletions sieve.h
Original file line number Diff line number Diff line change
Expand Up @@ -111,10 +111,7 @@ static UV prev_prime_in_sieve(const unsigned char* sieve, UV p) {
UV d_ = p/30; \
UV s_, mask_ = 2; \
UV lastd_ = l_/30; \
if (get_prime_cache(l_, &sieve_) < l_) { \
release_prime_cache(sieve_); \
croak("Could not generate sieve for %"UVuf, l_); \
} \
get_prime_cache(l_, &sieve_); \
s_ = sieve_[d_]; \
if (p <= 5) { \
p = (p <= 2) ? 2 : (p <= 3) ? 3 : 5; \
Expand Down
16 changes: 8 additions & 8 deletions t/13-primecount.t
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ plan tests => 0 + 1
+ $use64 * 3 * scalar(keys %pivals64)
+ scalar(keys %intervals)
+ 1
+ 8 + 2*$extra; # prime count specific methods
+ 4 + 2*$extra; # prime count specific methods

ok( eval { prime_count(13); 1; }, "prime_count in void context");

Expand Down Expand Up @@ -156,14 +156,14 @@ sub parse_range {

# Make sure each specific algorithm isn't broken.
SKIP: {
skip "Not XS -- skipping direct primecount tests", 6 unless $isxs;
skip "Not XS -- skipping direct primecount tests", 2 unless $isxs;
# This has to be above SIEVE_LIMIT in lehmer.c and lmo.c or nothing happens.
is(Math::Prime::Util::_XS_lehmer_pi (66123456), 3903023, "XS Lehmer count");
is(Math::Prime::Util::_XS_meissel_pi (66123456), 3903023, "XS Meissel count");
is(Math::Prime::Util::_XS_legendre_pi(66123456), 3903023,"XS Legendre count");
is(Math::Prime::Util::_XS_LMOS_pi (66123456), 3903023, "XS LMOS count");
is(Math::Prime::Util::_XS_LMO_pi (66123456), 3903023, "XS LMO count");
is(Math::Prime::Util::_XS_prime_count(66123456), 3903023, "XS sieve count");
#is(Math::Prime::Util::_XS_lehmer_pi (66123456),3903023,"XS Lehmer count");
#is(Math::Prime::Util::_XS_meissel_pi (66123456),3903023,"XS Meissel count");
#is(Math::Prime::Util::_XS_legendre_pi(66123456),3903023,"XS Legendre count");
#is(Math::Prime::Util::_XS_LMOS_pi (66123456),3903023,"XS LMOS count");
is(Math::Prime::Util::_XS_LMO_pi (66123456), 3903023,"XS LMO count");
is(Math::Prime::Util::_XS_prime_count(66123456), 3903023,"XS sieve count");
}
is(Math::Prime::Util::PP::_lehmer_pi (1456789), 111119, "PP Lehmer count");
is(Math::Prime::Util::PP::_sieve_prime_count(1456789), 111119, "PP sieve count");
Expand Down
2 changes: 1 addition & 1 deletion xt/legendre_phi.t
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ foreach my $a (0 .. $sqrtx+1) {
if ($a > $pcsqrtx) {
is ( $expect, $pcx - $a + 1, "sieved phi($x,$a) = Pi($x) - $a + 1" );
}
my $phixa = Math::Prime::Util::_XS_legendre_phi($x, $a);
my $phixa = Math::Prime::Util::legendre_phi($x, $a);
is( $phixa, $expect, "Legendre phi($x,$a) = $expect" );
}
done_testing();

0 comments on commit 9d3bad1

Please sign in to comment.