[libmath-prime-util-perl] 25/181: Mapes => table for small phi
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:51:03 UTC 2015
This is an automated email from the git hooks/post-receive script.
ppm-guest pushed a commit to annotated tag v0.36
in repository libmath-prime-util-perl.
commit fa9264cb736a994954fc53837799b008a50fce6f
Author: Dana Jacobsen <dana at acm.org>
Date: Thu Dec 19 18:16:39 2013 -0800
Mapes => table for small phi
---
Changes | 4 +++
lehmer.c | 122 ++++++++++++++++++++++++++-------------------------------------
lmo.c | 2 +-
3 files changed, 55 insertions(+), 73 deletions(-)
diff --git a/Changes b/Changes
index 547926b..49cf4b5 100644
--- a/Changes
+++ b/Changes
@@ -31,6 +31,10 @@ Revision history for Perl module Math::Prime::Util
is mostly in preparation for a alldivisors function, and moving a few
more functions to XS/Perl from Perl/XS to cut overhead.
+ - Switch from mapes to a cached primorial/totient small phi method.
+ Really only affects LMOS and Legendre, but it's pretty significant
+ for them. Thanks to Kim Walisch for pointing this out.
+
0.35 2013-12-08
diff --git a/lehmer.c b/lehmer.c
index 89edbc7..99ce156 100644
--- a/lehmer.c
+++ b/lehmer.c
@@ -273,67 +273,49 @@ static UV bs_prime_count(uint32_t n, uint32_t const* const primes, uint32_t last
#define FAST_DIV(x,y) \
( ((x) <= 4294967295U) ? (uint32_t)(x)/(uint32_t)(y) : (x)/(y) )
-/* Use Mapes' method to calculate phi(x,a) for small a. This is really
- * convenient and a little Perl script will spit this code out for whatever
- * limit we select. It gets unwieldy with large a values.
- */
-static UV mapes(UV x, UV a)
-{
- IV val;
- if (a == 0) return x;
- if (a == 1) return x-x/2;
- val = x-x/2-x/3+x/6;
- if (a >= 3) val += 0-x/5+x/10+x/15-x/30;
- if (a >= 4) val += 0-x/7+x/14+x/21-x/42+x/35-x/70-x/105+x/210;
- if (a >= 5) val += 0-x/11+x/22+x/33-x/66+x/55-x/110-x/165+x/330+x/77-x/154-x/231+x/462-x/385+x/770+x/1155-x/2310;
- if (a >= 6) val += 0-x/13+x/26+x/39-x/78+x/65-x/130-x/195+x/390+x/91-x/182-x/273+x/546-x/455+x/910+x/1365-x/2730+x/143-x/286-x/429+x/858-x/715+x/1430+x/2145-x/4290-x/1001+x/2002+x/3003-x/6006+x/5005-x/10010-x/15015+x/30030;
- if (a >= 7) val += 0-x/17+x/34+x/51-x/102+x/85-x/170-x/255+x/510+x/119-x/238-x/357+x/714-x/595+x/1190+x/1785-x/3570+x/187-x/374-x/561+x/1122-x/935+x/1870+x/2805-x/5610-x/1309+x/2618+x/3927-x/7854+x/6545-x/13090-x/19635+x/39270+x/221-x/442-x/663+x/1326-x/1105+x/2210+x/3315-x/6630-x/1547+x/3094+x/4641-x/9282+x/7735-x/15470-x/23205+x/46410-x/2431+x/4862+x/7293-x/14586+x/12155-x/24310-x/36465+x/72930+x/17017-x/34034-x/51051+x/102102-x/85085+x/170170+x/255255-x/510510;
- return (UV) val;
-}
-
-#define mapes7(x) (((x) <= 4294967295U) ? mapes7_32(x) : mapes7_64(x))
-static UV mapes7_64(UV x) { /* A tiny bit faster setup for a=7 */
- UV val = x-x/2-x/3-x/5+x/6-x/7+x/10-x/11-x/13+x/14+x/15-x/17+x/21+x/22+x/26
- -x/30+x/33+x/34+x/35+x/39-x/42+x/51+x/55+x/65-x/66-x/70+x/77-x/78
- +x/85+x/91-x/102-x/105-x/110+x/119-x/130+x/143-x/154-x/165-x/170
- -x/182+x/187-x/195+x/210+x/221-x/231-x/238-x/255-x/273-x/286+x/330
- -x/357-x/374-x/385+x/390-x/429-x/442-x/455+x/462+x/510+x/546-x/561
- -x/595-x/663+x/714;
- if (x >= 715) {
- val += 0-x/715+x/770+x/858+x/910-x/935-x/1001-x/1105+x/1122+x/1155+x/1190
- -x/1309+x/1326+x/1365+x/1430-x/1547+x/1785+x/1870+x/2002+x/2145
- +x/2210-x/2310-x/2431+x/2618-x/2730+x/2805+x/3003+x/3094+x/3315
- -x/3570+x/3927-x/4290+x/4641+x/4862+x/5005-x/5610-x/6006+x/6545
- -x/6630+x/7293+x/7735-x/7854;
- if (x >= 9282)
- val += 0-x/9282-x/10010+x/12155-x/13090-x/14586-x/15015-x/15470+x/17017
- -x/19635-x/23205-x/24310+x/30030-x/34034-x/36465+x/39270+x/46410
- -x/51051+x/72930-x/85085+x/102102+x/170170+x/255255-x/510510;
+/* static uint32_t sprime[] = {0,2, 3, 5, 7, 11, 13, 17, 19, 23}; */
+/* static uint32_t sprimorial[] = {1,2,6,30,210,2310,30030,510510}; */
+/* static uint32_t stotient[] = {1,1,2, 8, 48, 480, 5760, 92160}; */
+static const uint16_t _s0[ 1] = {0};
+static const uint16_t _s1[ 2] = {0,1};
+static const uint16_t _s2[ 6] = {0,1,1,1,1,2};
+static const uint16_t _s3[30] = {0,1,1,1,1,1,1,2,2,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,7,7,7,7,8};
+static uint16_t _s4[210];
+static uint16_t _s5[2310];
+static uint16_t _s6[30030];
+static const uint16_t* sphicache[7] = { _s0,_s1,_s2,_s3,_s4,_s5,_s6 };
+static int sphi_init = 0;
+
+static UV tablephi(UV x, uint32_t a) {
+ switch (a) {
+ case 0: return x;
+ case 1: return x-x/2;
+ case 2: return x-x/2-x/3+x/6;
+ case 3: return (x/ 30U) * 8U + sphicache[3][x % 30U];
+ case 4: return (x/ 210U) * 48U + sphicache[4][x % 210U];
+ case 5: return (x/ 2310U) * 480U + sphicache[5][x % 2310U];
+ case 6: return (x/ 30030U) * 5760U + sphicache[6][x % 30030U];
+ default: {
+ UV xp = x / 17U;
+ return ((x /30030U) * 5760U + sphicache[6][x % 30030U]) -
+ ((xp/30030U) * 5760U + sphicache[6][xp % 30030U]);
+ }
}
- return val;
}
-
-static uint32_t mapes7_32(uint32_t x) {
- uint32_t val = x-x/2-x/3-x/5+x/6-x/7+x/10-x/11-x/13+x/14+x/15-x/17+x/21
- +x/22+x/26-x/30+x/33+x/34+x/35+x/39-x/42+x/51+x/55+x/65-x/66-x/70
- +x/77-x/78+x/85+x/91-x/102-x/105-x/110+x/119-x/130+x/143-x/154
- -x/165-x/170-x/182+x/187-x/195+x/210+x/221-x/231-x/238-x/255-x/273
- -x/286+x/330-x/357-x/374-x/385+x/390-x/429-x/442-x/455+x/462+x/510
- +x/546-x/561-x/595-x/663+x/714;
- if (x >= 715) {
- val += 0-x/715+x/770+x/858+x/910-x/935-x/1001-x/1105+x/1122+x/1155+x/1190
- -x/1309+x/1326+x/1365+x/1430-x/1547+x/1785+x/1870+x/2002+x/2145
- +x/2210-x/2310-x/2431+x/2618-x/2730+x/2805+x/3003+x/3094+x/3315
- -x/3570+x/3927-x/4290+x/4641+x/4862+x/5005-x/5610-x/6006+x/6545
- -x/6630+x/7293+x/7735-x/7854;
- if (x >= 9282)
- val += 0-x/9282-x/10010+x/12155-x/13090-x/14586-x/15015-x/15470+x/17017
- -x/19635-x/23205-x/24310+x/30030-x/34034-x/36465+x/39270+x/46410
- -x/51051+x/72930-x/85085+x/102102+x/170170+x/255255-x/510510;
+static void phitableinit(void) {
+ if (sphi_init == 0) {
+ int x;
+ for (x = 0; x < 210; x++)
+ _s4[x] = ((x/ 30)* 8+_s3[x% 30])-(((x/ 7)/ 30)* 8+_s3[(x/ 7)% 30]);
+ for (x = 0; x < 2310; x++)
+ _s5[x] = ((x/ 210)* 48+_s4[x% 210])-(((x/11)/ 210)* 48+_s4[(x/11)% 210]);
+ for (x = 0; x < 30030; x++)
+ _s6[x] = ((x/2310)*480+_s5[x%2310])-(((x/13)/2310)*480+_s5[(x/13)%2310]);
+ sphi_init = 1;
}
- return val;
}
+
/* Max memory = 2*A*X bytes, e.g. 2*400*24000 = 18.3 MB */
#define PHICACHEA 257
#define PHICACHEX 32769
@@ -348,6 +330,7 @@ static void phicache_init(cache_t* cache) {
cache->val[a] = 0;
cache->max[a] = 0;
}
+ phitableinit();
}
static void phicache_free(cache_t* cache) {
int a;
@@ -385,13 +368,13 @@ static IV _phi3(UV x, UV a, int sign, const uint32_t* const primes, const uint32
{
IV sum;
- if (a < 3)
- return sign * ((a==0) ? x : (a==1) ? x-x/2 : x-x/2-x/3+x/6);
- else if (PHI_CACHE_POPULATED(x, a))
+ if (PHI_CACHE_POPULATED(x, a))
return sign * cache->val[a][x];
else if (x < primes[a+1])
sum = sign;
- else if (1 && x <= primes[lastidx] && x < primes[a]*primes[a])
+ else if (a <= 7)
+ sum = sign * tablephi(x,a);
+ else if (x <= primes[lastidx] && x < primes[a]*primes[a])
sum = sign * (bs_prime_count(x, primes, lastidx) - a + 1);
else {
UV a2;
@@ -401,18 +384,13 @@ static IV _phi3(UV x, UV a, int sign, const uint32_t* const primes, const uint32
sum = sign * (x - a + iters);
for (a2 = 1; a2 <= iters; a2++)
sum += _phi3( FAST_DIV(x, primes[a2]), a2-1, -sign, primes, lastidx, cache);
- } else if (a >= 7) {
+ } else {
if (PHI_CACHE_POPULATED(x, 7))
sum = sign * cache->val[7][x];
else
- sum = sign * mapes7(x);
+ sum = sign * tablephi(x, 7);
for (a2 = 8; a2 <= a; a2++)
sum += _phi3( FAST_DIV(x,primes[a2]), a2-1, -sign, primes, lastidx, cache);
- } else {
- if (PHI_CACHE_POPULATED(x, a))
- sum = sign * cache->val[a][x];
- else
- sum = sign * mapes(x, a);
}
}
if (a < PHICACHEA && x < PHICACHEX)
@@ -585,8 +563,9 @@ static UV phi(UV x, UV a)
vc_t* arr;
cache_t pcache; /* Cache for recursive phi */
+ phitableinit();
if (a == 1) return ((x+1)/2);
- if (a <= 7) return mapes(x, a);
+ if (a <= 7) return tablephi(x, a);
lastidx = a+1;
primes = generate_small_primes(lastidx);
@@ -659,7 +638,7 @@ static UV phi(UV x, UV a)
#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 );
+ sum += arr[i].c * tablephi( arr[i].v, 7 );
vcarray_destroy(&a1);
Safefree(primes);
return (UV) sum;
@@ -902,14 +881,13 @@ static const unsigned char primes_small[] =
UV _XS_legendre_phi(UV x, UV a) {
/* For small values, calculate directly */
- if (a < 3) return (a == 0) ? x : (a == 1) ? x-x/2 : x-x/2-x/3+x/6;
- if (a <= 7) return mapes(x, a);
+ if (a <= 7) return tablephi(x, a);
/* For large values, do our non-recursive phi */
if (a > NPRIMES_SMALL) return phi(x,a);
/* Otherwise, recurse */
{
UV i;
- UV sum = mapes7(x);
+ UV sum = tablephi(x, 7);
for (i = 8; i <= a; i++) {
uint32_t p = primes_small[i];
UV xp = x/p;
diff --git a/lmo.c b/lmo.c
index e903f24..e85a08c 100644
--- a/lmo.c
+++ b/lmo.c
@@ -453,7 +453,7 @@ UV _XS_LMO_pi(UV n)
uint16 *factor_table;
sieve_t ss;
- const uint32 c = 7; /* We can use out Mapes function for c <= 7 */
+ const uint32 c = 7; /* We can use our Mapes function for c <= 7 */
/* For "small" n, use our table+segment sieve. */
if (n < SIEVE_LIMIT || n < 10000) return _XS_prime_count(2, n);
--
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