This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.06
in repository libmath-prime-util-perl.

commit a17975bfbd8ba5a6d4fd2de04d721d563c648888
Author: Dana Jacobsen <d...@acm.org>
Date:   Mon Jun 11 21:15:57 2012 -0600

    Small fixes
---
 TODO                   | 10 +++++----
 XS.xs                  |  3 +++
 factor.c               |  2 +-
 lib/Math/Prime/Util.pm |  7 +++++++
 sieve.c                |  7 ++++++-
 util.c                 | 55 ++++++++------------------------------------------
 6 files changed, 31 insertions(+), 53 deletions(-)

diff --git a/TODO b/TODO
index 2661193..6f36135 100644
--- a/TODO
+++ b/TODO
@@ -17,8 +17,10 @@
 
 - speed up random_prime for large numbers
 
-- test x64 asm if Perl is built for 32-bit
-
-- Clean up prime_count ranges.  Write tests.
-
 - better prime count upper/lower bounds
+
+- Thread-safety.  Change all uses of sieve and segment to include a free call.
+  For segment we can:
+     1) mutex  (holy slowdown, batman)
+     2) use one cache, and malloc/free if the cache isn't available.
+  For sieve I think we'll need a counting semaphore.
diff --git a/XS.xs b/XS.xs
index 017db4f..eac8c57 100644
--- a/XS.xs
+++ b/XS.xs
@@ -20,6 +20,9 @@ prime_precalc(IN UV n)
 void
 prime_memfree()
 
+void
+_prime_memfreeall()
+
 UV
 prime_count(IN UV low, IN UV high = 0)
   CODE:
diff --git a/factor.c b/factor.c
index 70192a5..4f49763 100644
--- a/factor.c
+++ b/factor.c
@@ -189,7 +189,7 @@ int trial_factor(UV n, UV *factors, UV maxtrial)
  */
 static int is_perfect_square(UV n, UV* sqrtn)
 {
-  UV m, lm;
+  UV m;  /* lm */
   m = n & 127;
   if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a)  return 0;
   /* 82% of non-squares rejected here */
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 330c4fc..57122bd 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -35,6 +35,9 @@ BEGIN {
     croak "XS Code not available: $@";
   }
 }
+END {
+  _prime_memfreeall;
+}
 
 
 my $_maxparam  = (_maxbits == 32) ? 4294967295 : 18446744073709551615;
@@ -702,6 +705,10 @@ Perl versions earlier than 5.8.0 have issues with 64-bit.  
The test suite will
 try to determine if your Perl is broken.  This will show up in factoring tests.
 Perl 5.6.2 32-bit works fine, as do later versions with 32-bit and 64-bit.
 
+Because static caches are used, many functions are not threadsafe.  If you
+use C<prime_precalc> and all calls have inputs smaller than that number,
+then only C<nth_prime> is problematic.  This will be addressed in a later
+implementation.
 
 
 =head1 PERFORMANCE
diff --git a/sieve.c b/sieve.c
index 0636f4b..14092ce 100644
--- a/sieve.c
+++ b/sieve.c
@@ -67,7 +67,7 @@ void prime_precalc(UV n)
   /* TODO: should we prealloc the segment here? */
 }
 
-void prime_memfree(void)
+void _prime_memfreeall(void)
 {
   if (prime_cache_sieve != 0)
       free(prime_cache_sieve);
@@ -75,6 +75,11 @@ void prime_memfree(void)
   prime_cache_size = 0;
 
   free_prime_segment();
+}
+
+void prime_memfree(void)
+{
+  _prime_memfreeall();
 
   prime_precalc(0);
 }
diff --git a/util.c b/util.c
index b59ec1d..c7a704e 100644
--- a/util.c
+++ b/util.c
@@ -412,13 +412,11 @@ UV prime_count_lower(UV x)
   else if (x <  1100000)  { a = 2.20; }
   else if (x <  4500000)  { a = 2.31; }
   else if (x <233000000)  { a = 2.36; }
-  else if (x <240000000)  { a = 2.32; }
 #if BITS_PER_WORD == 32
   else a = 2.32;
 #else
   else if (x < UVCONST( 5433800000)) { a = 2.32; }
-  else if (x < UVCONST(10000000000)) { a = 2.15; }
-  else a = 1.80;
+  else if (x < UVCONST(50000000000)) { a = 2.15; }
 #endif
 
   return (UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) );
@@ -468,8 +466,7 @@ UV prime_count_upper(UV x)
 #if BITS_PER_WORD == 32
   else a = 2.362;
 #else
-  else if (x < UVCONST(10000000000)) { a = 2.362; }
-  else a = 2.51;
+  else if (x < UVCONST(50000000000)) { a = 2.362; }
 #endif
 
   /*
@@ -527,7 +524,7 @@ UV prime_count(UV low, UV high)
 {
   const unsigned char* cache_sieve;
   unsigned char* segment;
-  UV segbase, segment_size, low_d, high_d;
+  UV segment_size, low_d, high_d;
   UV count = 0;
 
   if ( (low <= 2) && (high < NPRIME_COUNT_SMALL) )
@@ -542,7 +539,6 @@ UV prime_count(UV low, UV high)
 
   low_d = low/30;
   high_d = high/30;
-  //printf("  our range %lu - %lu lies in bytes %lu - %lu\n", low, high, 
low_d, high_d);
 
   /* Count full bytes only -- no fragments from primary cache */
   segment_size = get_prime_cache(0, &cache_sieve) / 30;
@@ -551,12 +547,10 @@ UV prime_count(UV low, UV high)
     UV endp = (high_d >= (UV_MAX/30))  ?  UV_MAX-2  :  30*high_d+29;
     segment_size = get_prime_cache( sqrt(endp) + 1 , &cache_sieve) / 30;
   }
-  //printf("  primary cache is %lu bytes.  From %lu to %lu.  Size=%lu, 
count=%lu\n", segment_size, 0, (segment_size-1)*30+29, get_prime_cache_size(), 
count_zero_bits(cache_sieve, segment_size-1));
 
-  if (1 && low_d <= segment_size) {
+  if (low_d <= segment_size) {
     /* Count all the primes in the primary cache in our range */
     count += count_segment_ranged(cache_sieve, segment_size, low, high);
-    //printf("  counted in pcache, count (%lu,%lu) is %lu\n", low, high, 
count);
 
     if (high_d <= segment_size)
       return count;
@@ -578,15 +572,12 @@ UV prime_count(UV low, UV high)
     UV range_d = seghigh_d - low_d + 1;
     UV seglow  = (low_d*30 < low) ? low : low_d*30;
     UV seghigh = (seghigh_d == high_d) ? high : (seghigh_d*30+29);
-    //printf("  seg d: %lu - %lu   seg p: %lu - %lu\n", low_d, seghigh_d, 
seglow, seghigh);
 
-    //printf("  sieve from bytes %lu to %lu\n", low_d, seghigh_d);
     if (sieve_segment(segment, low_d, seghigh_d) == 0) {
-      croak("Could not segment sieve from %"UVuf" to %"UVuf, segbase+1, 
30*(segbase+segment_size)+29);
+      croak("Could not segment sieve from %"UVuf" to %"UVuf, low_d*30+1, 
30*seghigh_d+29);
       break;
     }
 
-    //printf("  count from %lu to %lu\n", seglow - low_d*30, seghigh - 
low_d*30);
     count += count_segment_ranged(segment, segment_size, seglow - low_d*30, 
seghigh - low_d*30);
 
     low_d += range_d;
@@ -822,7 +813,7 @@ static double const li2 = 
1.045163780117492784844588889194613136522615578151;
 
 double ExponentialIntegral(double x) {
   double const tol = 0.0000000001;
-  double val, term, x_to_the_n, fact_n;
+  double val, term, fact_n;
   double y, t;
   double sum = 0.0;
   double c = 0.0;
@@ -832,13 +823,6 @@ double ExponentialIntegral(double x) {
 
   if (x < 0) {
     /* continued fraction */
-#if 0
-    val = 0;
-    for (n = 50; n >= 1; n--) {
-      val = -n * n / (2 * n + 1 - x + val);
-    }
-    val = -exp(x) / (1 - x + val);
-#else
     double old;
     double lc = 0;
     double ld = 1 / (1 - x);
@@ -851,7 +835,6 @@ double ExponentialIntegral(double x) {
       if ( fabs(val-old) <= tol*fabs(val) )
         break;
     }
-#endif
   } else if (x < -log(tol)) {
     /* Convergent series */
     fact_n = 1;
@@ -867,29 +850,9 @@ double ExponentialIntegral(double x) {
       if (term < tol) break;
     }
     val = sum;
-    if (term > tol) printf("Tolerance not met after %d rounds for Ei(%lf)\n", 
n, x);
   } else {
     /* Asymptotic divergent series */
-#if 0
-    double last_term, sum;
-    last_term = 1.0;
-    x_to_the_n = 1;
-    fact_n = 1;
-
-    y = 1.0-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
-
-    for (n = 1; n <= 100; n++) {
-      x_to_the_n *= x;
-      fact_n *= n;
-      term = fact_n / x_to_the_n;
-      if (term > last_term) break;
-      last_term = term;
-      y = term-c;  t = sum+y;  c = (t-sum)-y;  sum = t;
-      /* printf("A  after adding %.8lf, sum = %.8lf\n", term, sum); */
-    }
-    val = (exp(x) / x) * sum - euler_mascheroni;
-#else
-    double last_term, sum;
+    double last_term;
 
     val = exp(x) / x;
     term = 1.0;
@@ -909,7 +872,6 @@ double ExponentialIntegral(double x) {
       }
     }
     val *= sum;
-#endif
   }
 
   return val;
@@ -1011,7 +973,7 @@ double RiemannR(double x) {
   part_term = 1;
 
   /* Do small k with zetas from table */
-  for (k = 1; k <= NPRECALC_ZETA; k++) {
+  for (k = 1; k <= (int)NPRECALC_ZETA; k++) {
     zeta = riemann_zeta_table[k+1-2];
     part_term *= flogx / k;
     term = part_term / (k * zeta);
@@ -1029,6 +991,5 @@ double RiemannR(double x) {
     /* printf("R  after adding %.8lf, sum = %.8lf\n", term, sum); */
   }
 
-
   return sum;
 }

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

Reply via email to