 Marcus - 3 years ago 199
C Question

# Counting primes up to 18 digits optimization

I have a task at school to count number of primes up to 10^18 in under 2 minutes and with no more than 2 GB of memory to use. For the first try I've implemented a segmented sieve with the following optimizations:

• used 32 bit integers to store the primes in a segmented as bits

• didn't store the odd numbers

• segmented the sieve into parts of size sqrt(n)

• counting the composite primes when marking them so I don't have to loop the sieve again

• used dynamic memory allocation to store the first primes up to sqrt(n) (in this case I've created a Queue in C to store the primes)

The problem is that for counting the primes up to 10^9 on my computer(which has pretty decent specs) it takes 13 seconds, therefore for 10^18 would take me some days.

My question is, is there some optimization that I am missing, or is there a better and faster way to count primes up to a number? The code:

``````#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>

typedef   signed          char  int8_t;
typedef   signed     short int int16_t;
typedef   signed           int int32_t;
typedef   signed long long int int64_t;

typedef unsigned          char  uint8_t;
typedef unsigned     short int uint16_t;
typedef unsigned           int uint32_t;
typedef unsigned long long int uint64_t;

#define  SIZE 32
#define DEBUG

#define  KRED "\x1B[31m"
#define  KGRN "\x1B[32m"
#define  KYEL "\x1B[33m"
#define  KBLU "\x1B[34m"
#define  KMAG "\x1B[35m"
#define  KCYN "\x1B[36m"
#define  KWHT "\x1B[37m"
#define RESET "\033[0m"

struct node {
uint64_t     data;
struct node* next;
};

struct queue {
struct node* first;
struct node* last;
uint32_t     size;
};

typedef struct node  Node;
typedef struct queue Queue;

/* Queue model */
uint8_t enqueue(Queue* queue, int64_t value) {
Node* node = (Node*)malloc(sizeof(Node*));

if (node == NULL)
return 0;

node->data = value;
if (queue->last)
queue->last->next = node;

queue->last = node;
if (queue->first == NULL)
queue->first = queue->last;

queue->size++;
return 1;
}

uint64_t dequeue(Queue* queue) {
Node*         node = queue->first;
uint64_t save_data = node->data;

if (queue->size == 0)
return 0;

queue->first = queue->first->next;
queue->size--;
free(node);

return save_data;
}

Node* queue_peek(Queue* queue) {
return queue->first;
}

uint32_t queue_size(Queue* queue) {
return queue->size;
}

Queue* init_queue() {
Queue* queue = (Queue*)malloc(sizeof(Queue*));

queue->first = queue->last = NULL;
queue->size  = 0;

return queue;
}

/* Working with bit arrays functions */
uint8_t count_set_bits(uint64_t nbr) {
uint8_t count = 0;

while (nbr) {
count++;
nbr &= (nbr - 1);
}

return count;
}

uint8_t get_bit(uint32_t array[], uint32_t position) {
const uint64_t mask = 1U << (position % SIZE);

return array[position / SIZE] & mask ? 1 : 0;
}

void clear_bit(uint32_t array[], uint32_t position) {
const uint64_t mask = ~(1U << (position % SIZE));

}

void set_bit(uint32_t array[], uint32_t position) {
array[position / SIZE] |= (1U << (position % SIZE));
}

/* Solve the problem */
Queue* initial_sieve(uint64_t limit) {
Queue*   queue   = init_queue();
uint64_t _sqrt   = (uint64_t)sqrt(limit);
uint32_t *primes = (uint32_t*)calloc(_sqrt / SIZE + 1, sizeof(uint32_t));

set_bit(primes, 0);
// working with reversed logic, otherwise primes should all me initialized to max uiint64_t

enqueue(queue, 2);
for (uint64_t number = 3; number <= _sqrt; number += 2) {
if (get_bit(primes, number / 2) == 0) {
enqueue(queue, number);

for (uint64_t position = number * number; position <= _sqrt; position += (number * 2)) {
set_bit(primes, position / 2);
}
}
else
set_bit(primes, number / 2);
}

return queue;
}

uint64_t count_primes(uint64_t limit) {
uint64_t start, end, delta;
uint64_t non_primes_counter, initial_size;
uint32_t *current_sieve;
Queue* queue;

queue = initial_sieve(limit);
initial_size = queue->size;
start = delta = (uint64_t)sqrt(limit);
end   = 2 * start;
non_primes_counter = 0;

printf("Limits: %llu -> %llu\n", start, end);
while (start < limit) {
Node*    prime = queue->first->next; // pass 2 since only odd maps are represented in the sieve
uint64_t count = 0;

current_sieve = (uint32_t*)calloc(delta / SIZE + 1, sizeof(uint32_t));
// memset(current_sieve, 0, sizeof(uint32_t) * delta);

while (prime != NULL) {
uint64_t first_composite = start / prime->data * prime->data;

// calculate the first multiple of the given prime in the interval
if (first_composite < start)
first_composite += prime->data;
if ((first_composite & 1) == 0)
first_composite += prime->data;

// set all the composites of the current prime in the given interval
for (uint64_t number = first_composite; number <= end; number += (prime->data) * 2) {
const uint64_t position = (number - start) / 2;

if (get_bit(current_sieve, position) == 0) {
set_bit(current_sieve, position);
count++;
}
}

// free(current_sieve);
prime = prime->next;
}

non_primes_counter += count;
start += delta;
end   += delta;

if (end > limit)
end = limit;
}

uint64_t total = (limit - delta + 1) / 2 - non_primes_counter;

printf("%sTotal composites and initial size: %llu %llu %s\n", KCYN, non_primes_counter, initial_size, RESET);
printf("%sTotal primes: %llu %s\n", KCYN, total, RESET);
return queue->size + (limit  - delta + 1) / 2 - non_primes_counter;
}

/* Main */
int main(int argc, char **argv) {
clock_t begin, end;
double  time;

if (argc < 2) {
printf("Invalid number of parameters\n");
printf("Program will exit now.\n");
return 0;
}

begin = clock();
printf("%sNumber of primes found up to %s%s: %s%llu.\n%s", KWHT, KCYN, argv, KYEL, count_primes(atoll(argv)), RESET);
end     = clock();
time    = (double)(end - begin) / CLOCKS_PER_SEC;

printf("%sTotal time elapsed since the starting of the program: %s%lf seconds.\n%s", KWHT, KYEL, time, RESET);
return 0;
}
``````

Thanks, Marcus DAle