[libmath-prime-util-perl] 13/16: Small fixes
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:44:27 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list