/* prime.c -- programming with prime numbers
* compile as: gcc -lgmp -o prime prime.c
* run as: ./prime -- takes no input */
#include <stdio.h>
#include <malloc.h>
#include <string.h>
#include <gmp.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;
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);
}
/* miller-rabin primality checker */
int is_even(mpz_t n)
{
mpz_t t;
mpz_init(t);
mpz_mod_ui(t, n, 2);
int r = mpz_cmp_ui(t, 0) ? 0 : 1;
mpz_clear(t);
return r;
}
void power_mod(mpz_t r, mpz_t b, mpz_t e, mpz_t m)
{
mpz_set_ui(r, 1);
while (mpz_cmp_ui(e, 0) > 0)
{
if (! is_even(e))
{
mpz_mul(r, r, b);
mpz_mod(r, r, m);
}
mpz_fdiv_q_ui(e, e, 2);
mpz_mul(b, b, b);
mpz_mod(b, b, m);
}
}
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 r = 0;
int s = 0;
while (is_even(d))
{
mpz_divexact_ui(d, d, 2);
s += 1;
}
power_mod(t, a, d, n);
if (mpz_cmp_ui(t, 1) == 0)
{
return 1;
}
while (r < s)
{
power_mod(t, a, d, n);
if (mpz_cmp(t, n1) == 0)
{
return 1;
}
r += 1;
mpz_mul_ui(d, d, 2);
}
return 0;
mpz_clears(d, n1, t, NULL);
}
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, 1234567890);
is_seeded = 1;
}
mpz_t a, n3;
mpz_inits(a, n3, NULL);
mpz_sub_ui(n3, n, 3);
int k = 25;
while (k > 0)
{
mpz_urandomm(a, gmpRandState, n3);
mpz_add_ui(a, a, 2);
if (! is_spsp(n, a))
{
return 0;
}
k -= 1;
}
return 1;
mpz_clears(a, n3, NULL);
}
/* integer factorization by trial division */
void td_factors(List **fs, mpz_t n)
{
mpz_t d, d2, q, r;
mpz_inits(d, d2, q, r, NULL);
mpz_set_ui(d, 2);
while (is_even(n))
{
mpz_t *f = malloc(sizeof(*f));
mpz_init_set(*f, d);
*fs = insert(*f, *fs);
mpz_divexact(n, n, d);
}
mpz_set_ui(d, 3);
mpz_set_ui(d2, 9);
while (mpz_cmp(n, d2) >= 0)
{
mpz_fdiv_qr(q, r, n, d);
if (mpz_cmp_ui(r, 0) > 0)
{
mpz_add_ui(d, d, 2);
mpz_mul(d2, d, d);
}
else
{
mpz_t *f = malloc(sizeof(*f));
mpz_init_set(*f, d);
*fs = insert(*f, *fs);
mpz_divexact(n, n, d);
}
}
if (mpz_cmp_ui(n, 1) != 0)
{
*fs = reverse(insert(n, *fs));
}
mpz_clears(d, d2, q, r, NULL);
}
/* integer factorization by pollard rho */
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;
}
}
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_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);
}
/* demonstration */
int main(int argc, char *argv[])
{
mpz_t n;
mpz_init(n);
List *ps = NULL;
List *fs = NULL;
printf("Primes less than a hundred:");
ps = primes(100);
while (ps != NULL)
{
printf(" %ld", (long) ps->data);
ps = ps->next;
}
printf("\n");
printf("Count of primes less than a million: %d\n",
length(primes(1000000)));
mpz_set_str(n, "600851475143", 10);
printf("Is_even(600851475143) is %d\n", is_even(n));
mpz_t b, e, m, t;
mpz_init_set_ui(b, 437);
mpz_init_set_ui(e, 13);
mpz_init_set_ui(m, 1741);
mpz_init(t);
power_mod(t, b, e, m);
printf("Power_mod of 437^13 (mod 1741) is %s\n", mpz_get_str(NULL, 10, t));
mpz_set_str(n, "600851475143", 10); /* composite */
printf("By 25 iterations of Miller-Random, 600851475143 is %s\n",
is_prime(n) ? "prime" : "composite");
mpz_set_str(n, "2305843009213693951", 10); /* prime */
printf("By 25 iterations of Miller-Random, 2305843009213693951 is %s\n",
is_prime(n) ? "prime" : "composite");
printf("Factors of 600851475143 by trial division:");
mpz_set_str(n, "600851475143", 10);
td_factors(&fs, n);
while (fs != NULL) {
printf(" %s", mpz_get_str(NULL, 10, fs->data));
fs = fs->next;
}
printf("\n");
printf("Factors of 600851475143 by Pollard rho:");
mpz_set_str(n, "600851475143", 10);
rho_factors(&fs, n);
while (fs != NULL) {
printf(" %s", mpz_get_str(NULL, 10, fs->data));
fs = fs->next;
}
printf("\n");
}