`/* Programming With Prime Numbers */`

/* http://programmingpraxis.files.wordpress.com/

2012/09/programmingwithprimenumbers.pdf */

```
```#include <stdlib.h>

#include <stdio.h>

#include <string.h>

#include <gmp.h>

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

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 *insert_in_order(void *x, List *xs)

{

if (xs == NULL || mpz_cmp(x, xs->data) < 0)

{

return insert(x, xs);

}

else

{

List *head = xs;

while (xs->next != NULL &&

mpz_cmp(x, xs->next->data) > 0)

{

xs = xs->next;

}

xs->next = insert(x, xs->next);

return head;

}

}

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;

}

List *sieve(long n)

{

char b[(n+1)/8+1];

int i, p;

List *ps = NULL;

memset(b, 255, sizeof(b));

for (p=2; p<=n; p++)

{

if (ISBITSET(b,p))

{

ps = insert((void *) p, ps);

for (i=p; i<=n; i+=p)

{

CLEARBIT(b,i);

}

}

}

return reverse(ps);

}

List *primes(long n)

{

int m = (n-1) / 2;

char b[m/8+1];

int i = 0;

int p = 3;

List *ps = NULL;

int j;

ps = insert((void *) 2, ps);

memset(b, 255, 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);

}

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;

}

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);

}

int is_spsp(mpz_t n, mpz_t a)

{

mpz_t d, n1, t;

mpz_inits(d, n1, t, NULL);

mpz_sub_ui(n1, n, 1);

mpz_set(d, n1);

int s = 0;

while (mpz_even_p(d))

{

mpz_divexact_ui(d, d, 2);

s += 1;

}

mpz_powm(t, a, d, n);

if (mpz_cmp_ui(t, 1) == 0 ||

mpz_cmp(t, n1) == 0)

{

mpz_clears(d, n1, t, NULL);

return 1;

}

while (--s > 0)

{

mpz_mul(t, t, t);

mpz_mod(t, t, n);

if (mpz_cmp(t, n1) == 0)

{

mpz_clears(d, n1, t, NULL);

return 1;

}

}

mpz_clears(d, n1, t, NULL);

return 0;

}

int is_prime(mpz_t n)

{

static gmp_randstate_t gmpRandState;

static int is_seeded = 0;

if (! is_seeded)

{

gmp_randinit_default(gmpRandState);

gmp_randseed_ui(gmpRandState, time(NULL));

is_seeded = 1;

}

mpz_t a, n3, t;

mpz_inits(a, n3, t, NULL);

mpz_sub_ui(n3, n, 3);

int i;

int k = 25;

long unsigned ps[] = { 2, 3, 5, 7 };

if (mpz_cmp_ui(n, 2) < 0)

{

mpz_clears(a, n3, t, NULL);

return 0;

}

for (i = 0; i < sizeof(ps) /

sizeof(long unsigned); i++)

{

mpz_mod_ui(t, n, ps[i]);

if (mpz_cmp_ui(t, 0) == 0)

{

mpz_clears(a, n3, t, NULL);

return mpz_cmp_ui(n, ps[i]) == 0;

}

}

while (k > 0)

{

mpz_urandomm(a, gmpRandState, n3);

mpz_add_ui(a, a, 2);

if (! is_spsp(n, a))

{

mpz_clears(a, n3, t, NULL);

return 0;

}

k -= 1;

}

mpz_clears(a, n3, t, NULL);

return 1;

}

void rho_factor(mpz_t f, mpz_t n, long long unsigned c)

{

mpz_t t, h, d, r;

mpz_init_set_ui(t, 2);

mpz_init_set_ui(h, 2);

mpz_init_set_ui(d, 1);

mpz_init_set_ui(r, 0);

while (mpz_cmp_si(d, 1) == 0)

{

mpz_mul(t, t, t);

mpz_add_ui(t, t, c);

mpz_mod(t, t, n);

mpz_mul(h, h, h);

mpz_add_ui(h, h, c);

mpz_mod(h, h, n);

mpz_mul(h, h, h);

mpz_add_ui(h, h, c);

mpz_mod(h, h, n);

mpz_sub(r, t, h);

mpz_gcd(d, r, n);

}

if (mpz_cmp(d, n) == 0) /* cycle */

{

rho_factor(f, n, c+1);

}

else if (mpz_probab_prime_p(d, 25)) /* success */

{

mpz_set(f, d);

}

else /* found composite factor */

{

rho_factor(f, d, c+1);

}

mpz_clears(t, h, d, r, NULL);

}

void rho_factors(List **fs, mpz_t n)

{

while (mpz_even_p(n))

{

mpz_t *f = malloc(sizeof(*f));

mpz_init_set_ui(*f, 2);

*fs = insert(*f, *fs);

mpz_divexact_ui(n, n, 2);

}

if (mpz_cmp_ui(n, 1) == 0) return;

while (! (mpz_probab_prime_p(n, 25)))

{

mpz_t *f = malloc(sizeof(*f));

mpz_init_set_ui(*f, 0);

rho_factor(*f, n, 1);

*fs = insert_in_order(*f, *fs);

mpz_divexact(n, n, *f);

}

*fs = insert_in_order(n, *fs);

}

int main(int argc, char *argv[])

{

mpz_t n;

mpz_init(n);

List *ps = NULL;

List *fs = NULL;

ps = sieve(100);

while (ps != NULL)

{

printf("%d%s", (int) ps->data,

(ps->next == NULL) ? "\n" : " ");

ps = ps->next;

}

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(600851475143LL));

fs = td_factors(600851475143LL);

while (fs != NULL)

{

printf("%llu%s",

(unsigned long long int) fs->data,

(fs->next == NULL) ? "\n" : " ");

fs = fs->next;

}

mpz_t a;

mpz_init(a);

mpz_set_str(n, "2047", 10);

mpz_set_str(a, "2", 10);

printf("%d\n", is_spsp(n, a));

mpz_set_str(n, "600851475143", 10);

printf("%d\n", is_prime(n));

mpz_set_str(n, "2305843009213693951", 10);

printf("%d\n", is_prime(n));

` mpz_set_str(n, "600851475143", 10);`

rho_factors(&fs, n);

while (fs != NULL) {

printf("%s%s",

mpz_get_str(NULL, 10, fs->data),

(fs->next == NULL) ? "\n" : " ");

fs = fs->next;

}

}