codepad

Language:

/* prime.c -- programming with prime numbers
 * compile as: gcc -O3 prime.c -lgmp -o prime
 * run as: ./prime */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <time.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;
    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(unsigned long long int n)
{
    if (n % 2 == 0)
    {
        return n == 2;
    }

unsigned long long int d = 3;

while (d*d <= n)
    {
        if (n % d == 0)
        {
            return 0;
        }
        d += 2;
    }

return 1;
}

/* factorization by trial division */

List *td_factors(unsigned long long int n)
{
    List *fs = NULL;

while (n % 2 == 0)
    {
        fs = insert((void *) 2, fs);
        n /= 2;
    }

if (n == 1)
    {
        return reverse(fs);
    }

unsigned long long int 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);
}

/* miller-rabin primality checker */

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

/* 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 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(r);

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

/* demonstration */

int main(int argc, char *argv[])
{
    mpz_t n;
    mpz_init(n);
    List *ps = NULL;
    List *fs = NULL;

ps = primes(100); /* 2 3 5 7 11 13 17 19 23 29 31 37 41 */
    while (ps != NULL) /* 43 47 53 59 61 67 71 73 79 83 89 97 */
    {
        printf("%ld%s", (long) ps->data,
            (ps->next == NULL) ? "\n" : " ");
        ps = ps->next;
    }

printf("%d\n", length(primes(1000000))); /* 78498 */

printf("%d\n", td_prime(600851475143LL)); /* composite */

fs = td_factors(600851475143LL); /* 71 839 1471 6857 */
    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)); /* pseudo-prime */

mpz_set_str(n, "600851475143", 10); /* composite */
    printf("%d\n", is_prime(n));

mpz_set_str(n, "2305843009213693951", 10); /* prime */
    printf("%d\n", is_prime(n));

mpz_set_str(n, "600851475143", 10);
    rho_factors(&fs, n); /* 71 839 1471 6857 */
    while (fs != NULL) {
        printf("%s%s", mpz_get_str(NULL, 10, fs->data),
            (fs->next == NULL) ? "\n" : " ");
        fs = fs->next;
    }
}

Private [?]

Run code