#include <stddef.h>
#include <stdlib.h>
#include <inttypes.h>
#include <stdio.h>
#include <stdbool.h>
#include <string.h>

bool bit_read(uint64_t const* const buf, size_t const idx)
{
    size_t const offset = idx >> 6; // log2(64)
    uint64_t const shift = idx & 63; // 64 - 1

    return buf[offset] >> shift & 0x01;
}

void bit_set(uint64_t* const buf, size_t const idx)
{
    size_t const offset = idx >> 6; // log2(64)
    uint64_t const shift = idx & 63; // 64 - 1

    buf[offset] = buf[offset] ^ (1ull << shift);
}

void bit_unset(uint64_t* const buf, size_t const idx)
{
    size_t const offset = idx >> 6; // log2(64)
    uint64_t const shift = idx & 63; // 64 - 1

    buf[offset] = buf[offset] & ~(1ull << shift);
}

/**
 * in: bit-buffer with len many bits.
 * out: bit-buffer is sieved with prime bit-indexes 2n+1 set.
 */
void eratosthenes(uint64_t* buf, size_t len)
{
    memset(buf, 0xAA, len >> 3); // log2(8) from bits->bytes
    bit_unset(buf, 0);
    bit_unset(buf, 1);
    bit_unset(buf, 2);

    for(size_t i = 3; i < len; i += 2)
    {
        if(bit_read(buf, i))
        {
            for(size_t j = i * 2; j < len; j += i)
            {
                bit_unset(buf, j);
            }
        }
    }
}

typedef struct
{
    size_t* list;
    size_t len;
    size_t cap;
} List;

void list_init(List* const list)
{
    list->list = malloc(8 * sizeof(size_t));
    list->len = 0;
    list->cap = 8;
}

void list_append(List* const list, size_t const item)
{
    if(list->len + 1 > list->cap)
    {
        list->cap = list->cap + (list->cap >> 2); // growth by 1.5
        list->list = realloc(list->list, list->cap * sizeof(size_t));
    }

    list->list[list->len] = item;
    list->len++;
}

void list_print(List const* const list)
{
    char const* comma = "";
    for(size_t i = 0; i < list->len; i++)
    {
        printf("%s%zu", comma, list->list[i]);
        comma = ", ";
    }
}

void factorize(size_t const num, List* const factors, List const* const primes)
{
    size_t prime_idx = 0;
    size_t n = num;

    while(n > 1)
    {
        size_t const prime = primes->list[prime_idx];
        if(n % prime == 0)
        {
            n = n / prime;
            list_append(factors, prime);
        }
        else
        {
            prime_idx++;
        }

        if(prime_idx > primes->len)
        {
            printf("out of primes\n");
            abort();
        }
    }
}

int main()
{
    size_t const NUMBER = 1000;
    uint64_t* buf = malloc(NUMBER >> 3);
    eratosthenes(buf, NUMBER);

    List factors; list_init(&factors);
    List primes;  list_init(&primes);

    for(size_t i = 0; i < NUMBER; i++)
    {
        if(bit_read(buf, i))
        {
            printf("prime: %zu", i);
            list_append(&primes, i);

            size_t const square = i * i;
            printf(", square: %zu", square);

            factors.len = 0;
            factorize(square, &factors, &primes);
            printf(", square factors: ");
            list_print(&factors);
            printf("\n");
        }
    }

    return 0;
}