/* prime.c -- programming with prime numbers
 * compile as: gcc -O3 prime.c -o prime
 * run as: ./prime */
 
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
 
/* bit arrays */
 
#define ISBITSET(x, i) (( x[i>>3] & (1<<(i&7)) ) != 0)
#define SETBIT(x, i)   x[i>>3] |= (1<<(i&7));
#define CLEARBIT(x, i) x[i>>3] &= (1<<(i&7)) ^ 0xFF;
 
/* basic list library */
 
typedef struct list {
    void *data;
    struct list *next;
} List;
 
List *insert(void *data, List *next)
{
    List *new;
 
    new = malloc(sizeof(List));
    new->data = data;
    new->next = next;
    return new;
}
 
List *reverse(List *list) {
    List *new = NULL;
    List *next;
 
    while (list != NULL)
    {
        next = list->next;
        list->next = new;
        new = list;
        list = next;
    }
 
    return new;
}
 
int length(List *xs)
{
    int len = 0;
    while (xs != NULL)
    {
        len += 1;
        xs = xs->next;
    }
    return len;
}
 
/* sieve of eratosthenes */
 
List *primes(long n)
{
    int m = (n-1) / 2;
    unsigned char b[(m+7)/8];
    int i = 0;
    int p = 3;
    List *ps = NULL;
    int j;
 
    ps = insert((void *) 2, ps);
 
    memset(b, 0xFF, sizeof(b));
 
    while (p*p <= n)
    {
        if (ISBITSET(b,i))
        {
            ps = insert((void *) p, ps);
            j = (p*p - 3) / 2;
            while (j < m)
            {
                CLEARBIT(b, j);
                j += p;
            }
        }
        i += 1; p += 2;
    }
 
    while (i < m)
    {
        if (ISBITSET(b,i))
        {
            ps = insert((void *) p, ps);
        }
        i += 1; p += 2;
    }
 
    return reverse(ps);
}
 
/* test primality by trial division */
 
int td_prime(long long unsigned n)
{
    if (n % 2 == 0)
    {
        return n == 2;
    }
    
    long long unsigned d = 3;
    
    while (d*d <= n)
    {
        if (n % d == 0)
        {
            return 0;
        }
        d += 2;
    }
    
    return 1;
}
 
/* factorization by trial division */
 
List *td_factors(long long unsigned n)
{
    List *fs = NULL;
    
    while (n % 2 == 0)
    {
        fs = insert((void *) 2, fs);
        n /= 2;
    }
    
    if (n == 1)
    {
        return reverse(fs);
    }
    
    long long unsigned f = 3;
    
    while (f*f <= n)
    {
        if (n % f == 0)
        {
            fs = insert((void *) f, fs);
            n /= f;
        }
        else
        {
            f += 2;
        }
    }
    
    fs = insert((void *) n, fs);
    return reverse(fs);
}
 
/* demonstration */
 
int main(int argc, char *argv[])
{
    List *ps = NULL;
    List *fs = NULL;
 
    ps = primes(100);
    while (ps != NULL)
    {
        printf("%ld%s", (long) ps->data,
            (ps->next == NULL) ? "\n" : " ");
        ps = ps->next;
    }
 
    printf("%d\n", length(primes(1000000)));
    
    printf("%d\n", td_prime(2));
 
    fs = td_factors(600851475143);
    while (fs != NULL)
    {
        printf("%llu%s", (long long unsigned) fs->data,
            (fs->next == NULL) ? "\n" : " ");
        fs = fs->next;
    }
 
    return 0;
}

2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
78498
1
71 839 1471 6857