Kimberlee Integer Model revised this on 21 Aug 2022

prime_square_sieve.c

#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; }