/* 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->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);
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)));
fs = td_factors(600851475143);
while (fs != NULL)
{
printf("%llu%s", (long long unsigned) fs
->data
, (fs->next == NULL) ? "\n" : " ");
fs = fs->next;
}
return 0;
}
LyogcHJpbWUuYyAtLSBwcm9ncmFtbWluZyB3aXRoIHByaW1lIG51bWJlcnMKICogY29tcGlsZSBhczogZ2NjIC1PMyBwcmltZS5jIC1vIHByaW1lCiAqIHJ1biBhczogLi9wcmltZSAqLwogCiNpbmNsdWRlIDxzdGRpby5oPgojaW5jbHVkZSA8c3RkbGliLmg+CiNpbmNsdWRlIDxzdGRpbnQuaD4KI2luY2x1ZGUgPHN0cmluZy5oPgogCi8qIGJpdCBhcnJheXMgKi8KIAojZGVmaW5lIElTQklUU0VUKHgsIGkpICgoIHhbaT4+M10gJiAoMTw8KGkmNykpICkgIT0gMCkKI2RlZmluZSBTRVRCSVQoeCwgaSkgICB4W2k+PjNdIHw9ICgxPDwoaSY3KSk7CiNkZWZpbmUgQ0xFQVJCSVQoeCwgaSkgeFtpPj4zXSAmPSAoMTw8KGkmNykpIF4gMHhGRjsKIAovKiBiYXNpYyBsaXN0IGxpYnJhcnkgKi8KIAp0eXBlZGVmIHN0cnVjdCBsaXN0IHsKICAgIHZvaWQgKmRhdGE7CiAgICBzdHJ1Y3QgbGlzdCAqbmV4dDsKfSBMaXN0OwogCkxpc3QgKmluc2VydCh2b2lkICpkYXRhLCBMaXN0ICpuZXh0KQp7CiAgICBMaXN0ICpuZXc7CiAKICAgIG5ldyA9IG1hbGxvYyhzaXplb2YoTGlzdCkpOwogICAgbmV3LT5kYXRhID0gZGF0YTsKICAgIG5ldy0+bmV4dCA9IG5leHQ7CiAgICByZXR1cm4gbmV3Owp9CiAKTGlzdCAqcmV2ZXJzZShMaXN0ICpsaXN0KSB7CiAgICBMaXN0ICpuZXcgPSBOVUxMOwogICAgTGlzdCAqbmV4dDsKIAogICAgd2hpbGUgKGxpc3QgIT0gTlVMTCkKICAgIHsKICAgICAgICBuZXh0ID0gbGlzdC0+bmV4dDsKICAgICAgICBsaXN0LT5uZXh0ID0gbmV3OwogICAgICAgIG5ldyA9IGxpc3Q7CiAgICAgICAgbGlzdCA9IG5leHQ7CiAgICB9CiAKICAgIHJldHVybiBuZXc7Cn0KIAppbnQgbGVuZ3RoKExpc3QgKnhzKQp7CiAgICBpbnQgbGVuID0gMDsKICAgIHdoaWxlICh4cyAhPSBOVUxMKQogICAgewogICAgICAgIGxlbiArPSAxOwogICAgICAgIHhzID0geHMtPm5leHQ7CiAgICB9CiAgICByZXR1cm4gbGVuOwp9CiAKLyogc2lldmUgb2YgZXJhdG9zdGhlbmVzICovCiAKTGlzdCAqcHJpbWVzKGxvbmcgbikKewogICAgaW50IG0gPSAobi0xKSAvIDI7CiAgICB1bnNpZ25lZCBjaGFyIGJbKG0rNykvOF07CiAgICBpbnQgaSA9IDA7CiAgICBpbnQgcCA9IDM7CiAgICBMaXN0ICpwcyA9IE5VTEw7CiAgICBpbnQgajsKIAogICAgcHMgPSBpbnNlcnQoKHZvaWQgKikgMiwgcHMpOwogCiAgICBtZW1zZXQoYiwgMHhGRiwgc2l6ZW9mKGIpKTsKIAogICAgd2hpbGUgKHAqcCA8PSBuKQogICAgewogICAgICAgIGlmIChJU0JJVFNFVChiLGkpKQogICAgICAgIHsKICAgICAgICAgICAgcHMgPSBpbnNlcnQoKHZvaWQgKikgcCwgcHMpOwogICAgICAgICAgICBqID0gKHAqcCAtIDMpIC8gMjsKICAgICAgICAgICAgd2hpbGUgKGogPCBtKQogICAgICAgICAgICB7CiAgICAgICAgICAgICAgICBDTEVBUkJJVChiLCBqKTsKICAgICAgICAgICAgICAgIGogKz0gcDsKICAgICAgICAgICAgfQogICAgICAgIH0KICAgICAgICBpICs9IDE7IHAgKz0gMjsKICAgIH0KIAogICAgd2hpbGUgKGkgPCBtKQogICAgewogICAgICAgIGlmIChJU0JJVFNFVChiLGkpKQogICAgICAgIHsKICAgICAgICAgICAgcHMgPSBpbnNlcnQoKHZvaWQgKikgcCwgcHMpOwogICAgICAgIH0KICAgICAgICBpICs9IDE7IHAgKz0gMjsKICAgIH0KIAogICAgcmV0dXJuIHJldmVyc2UocHMpOwp9CgovKiB0ZXN0IHByaW1hbGl0eSBieSB0cmlhbCBkaXZpc2lvbiAqLwoKaW50IHRkX3ByaW1lKGxvbmcgbG9uZyB1bnNpZ25lZCBuKQp7CiAgICBpZiAobiAlIDIgPT0gMCkKICAgIHsKICAgICAgICByZXR1cm4gbiA9PSAyOwogICAgfQogICAgCiAgICBsb25nIGxvbmcgdW5zaWduZWQgZCA9IDM7CiAgICAKICAgIHdoaWxlIChkKmQgPD0gbikKICAgIHsKICAgICAgICBpZiAobiAlIGQgPT0gMCkKICAgICAgICB7CiAgICAgICAgICAgIHJldHVybiAwOwogICAgICAgIH0KICAgICAgICBkICs9IDI7CiAgICB9CiAgICAKICAgIHJldHVybiAxOwp9CgovKiBmYWN0b3JpemF0aW9uIGJ5IHRyaWFsIGRpdmlzaW9uICovCgpMaXN0ICp0ZF9mYWN0b3JzKGxvbmcgbG9uZyB1bnNpZ25lZCBuKQp7CiAgICBMaXN0ICpmcyA9IE5VTEw7CiAgICAKICAgIHdoaWxlIChuICUgMiA9PSAwKQogICAgewogICAgICAgIGZzID0gaW5zZXJ0KCh2b2lkICopIDIsIGZzKTsKICAgICAgICBuIC89IDI7CiAgICB9CiAgICAKICAgIGlmIChuID09IDEpCiAgICB7CiAgICAgICAgcmV0dXJuIHJldmVyc2UoZnMpOwogICAgfQogICAgCiAgICBsb25nIGxvbmcgdW5zaWduZWQgZiA9IDM7CiAgICAKICAgIHdoaWxlIChmKmYgPD0gbikKICAgIHsKICAgICAgICBpZiAobiAlIGYgPT0gMCkKICAgICAgICB7CiAgICAgICAgICAgIGZzID0gaW5zZXJ0KCh2b2lkICopIGYsIGZzKTsKICAgICAgICAgICAgbiAvPSBmOwogICAgICAgIH0KICAgICAgICBlbHNlCiAgICAgICAgewogICAgICAgICAgICBmICs9IDI7CiAgICAgICAgfQogICAgfQogICAgCiAgICBmcyA9IGluc2VydCgodm9pZCAqKSBuLCBmcyk7CiAgICByZXR1cm4gcmV2ZXJzZShmcyk7Cn0KCi8qIGRlbW9uc3RyYXRpb24gKi8KIAppbnQgbWFpbihpbnQgYXJnYywgY2hhciAqYXJndltdKQp7CiAgICBMaXN0ICpwcyA9IE5VTEw7CiAgICBMaXN0ICpmcyA9IE5VTEw7CiAKICAgIHBzID0gcHJpbWVzKDEwMCk7CiAgICB3aGlsZSAocHMgIT0gTlVMTCkKICAgIHsKICAgICAgICBwcmludGYoIiVsZCVzIiwgKGxvbmcpIHBzLT5kYXRhLAogICAgICAgICAgICAocHMtPm5leHQgPT0gTlVMTCkgPyAiXG4iIDogIiAiKTsKICAgICAgICBwcyA9IHBzLT5uZXh0OwogICAgfQogCiAgICBwcmludGYoIiVkXG4iLCBsZW5ndGgocHJpbWVzKDEwMDAwMDApKSk7CiAgICAKICAgIHByaW50ZigiJWRcbiIsIHRkX3ByaW1lKDIpKTsKCiAgICBmcyA9IHRkX2ZhY3RvcnMoNjAwODUxNDc1MTQzKTsKICAgIHdoaWxlIChmcyAhPSBOVUxMKQogICAgewogICAgICAgIHByaW50ZigiJWxsdSVzIiwgKGxvbmcgbG9uZyB1bnNpZ25lZCkgZnMtPmRhdGEsCiAgICAgICAgICAgIChmcy0+bmV4dCA9PSBOVUxMKSA/ICJcbiIgOiAiICIpOwogICAgICAgIGZzID0gZnMtPm5leHQ7CiAgICB9CgogICAgcmV0dXJuIDA7Cn0=