This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.23 in repository libmath-prime-util-perl.
commit 7be91a6ad6e358ec9de02359c8c50397855e8f2e Author: Dana Jacobsen <d...@acm.org> Date: Sun Mar 3 18:19:53 2013 -0800 Change type of mobius return. Some different ranged algorithms in comments --- XS.xs | 2 +- util.c | 108 ++++++++++++++++++++++++++++++++++++++++++++++++----------------- util.h | 4 +-- 3 files changed, 83 insertions(+), 31 deletions(-) diff --git a/XS.xs b/XS.xs index 8892164..bd60751 100644 --- a/XS.xs +++ b/XS.xs @@ -415,7 +415,7 @@ _XS_moebius(IN UV lo, IN UV hi = 0) UV i; PPCODE: if (hi != lo && hi != 0) { /* mobius in a range */ - IV* mu = _moebius_range(lo, hi); + char* mu = _moebius_range(lo, hi); MPUassert( mu != 0, "_moebius_range returned 0" ); EXTEND(SP, hi-lo+1); for (i = lo; i <= hi; i++) diff --git a/util.c b/util.c index b041fb6..66525c1 100644 --- a/util.c +++ b/util.c @@ -651,46 +651,98 @@ UV _XS_nth_prime(UV n) /* Return an IV array with lo-hi+1 elements. mu[k-lo] = µ(k) for k = lo .. hi. * It is the callers responsibility to call Safefree on the result. */ -IV* _moebius_range(UV lo, UV hi) +#define PGTLO(p,lo) (p >= lo) ? p : (p*(lo/p) + ((lo%p)?p:0)) +char* _moebius_range(UV lo, UV hi) { - IV* mu; - UV i, p, sqrtn; + char* mu; + UV i, p; + UV sqrtn = isqrt(hi); +#if 1 + IV* A; /* This implementation follows that of Deléglise & Rivat (1996), which is - * a segmented version of Lioen & van de Lune (1994). + * a segmented version of Lioen & van de Lune (1994). Downside is that it + * uses an array of IV's, so too much memory. Pawleicz (2011) shows a small + * variation but seems to just be a rearrangement -- there is no time or + * space difference on my machines. (TODO: but maybe for hi > 2^63?) */ - sqrtn = isqrt(hi); - New(0, mu, hi-lo+1, IV); - if (mu == 0) + New(0, A, hi-lo+1, IV); + if (A == 0) croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1); for (i = lo; i <= hi; i++) - mu[i-lo] = 1; - if (lo == 0) mu[0] = 0; + A[i-lo] = 1; prime_precalc(sqrtn); for (p = 2; p <= sqrtn; p = _XS_next_prime(p)) { - i = p*p; - if (i < lo) - i = i*(lo/i) + ( (lo%i) ? i : 0 ); - while (i <= hi) { + UV p2 = p*p; + for (i = PGTLO(p2, lo); i <= hi; i += p2) + A[i-lo] = 0; + for (i = PGTLO(p, lo); i <= hi; i += p) + A[i-lo] *= -(IV)p; + } + New(0, mu, hi-lo+1, char); + if (mu == 0) + croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1); + memset(mu, 0, hi-lo+1); + for (i = lo; i <= hi; i++) { + IV a = A[i-lo]; + if (a != 0) + mu[i-lo] = (a != i && -a != i) ? (a<0) - (a>0) + : (a>0) - (a<0); + } + Safefree(A); +#endif +#if 0 + /* Simple char method, Needs way too many primes */ + New(0, mu, hi-lo+1, char); + if (mu == 0) + croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1); + memset(mu, 1, hi-lo+1); + if (lo == 0) mu[0] = 0; + prime_precalc( _XS_nth_prime_upper(hi) ); + for (p = 2; p <= hi; p = _XS_next_prime(p)) { + UV p2 = p*p; + for (i = PGTLO(p2, lo); i <= hi; i += p2) mu[i-lo] = 0; - i += p*p; - } - i = p; - if (i < lo) - i = i*(lo/i) + ( (lo%i) ? i : 0 ); - while (i <= hi) { - mu[i-lo] *= -(IV)p; - i += p; - } + for (i = PGTLO(p, lo); i <= hi; i += p) + mu[i-lo] = -mu[i-lo]; } +#endif +#if 0 + /* Kuznetsov's transform of Deléglise & Rivat (1996) into logs. + * (1) I'm using the log function, which should be fixed (easy). + * (2) it doesn't work. try 64-101 vs. 64 vs. 100. + */ + unsigned char* A; + New(0, A, hi-lo+1, unsigned char); + if (A == 0) + croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1); + memset(A, 0, hi-lo+1); + if (sqrtn*sqrtn != hi) sqrtn++; /* ceil sqrtn */ + prime_precalc(sqrtn); + for (p = 2; p <= sqrtn; p = _XS_next_prime(p)) { + UV p2 = p*p; + unsigned char l = 1 | (unsigned char) ( log(p)/log(2) ); + for (i = PGTLO(p, lo); i <= hi; i += p) + A[i-lo] += l; + for (i = PGTLO(p2, lo); i <= hi; i += p2) + A[i-lo] |= 0x80; + } + + New(0, mu, hi-lo+1, char); + if (mu == 0) + croak("Could not get memory for %"UVuf" moebius results\n", hi-lo+1); for (i = lo; i <= hi; i++) { - IV m = mu[i-lo]; - if (m != 0) { - if (m != i && -m != i) m *= -1; - mu[i-lo] = (m>0) - (m<0); - } + unsigned char a = A[i-lo]; + unsigned char log2i = (unsigned char) ( log(i)/log(2) ); + //printf("i = %lu a = %lu log2i = %lu\n", i, a, log2i); + if (a & 0x80) { mu[i-lo] = 0; } + else if (a >= log2i) { mu[i-lo] = 1 - 2*(a&1); } + else { mu[i-lo] = -1 + 2*(a&1); } } + if (lo == 0) mu[0] = 0; + Safefree(A); +#endif return mu; } @@ -720,7 +772,7 @@ IV _XS_mertens(UV n) { /* Deléglise and Rivat (1996) using u = n^1/2 and unsegmented. */ /* Very simple, but they use u = n^1/3 and segment */ UV u, i, m, nmk; - IV* mu; + char* mu; IV* M; IV sum; diff --git a/util.h b/util.h index 9d9821e..a7b56e4 100644 --- a/util.h +++ b/util.h @@ -14,8 +14,8 @@ extern UV _XS_prev_prime(UV x); extern UV _XS_prime_count(UV low, UV high); extern UV _XS_nth_prime(UV x); -extern IV* _moebius_range(UV low, UV high); -extern IV _XS_mertens(UV n); +extern char* _moebius_range(UV low, UV high); +extern IV _XS_mertens(UV n); extern double _XS_chebyshev_theta(UV n); extern double _XS_chebyshev_psi(UV n); -- 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