cprime/src/sieve.c

#include "sieve.h"

#include <stdlib.h>
#include <math.h>

static int isprime(sieve *s, int64_t num);
static void setnotprime(sieve *s, int64_t num);
static void setnotprimes(sieve *s, int64_t min, int64_t max);

void sieve_init(sieve *s, int64_t max)
{
	s->nbytes = (max / 16) + 1;
	s->bytes = calloc(s->nbytes, 1);

	s->maxNum = max;

	s->primes = NULL;
	s->nprimes = 0;
}

void sieve_free(sieve *s)
{
	free(s->bytes);
	if (s->primes != NULL)
		free(s->primes);
}

void sieve_generate(sieve *s)
{
	// Generate primes
	int64_t max = s->maxNum;
	int64_t sqrtnum = sqrt(max);
	for (int64_t i = 3; i <= sqrtnum; i += 2)
	{
		if (isprime(s, i)) {
			setnotprimes(s, i, max);
		}
	}

	// Generate primes array
	s->nprimes = max / ((int)log(max) - 1) + 1;
	s->primes = calloc(s->nprimes, sizeof(int64_t));

	int i = 0;
	s->primes[i++] = 2;
	for (int64_t p = 3; p < max; p += 2) {
		if (isprime(s, p))
			s->primes[i++] = p;
	}
}

void sieve_factor(sieve *s, int64_t num, int64_t *buf)
{
	int j = 0;
	int64_t sqrtnum = sqrt(num);
	for (int i = 0; i < s->nprimes; ++i)
	{
		int64_t p = s->primes[i];
		if (p == 0)
			break;
		if (p > sqrtnum)
			break;

		if (num % p == 0)
		{
			do
			{
				buf[j++] = p;
				num /= p;
			} while (num % p == 0);
			sqrtnum = sqrt(num);
		}
	}

	if (num != 1 || j == 0)
		buf[j++] = num;
	buf[j++] = 0;
}

/*
 * See whether a number is prime or not
 */
static int isprime(sieve *s, int64_t num)
{
	if (num == 2)
		return 1;
	if ((num & 1) == 0)
		return 0;

	unsigned char c = s->bytes[num >> 4];
	return (c & (1 << ((num % 16) >> 1))) == 0;
}

/*
 * Set number to be not prime (aka set its bit to 1)
 */
static void setnotprime(sieve *s, int64_t num)
{
	if ((num & 1) == 0)
		return;

	int bi = num >> 4;
	s->bytes[bi] = s->bytes[bi] | (1 << ((num % 16) >> 1));
}

static void setnotprimes(sieve *s, int64_t min, int64_t max)
{
	int64_t j = min;
	min *= 2;
	int64_t mul = 3;
	while (min <= max) {
		setnotprime(s, min);
		min = j * mul;
		mul += 2;
	}
}