primes.cpp

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

#define gbit(a, i) (((a)[(i) / 64] >> ((i) % 64)) & 1)
#define tbit(a, i) (a)[(i) / 64] ^= UINT64_C(1) << ((i) % 64)
#define cbit(a, i) (a)[(i) / 64] &= ~(UINT64_C(1) << ((i) % 64))

//extern "C"
uint64_t* atkin(uint64_t max, size_t* lenx)
{
    size_t len;
    if (max <= 5) {
        if (max < 2) {
            len = 0;
            return (uint64_t *)malloc(0);
        } else if (max < 3) {
            len = 1;
            uint64_t* primes = (uint64_t *)malloc(1 * sizeof(uint64_t));
            primes[0] = 2;
            return primes;
        } else if (max < 5) {
            len = 2;
            uint64_t* primes = (uint64_t *)malloc(2 * sizeof(uint64_t));
            primes[0] = 2;
            primes[1] = 3;
            return primes;
        } else {
            len = 3;
            uint64_t* primes = (uint64_t *)malloc(3 * sizeof(uint64_t));
            primes[0] = 2;
            primes[1] = 3;
            primes[2] = 5;
            return primes;
        }
    }
    uint64_t* sieve =
        (uint64_t *)calloc(((max + 1) + 63) / 64, sizeof(*sieve));

    uint64_t limit = sqrt(max + 1);
    for (uint64_t x = limit; x--;)
    {
        for (uint64_t y = limit; y--;)
        {
            uint64_t index = 4 * x * x + y * y;
            if (index <= max) {
                int mod = index % 60;
                if (mod == 1 || mod == 13 || mod == 17 || mod == 29 ||
                    mod == 37 || mod == 41 || mod == 49 || mod == 53) {
                    tbit(sieve, index);
                }
            }

            index = 3 * x * x + y * y;
            if (index <= max) {
                int mod = index % 60;
                if (mod == 7 || mod == 19 || mod == 31 || mod == 43) {
                    tbit(sieve, index);
                }
            }

            if (x <= y)
                continue;

            index = 3 * x * x - y * y;
            if (index <= max) {
                int mod = index % 60;
                if (mod == 11 || mod == 23 || mod == 47 || mod == 59) {
                    tbit(sieve, index);
                }
            }
        }
    }

    len = max / 6;
    uint64_t* primes = (uint64_t *)malloc(len * sizeof(uint64_t));
    primes[0] = 2;
    primes[1] = 3;
    primes[2] = 5;
    len = 3;
    for (uint64_t i = 7; i <= limit; i++)
    {
        if (gbit(sieve, i)) {
            uint64_t val = i * i;
            for (uint64_t k = val; k <= max; k += val)
                cbit(sieve, k);
        }
    }

    for (uint64_t i = 7; i <= max; i++)
    {
        if (gbit(sieve, i))
            primes[len++] = i;
    }

    free(sieve);
    *lenx = len;
    return primes;
}