This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.10
in repository libmath-prime-util-perl.
commit e7ecf6f95e8ccae03b9aecdf2ad542f8fcaf0a42
Author: Dana Jacobsen <d...@acm.org>
Date: Mon Jul 16 17:29:00 2012 -0600
Strip out the prime_count and nth_prime bounds and approx from C code
---
XS.xs | 6 +-
util.c | 255 ++++-------------------------------------------------------------
util.h | 8 ---
3 files changed, 17 insertions(+), 252 deletions(-)
diff --git a/XS.xs b/XS.xs
index fb0d1a0..d5d52ac 100644
--- a/XS.xs
+++ b/XS.xs
@@ -231,7 +231,7 @@ _XS_factor(IN UV n)
if (!split_success) {
split_success = pbrent_factor(n, tofac_stack+ntofac, 1500)-1;
- if (verbose) { if (split_success) printf("pbrent 1: %lu %lu\n",
tofac_stack[ntofac], tofac_stack[ntofac+1]); else printf("pbrent 0\n"); }
+ if (verbose) { if (split_success) printf("pbrent 1: %"UVuf"
%"UVuf"\n", tofac_stack[ntofac], tofac_stack[ntofac+1]); else printf("pbrent
0\n"); }
}
if (!split_success) {
split_success = squfof_factor(n, tofac_stack+ntofac, 256*1024)-1;
@@ -240,7 +240,7 @@ _XS_factor(IN UV n)
/* Time to do expensive primality check and exit now if it is */
if (!split_success && _XS_is_prime(n)) {
- if (verbose) printf("oops, %lu is prime\n", n);
+ if (verbose) printf("oops, %"UVuf" is prime\n", n);
factored_stack[nfactored++] = n;
n = 1;
break;
@@ -269,7 +269,7 @@ _XS_factor(IN UV n)
n = tofac_stack[ntofac]; /* Set n to the other one */
} else {
/* trial divisions */
- if (verbose) printf("doing trial on %lu\n", (unsigned long)n);
+ if (verbose) printf("doing trial on %"UVuf"\n", n);
UV f = tlim;
UV m = tlim % 30;
UV limit = (UV) (sqrt(n)+0.1);
diff --git a/util.c b/util.c
index 9114315..a636a39 100644
--- a/util.c
+++ b/util.c
@@ -392,148 +392,6 @@ static const unsigned char prime_count_small[] =
16,16,16,16,16,16,17,17,18,18,18,18,18,18,19};
#define NPRIME_COUNT_SMALL
(sizeof(prime_count_small)/sizeof(prime_count_small[0]))
-static const double F1 = 1.0;
-UV _XS_prime_count_lower(UV x)
-{
- double fx, flogx;
- double a = 1.80; /* Dusart 1999, page 14 */
-
- if (x < NPRIME_COUNT_SMALL)
- return prime_count_small[x];
-
- fx = (double)x;
- flogx = log(x);
-
- if (x < 599)
- return (UV) (fx / (flogx-0.7));
-
- if (x < 2700) { a = 0.30; }
- else if (x < 5500) { a = 0.90; }
- else if (x < 19400) { a = 1.30; }
- else if (x < 32299) { a = 1.60; }
- else if (x < 176000) { a = 1.80; }
- else if (x < 315000) { a = 2.10; }
- else if (x < 1100000) { a = 2.20; }
- else if (x < 4500000) { a = 2.31; }
- else if (x <233000000) { a = 2.36; }
-#if BITS_PER_WORD == 32
- else a = 2.32;
-#else
- else if (x < UVCONST( 5433800000)) { a = 2.32; }
- else if (x < UVCONST(60000000000)) { a = 2.15; }
-#endif
-
- return (UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) );
-}
-
-
-UV _XS_prime_count_upper(UV x)
-{
- double fx, flogx;
- double a = 2.51; /* Dusart 1999, page 14 */
-
- if (x < NPRIME_COUNT_SMALL)
- return prime_count_small[x];
-
- fx = (double)x;
- flogx = log(x);
-
- /* This function is unduly complicated. */
-
- if (x < 1621) return (UV) (fx / (flogx-1.048) + F1);
- if (x < 5000) return (UV) (fx / (flogx-1.071) + F1);
- if (x < 15900) return (UV) (fx / (flogx-1.098) + F1);
-
- if (x < 24000) { a = 2.30; }
- else if (x < 59000) { a = 2.48; }
- else if (x < 350000) { a = 2.52; }
- else if (x < 355991) { a = 2.54; }
- else if (x < 356000) { a = 2.51; }
- else if (x < 3550000) { a = 2.50; }
- else if (x < 3560000) { a = 2.49; }
- else if (x < 5000000) { a = 2.48; }
- else if (x < 8000000) { a = 2.47; }
- else if (x < 13000000) { a = 2.46; }
- else if (x < 18000000) { a = 2.45; }
- else if (x < 31000000) { a = 2.44; }
- else if (x < 41000000) { a = 2.43; }
- else if (x < 48000000) { a = 2.42; }
- else if (x <119000000) { a = 2.41; }
- else if (x <182000000) { a = 2.40; }
- else if (x <192000000) { a = 2.395; }
- else if (x <213000000) { a = 2.390; }
- else if (x <271000000) { a = 2.385; }
- else if (x <322000000) { a = 2.380; }
- else if (x <400000000) { a = 2.375; }
- else if (x <510000000) { a = 2.370; }
- else if (x <682000000) { a = 2.367; }
-#if BITS_PER_WORD == 32
- else a = 2.362;
-#else
- else if (x < UVCONST(60000000000)) { a = 2.362; }
-#endif
-
- /*
- * An alternate idea:
- * float alog[23] = { 2.30,2.30,2.30,2.30,2.30,2.30,2.30 ,2.30,2.30,2.30,
- * 2.47,2.49,2.53,2.50,2.49,2.49,2.456,2.44,2.40,2.370,
- * 2.362,2.362,2.362,2.362};
- * float clog[23] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
- * 3, 1, 2, 1, 3, 2, 5, -6, 1, 1,
- * 1, 1, 1, 1};
- * if ((int)flogx < 23) {
- * a = alog[(int)flogx];
- * return ((UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) ) +
clog[(int)flogx] + 0.01);
- * }
- *
- * Another thought is to use more terms in the Li(x) expansion along with
- * a subtraction [Li(x) > Pi(x) for x < 10^316 or so, so for our 64-bit
- * version we should be fine].
- */
-
- return (UV) ( (fx/flogx) * (F1 + F1/flogx + a/(flogx*flogx)) + F1 );
-}
-
-
-UV _XS_prime_count_approx(UV x)
-{
- /*
- * A simple way:
- * return ((prime_count_lower(x) + prime_count_upper(x)) / 2);
- * With the current bounds, this is within ~131k at 10^10 and 436B at 10^19.
- *
- * The logarithmic integral works quite well, with absolute errors of
- * ~3100 at 10^10 and ~100M at 10^19
- *
- * Riemann's R function works even better, with errors of ~1828 at 10^10
- * and 24M at 10^19.
- *
- * Method 10^10 %error 10^19 %error
- * ----------------- ------------ ------------
- * average bounds .01% .0002%
- * li(n) .0007% .00000004%
- * li(n)-li(n^.5)/2 .0004% .00000001%
- * R(n) .0004% .00000001%
- *
- * Getting fancier, one try using Riemann's pi formula:
- * http://trac.sagemath.org/sage_trac/ticket/8135
- */
- double R;
- if (x < NPRIME_COUNT_SMALL)
- return prime_count_small[x];
-
- //R = _XS_LogarithmicIntegral(x) - _XS_LogarithmicIntegral(sqrt(x))/2;
- R = _XS_RiemannR(x);
-
- /* We could add the additional factor:
- * R = R - (1.0 / log(x)) + (M_1_PI * atan(M_PI/log(x)))
- * but it's extraordinarily small, so not worth calculating here.
- */
-
- return (UV)(R+0.5);
-}
-
-
UV _XS_prime_count(UV low, UV high)
{
const unsigned char* cache_sieve;
@@ -614,37 +472,9 @@ static const unsigned short primes_small[] =
409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499};
#define NPRIMES_SMALL (sizeof(primes_small)/sizeof(primes_small[0]))
-/* The nth prime will be greater than or equal to this number */
-UV _XS_nth_prime_lower(UV n)
-{
- double fn = (double) n;
- double flogn, flog2n, lower;
-
- if (n < NPRIMES_SMALL)
- return (n==0) ? 0 : primes_small[n];
-
- flogn = log(n);
- flog2n = log(flogn); /* Note distinction between log_2(n) and log^2(n) */
-
- /* Dusart 1999 page 14, for all n >= 2 */
- lower = fn * (flogn + flog2n - 1.0 + ((flog2n-2.25)/flogn));
-
- /* Watch out for overflow */
- if (lower >= (double)UV_MAX) {
-#if BITS_PER_WORD == 32
- if (n <= UVCONST(203280221)) return UVCONST(4294967291);
-#else
- if (n <= UVCONST(425656284035217743)) return UVCONST(18446744073709551557);
-#endif
- croak("nth_prime_lower(%"UVuf") overflow", n);
- }
-
- return (UV) lower;
-}
-
-
+/* Note: We're keeping this here because we use it for nth_prime */
/* The nth prime will be less or equal to this number */
-UV _XS_nth_prime_upper(UV n)
+static UV _XS_nth_prime_upper(UV n)
{
double fn = (double) n;
double flogn, flog2n, upper;
@@ -655,75 +485,17 @@ UV _XS_nth_prime_upper(UV n)
flogn = log(n);
flog2n = log(flogn); /* Note distinction between log_2(n) and log^2(n) */
- if (n >= 39017)
- upper = fn * ( flogn + flog2n - 0.9484 ); /* Dusart 1999 page 14*/
- else if (n >= 27076)
- upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.80)/flogn)); /*Dusart
1999*/
- else if (n >= 7022)
- upper = fn * ( flogn + 0.9385 * flog2n ); /* Robin 1983 */
+ if (n >= 688383) /* Dusart 2010 page 2 */
+ upper = fn * (flogn + flog2n - 1.0 + ((flog2n-2.00)/flogn));
+ else if (n >= 178974) /* Dusart 2010 page 7 */
+ upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.95)/flogn));
+ else if (n >= 27076) /* Dusart 1999 page 14 */
+ upper = fn * (flogn + flog2n - 1.0 + ((flog2n-1.80)/flogn));
+ else if (n >= 6) /* Modified from Robin 1983 for 6-27075 _only_ */
+ upper = fn * ( flogn + 0.6000 * flog2n );
else
upper = fn * ( flogn + flog2n );
- /* Watch out for overflow */
- if (upper >= (double)UV_MAX) {
-#if BITS_PER_WORD == 32
- if (n <= UVCONST(203280221)) return UVCONST(4294967291);
-#else
- if (n <= UVCONST(425656284035217743)) return UVCONST(18446744073709551557);
-#endif
- croak("nth_prime_upper(%"UVuf") overflow", n);
- }
-
- return (UV) ceil(upper);
-}
-
-
-UV _XS_nth_prime_approx(UV n)
-{
- double fn, flogn, flog2n, order, approx;
-
- if (n < NPRIMES_SMALL)
- return primes_small[n];
-
- /* This isn't too bad:
- * return ((nth_prime_lower(n) + nth_prime_upper(n)) / 2);
- */
-
- fn = (double) n;
- flogn = log(n);
- flog2n = log(flogn);
-
- /* Cipolla 1902:
- * m=0 fn * ( flogn + flog2n - 1 );
- * m=1 + ((flog2n - 2)/flogn) );
- * m=2 - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
- * + O((flog2n/flogn)^3)
- *
- * Shown in Dusart 1999 page 12, as well as other sources such as:
- * http://www.emis.de/journals/JIPAM/images/153_02_JIPAM/153_02.pdf
- * where the main issue you run into is that you're doing polynomial
- * interpolation, so it oscillates like crazy with many high-order terms.
- * Hence I'm leaving it at m=2.
- */
- approx = fn * ( flogn + flog2n - 1
- + ((flog2n - 2)/flogn)
- - (((flog2n*flog2n) - 6*flog2n + 11) / (2*flogn*flogn))
- );
-
- /* Apply a correction to help keep values close */
- order = flog2n/flogn;
- order = order*order*order * fn;
-
- if (n < 259) approx += 10.4 * order;
- else if (n < 775) approx += 7.52 * order;
- else if (n < 1271) approx += 5.6 * order;
- else if (n < 2000) approx += 5.2 * order;
- else if (n < 4000) approx += 4.3 * order;
- else if (n < 12000) approx += 3.0 * order;
- else if (n <150000) approx += 2.1 * order;
- else if (n <200000000) approx += 0.0 * order;
- else approx += -0.010 * order; /* -0.025 is better */
-
/* For all three analytical functions, it is possible that for a given valid
* input, we will not be able to return an output that fits in the UV type.
* For example, if they ask for the 203280222nd prime, we should return
@@ -734,16 +506,17 @@ UV _XS_nth_prime_approx(UV n)
* This should maintain the invariant:
* nth_prime_lower(n) <= nth_prime(n) <= nth_prime_upper(n)
*/
- if (approx >= (double)UV_MAX) {
+ /* Watch out for overflow */
+ if (upper >= (double)UV_MAX) {
#if BITS_PER_WORD == 32
if (n <= UVCONST(203280221)) return UVCONST(4294967291);
#else
if (n <= UVCONST(425656284035217743)) return UVCONST(18446744073709551557);
#endif
- croak("nth_prime_approx(%"UVuf") overflow", n);
+ croak("nth_prime_upper(%"UVuf") overflow", n);
}
- return (UV) (approx + 0.5);
+ return (UV) ceil(upper);
}
diff --git a/util.h b/util.h
index 469a07a..69873c4 100644
--- a/util.h
+++ b/util.h
@@ -12,14 +12,6 @@ extern UV _XS_prev_prime(UV x);
extern UV _XS_prime_count(UV low, UV high);
extern UV _XS_nth_prime(UV x);
-/* These have been moved into the main Util.pm */
-extern UV _XS_prime_count_lower(UV x);
-extern UV _XS_prime_count_upper(UV x);
-extern UV _XS_prime_count_approx(UV x);
-extern UV _XS_nth_prime_lower(UV n);
-extern UV _XS_nth_prime_upper(UV x);
-extern UV _XS_nth_prime_approx(UV x);
-
extern double _XS_ExponentialIntegral(double x);
extern double _XS_LogarithmicIntegral(double x);
extern double _XS_RiemannR(double x);
--
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
