ketrptr
parent
ecbbcc1f72
commit
1843315b51
|
|
@ -1,8 +0,0 @@
|
|||
CFALGSD=-Wall -Wextra -g -ggdb
|
||||
CFALGSF=-Wall -Wextra -O3 -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
|
||||
|
|
@ -1,196 +0,0 @@
|
|||
#include <assert.h>
|
||||
#include <stdint.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <time.h>
|
||||
#include <unistd.h>
|
||||
|
||||
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 = 1e8;
|
||||
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 sieve() {
|
||||
|
||||
clock_t t0 = clock();
|
||||
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);
|
||||
clock_t t1 = clock();
|
||||
clock_t delta = t1 - t0;
|
||||
printf("clocks: %ld, clocks per sec: %ld, s: %ld ", delta, CLOCKS_PER_SEC,
|
||||
delta / CLOCKS_PER_SEC);
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue