## Programming With Prime Numbers: Source Code In C

```/* 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;     } }```