[libmath-prime-util-perl] 18/72: Add Pk_2 function, simplify some code

```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 040a785a0e86a0edfac00ed5dffb6816a6f32cb5
Author: Dana Jacobsen <d...@acm.org>
Date:   Sat Aug 31 10:05:31 2013 -0700

Add Pk_2 function, simplify some code
---
lehmer.c | 129 ++++++++++++++++++++-------------------------------------------
1 file changed, 40 insertions(+), 89 deletions(-)

diff --git a/lehmer.c b/lehmer.c
index 48958d0..3cfa8a7 100644
--- a/lehmer.c
+++ b/lehmer.c
@@ -5,7 +5,7 @@

/* Below this size, just sieve. */
#define SIEVE_LIMIT  60000000
-#define MAX_PHI_MEM  (128*1024*1024)
+#define MAX_PHI_MEM  (256*1024*1024)

/*****************************************************************************
*
@@ -650,15 +650,13 @@ static UV phi(UV x, UV a)
if (a > 7) {
if (verbose > 0) printf("CUT TO SMALL PHI AT a = %lu\n", a);
#ifdef _OPENMP
schedule(dynamic, 16)
+    #pragma omp parallel for reduction(+: sum) num_threads(2)
schedule(dynamic, 16)
#endif
for (i = 0; i < a1.n; i++)
sum += arr[i].c * phi_small( arr[i].v, a, primes );
} else {
#ifdef _OPENMP
schedule(dynamic, 16)
+    #pragma omp parallel for reduction(+: sum) schedule(dynamic, 16)
#endif
for (i = 0; i < a1.n; i++)
sum += arr[i].c * mapes7( arr[i].v );
@@ -669,7 +667,32 @@ static UV phi(UV x, UV a)
}

+extern UV _XS_meissel_pi(UV n);
+static UV Pk_2(UV n, UV a)
+{
+  UV b, lastprime, lastpc, lastw, lastwpc, i, P2;
+  const UV* primes; /* small prime cache, first b=pi(z)=pi(sqrt(n)) */

+  b = _XS_meissel_pi(isqrt(n));        /* b = pi(floor(n^1/2)) */
+  lastprime = b*SIEVE_MULT+1;
+  primes = generate_small_primes(lastprime);
+  lastpc = primes[lastprime];
+
+  /* Ensure we have a large enough base sieve */
+  prime_precalc(isqrt(n / primes[a+1]));
+
+  P2 = lastw = lastwpc = 0;
+  for (i = b; i > a; i--) {
+    UV w = n / primes[i];
+    lastwpc = (w <= lastpc) ? bs_prime_count(w, primes, lastprime)
+                            : lastwpc + _XS_prime_count(lastw+1, w);
+    lastw = w;
+    P2 += lastwpc;
+  }
+  P2 -= ((b+a-2) * (b-a+1) / 2) - a + 1;
+  Safefree(primes);
+  return P2;
+}

/* Legendre's method.  Interesting and a good test for phi(x,a), but Lehmer's
@@ -677,21 +700,12 @@ static UV phi(UV x, UV a)
UV _XS_legendre_pi(UV n)
{
UV a, phina;
-  DECLARE_TIMING_VARIABLES;
if (n < SIEVE_LIMIT)
return _XS_prime_count(2, n);

-  if (verbose > 0) printf("legendre %lu stage 1: calculate a \n", n);
-  TIMING_START;
a = _XS_legendre_pi(isqrt(n));
-  TIMING_END_PRINT("stage 1")
-
-  if (verbose > 0) printf("legendre %lu stage 2: phi(x,a) (a=%lu)\n", n, a);
-  TIMING_START;
phina = phi(n, a);
-  if (verbose > 0) printf("phicache size: %lu\n", (unsigned
long)phicache_size());
phicache_free();
-  TIMING_END_PRINT("stage 2")
return phina + a - 1;
}

@@ -699,52 +713,14 @@ UV _XS_legendre_pi(UV n)
/* Meissel's method. */
UV _XS_meissel_pi(UV n)
{
-  UV a, b, c, sum, i, lastprime, lastpc, lastw, lastwpc;
-  const UV* primes = 0;  /* small prime cache */
-  DECLARE_TIMING_VARIABLES;
+  UV a, sum;
if (n < SIEVE_LIMIT)
return _XS_prime_count(2, n);

-  if (verbose > 0) printf("meissel %lu stage 1: calculate a,b,c \n", n);
-  TIMING_START;
a = _XS_meissel_pi(icbrt(n));        /* a = floor(n^1/3) */
-  b = _XS_meissel_pi(isqrt(n));        /* b = floor(n^1/2) */
-  c = a;                               /* c = a            */
-  TIMING_END_PRINT("stage 1")

-  if (verbose > 0) printf("meissel %lu stage 2: phi(x,a) (a=%lu b=%lu
c=%lu)\n", n, a, b, c);
-  TIMING_START;
-  sum = phi(n, a) + ((b+a-2) * (b-a+1) / 2);
-  if (verbose > 1) printf("phicache size: %lu\n", (unsigned
long)phicache_size());
+  sum = phi(n, a) + a - 1 - Pk_2(n, a);
phicache_free();
-  if (verbose > 0) printf("phi(%lu,%lu) = %lu.  sum = %lu\n", n, a, sum -
((b+a-2) * (b-a+1) / 2), sum);
-  TIMING_END_PRINT("phi(x,a)")
-
-  lastprime = b*SIEVE_MULT+1;
-  if (verbose > 0) printf("meissel %lu stage 3: %lu small primes\n", n,
lastprime);
-  TIMING_START;
-  primes = generate_small_primes(lastprime);
-  if (primes == 0) croak("Error generating primes.\n");
-  lastpc = primes[lastprime];
-  TIMING_END_PRINT("small primes")
-
-  prime_precalc(isqrt(n / primes[a+1]));
-  prime_precalc( (UV) pow(n, 2.0/5.0) );  /* Sieve more for speed */
-
-  if (verbose > 0) printf("meissel %lu stage 4: loop %lu to %lu, pc to %lu\n",
n, a+1, b, n/primes[a+1]);
-  TIMING_START;
-  /* Reverse the i loop so w increases.  Count w in segments. */
-  lastw = 0;
-  lastwpc = 0;
-  for (i = b; i > a; i--) {
-    UV w = n / primes[i];
-    lastwpc = (w <= lastpc) ? bs_prime_count(w, primes, lastprime)
-                            : lastwpc + _XS_prime_count(lastw+1, w);
-    lastw = w;
-    sum = sum - lastwpc;
-  }
-  TIMING_END_PRINT("stage 4")
-  Safefree(primes);
return sum;
}

@@ -831,7 +807,7 @@ UV _XS_lehmer_pi(UV n)
* particularly fast. */
UV _XS_LMO_pi(UV n)
{
-  UV n12, n13, a, b, sum, i, j, lastprime, lastpc, lastw, lastwpc, P2, S1, S2;
+  UV n12, n13, a, b, sum, i, j, k, lastprime, P2, S1, S2;
const UV* primes = 0;  /* small prime cache */
signed char* mu = 0;   /* moebius to n^1/3 */
UV*   lpf = 0;         /* least prime factor to n^1/3 */
@@ -839,27 +815,18 @@ UV _XS_LMO_pi(UV n)
if (n < SIEVE_LIMIT)
return _XS_prime_count(2, n);

-  if (verbose > 0) printf("LMO %lu stage 1: calculate pi(n^1/3) \n", n);
-  TIMING_START;
n13 = icbrt(n);
n12 = isqrt(n);
a = _XS_lehmer_pi(n13);
b = _XS_lehmer_pi(n12);
-  TIMING_END_PRINT("stage 1")

lastprime = b*SIEVE_MULT+1;
if (lastprime < n13) lastprime = n13;
-  if (verbose > 0) printf("LMO %lu stage 2: %lu small primes\n", n, lastprime);
-  TIMING_START;
primes = generate_small_primes(lastprime);
if (primes == 0) croak("Error generating primes.\n");
-  lastpc = primes[lastprime];
-  TIMING_END_PRINT("small primes")

-  if (verbose > 0) printf("LMO %lu stage 3: moebius and lpf\n", n);
-  TIMING_START;
New(0, mu, n13+1, signed char);
-  memset(mu, 1, sizeof(char) * (n13+1));
+  memset(mu, 1, sizeof(signed char) * (n13+1));
Newz(0, lpf, n13+1, UV);
mu[0] = 0;
for (i = 1; i <= n13; i++) {
@@ -874,41 +841,25 @@ UV _XS_LMO_pi(UV n)
mu[j] = 0;
}
lpf[1] = UV_MAX;  /* Set lpf[1] to max */
-  TIMING_END_PRINT("moebius and lpf")

/* Thanks to Kim Walisch for help with the S1+S2 calculations. */
-  if (verbose > 0) printf("LMO %lu stage 4: S1 and S2\n", n);
-  TIMING_START;
-  UV k = (a < 7) ? a : 7;
+  k = (a < 7) ? a : 7;
S1 = 0;
S2 = 0;
for (i = 1; i <= n13; i++)
if (lpf[i] > primes[k])
-      S1 += mu[i] * mapes(n/i, k);
+      S1 += mu[i] * phi_small(n/i, k, primes);
for (i = k; i+1 < a; i++)
+#ifdef _OPENMP
+    #pragma omp parallel for reduction(+: S2) schedule(dynamic, 16)
+#endif
for (j = (n13/primes[i+1])+1; j <= n13; j++)
if (lpf[j] > primes[i+1])
-        S2 -= mu[j] * phi_small(n / (j*primes[i+1]), i, primes);
-  if (verbose > 0) printf("phicache size: %lu\n", (unsigned
long)phicache_size());
+        S2 += -mu[j] * phi_small(n / (j*primes[i+1]), i, primes);
phicache_free();
-  TIMING_END_PRINT("S1 and S2")

-  if (verbose > 0) printf("LMO %lu stage 5: P2\n", n);
-  TIMING_START;
-  prime_precalc(isqrt(n / primes[a+1]));
-  /* Reverse the i loop so w increases.  Count w in segments. */
-  lastw = 0;
-  lastwpc = 0;
-  P2 = 0;
-  for (i = b; i > a; i--) {
-    UV w = n / primes[i];
-    lastwpc = (w <= lastpc) ? bs_prime_count(w, primes, lastprime)
-                            : lastwpc + _XS_prime_count(lastw+1, w);
-    lastw = w;
-    P2 += lastwpc;
-  }
-  P2 -= ((b+a-2) * (b-a+1) / 2) - a + 1;
-  TIMING_END_PRINT("stage 5 (P2)")
+  P2 = Pk_2(n, a);
+
/* printf("S1 = %lu\nS2 = %lu\na  = %lu\nP2 = %lu\n", S1, S2, a, P2); */
sum = (S1 + S2) + a - 1 - P2;
Safefree(primes);

--
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
```