diff --git a/gitignore b/gitignore new file mode 100644 index 0000000..86a61fb --- /dev/null +++ b/gitignore @@ -0,0 +1,6 @@ +* +!**/*.c +!**/*.h +# !**/**/*.h +# !**/**/*.c +!**/Makefile diff --git a/primes/Makefile b/primes/Makefile new file mode 100644 index 0000000..d8e1eb9 --- /dev/null +++ b/primes/Makefile @@ -0,0 +1,8 @@ +CFALGSD=-Wall -Wextra -g -ggdb +CFALGSF=-Wall -Wextra -Ofast -march=native -mtune=native -funroll-all-loops +CC = gcc + +fast: src/sive.c + $(CC) $(CFALGSF) src/sive.c -o fastOut +debug: src/sive.c + $(CC) $(CFALGSD) src/sive.c -o debugOut diff --git a/primes/src/sive.c b/primes/src/sive.c new file mode 100644 index 0000000..9278af3 --- /dev/null +++ b/primes/src/sive.c @@ -0,0 +1,189 @@ +#include +#include +#include +#include +#include + +typedef uint64_t WORD; +typedef uint8_t BOOL; +const BOOL false = 0; +const BOOL true = ~false; + +// const WORD limit = 542; +// const WORD limit = 1e2; +const WORD limit = 1e9; +const WORD nr = 78498; + +// Sive o Atkins hit wheel +const WORD s[] = {1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59}; +const WORD wheelLen = 16; + +int main() { + BOOL *isPrime = calloc(limit, sizeof(WORD)); + assert(isPrime != NULL); + + /* + // Initialize the sieve with enough wheels to include limit: +for n ← 60 × w + x where w ∈ {0,1,...,limit ÷ 60}, x ∈ s: + is_prime(n) ← false + */ + for (WORD w = 0; w < limit / 60; w++) { + for (WORD i = 0; i < wheelLen; i++) { + WORD n = 60 * w + s[i]; + if (n >= limit) + break; + isPrime[n] = false; + } + } + + /* + * // Put in candidate primes: +// integers which have an odd number of +// representations by certain quadratic forms. +// Algorithm step 3.1: +for n ≤ limit, n ← 4x²+y² where x ∈ {1,2,...} and y ∈ {1,3,...} // all x's odd +y's if n mod 60 ∈ {1,13,17,29,37,41,49,53}: is_prime(n) ← ¬is_prime(n) // +toggle state + */ + WORD mask1[] = {1, 13, 17, 29, 37, 41, 49, 53}; + WORD maskLen1 = 8; + for (WORD x = 1; 4 * x * x < limit; x += 1) { + for (WORD y = 1; 1; y += 2) { + WORD n = 4 * x * x + y * y; + if (n >= limit) { + break; + } + WORD mod = n % 60; + for (WORD i = 0; i < maskLen1; i++) { + if (mod == mask1[i]) { + isPrime[n] = ~isPrime[n]; + // printf("3.1 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not"); + break; + } + } + } + } + + /* + // Algorithm step 3.2: +for n ≤ limit, n ← 3x²+y² where x ∈ {1,3,...} and y ∈ {2,4,...} // only odd x's + if n mod 60 ∈ {7,19,31,43}: // and even y's + is_prime(n) ← ¬is_prime(n) // toggle state + */ + WORD mask2[] = {7, 19, 31, 43}; + WORD maskLen2 = 4; + for (WORD x = 1; x * x * 3 < limit; x += 2) { + for (WORD y = 2; y * y < limit; y += 2) { + + WORD n = 3 * x * x + y * y; + if (n >= limit) { + break; + } + WORD mod = n % 60; + for (WORD i = 0; i < maskLen2; i++) { + if (mod == mask2[i]) { + isPrime[n] = ~isPrime[n]; + // printf("3.2 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not"); + break; + } + } + } + } + + /* + // Algorithm step 3.3: + for n ≤ limit, n ← 3x²-y² where x ∈ {2,3,...} and y ∈ {x-1,x-3,...,1} //all + even/odd if n mod 60 ∈ {11,23,47,59}: // + odd/even combos is_prime(n) ← ¬is_prime(n) // toggle state + */ + WORD mask3[] = {11, 23, 47, 59}; + WORD maskLen3 = 4; + for (WORD x = 2; (3 * x * x) - ((x - 1) * (x - 1)) < limit; x += 1) { + WORD y = x - 1; + while (y >= 1) { + // printf("x: %lu, y:%lu\n", x, y); + // sleep(1); + + WORD n = 3 * x * x - y * y; + + if (n >= limit) { + // breakOuter = true; + if (y < 2) { + break; + } + y -= 2; + continue; + } + WORD mod = n % 60; + for (WORD i = 0; i < maskLen3; i++) { + if (mod == mask3[i]) { + isPrime[n] = ~isPrime[n]; + // printf("3.3 %lu is %s prime\n", n, (isPrime[n]) ? "" : "not"); + break; + } + } + y = (y > 2) ? y - 2 : 0; + } + } + + /* + * +// Eliminate composites by sieving, only for those occurrences on the wheel: +for n² ≤ limit, n ← 60 × w + x where w ∈ {0,1,...}, x ∈ s, n ≥ 7: + if is_prime(n): + // n is prime, omit multiples of its square; this is sufficient + // because square-free composites can't get on this list + for c ≤ limit, c ← n² × (60 × w + x) where w ∈ {0,1,...}, x ∈ s: + is_prime(c) ← false + + */ + + for (WORD w = 0; w / 60 < limit; w++) { + WORD n = ~0; + for (WORD i = 0; i < wheelLen; i++) { + n = 60 * w + s[i]; + if (n * n >= limit) { + break; + } + // or break;? + + if (n < limit && isPrime[n]) { + for (WORD ww = 0; n * n * (60 * ww) < limit; ww++) { + WORD c; + for (WORD ii = 0; ii < wheelLen; ii++) { + c = n * n * (60 * ww + s[ii]); + if (c < limit) { + isPrime[c] = false; + // printf("elim: %lu is %s prime\n", c, (isPrime[c]) ? "" : + // "not"); + } + } + } + } + } + } + + // isPrime[2] = true; + // isPrime[3] = true; + // isPrime[5] = true; + printf("2 3 5"); + WORD sum = 3; + for (WORD w = 0; w / 60 < limit; w++) { + WORD n = ~0; + for (WORD i = 0; i < wheelLen; i++) { + n = 60 * w + s[i]; + if (n < 7) + continue; + if (n >= limit) + break; + if (isPrime[n]) { + // printf(" %lu ", n); + sum += 1; + } + } + } + printf("\n"); + printf("%lu\n", sum); + + return 0; +} diff --git a/ub10/a.out b/ub10/a.out new file mode 100755 index 0000000..72e439d Binary files /dev/null and b/ub10/a.out differ