fork(1) download
copy 
  1. /* prime.c -- programming with prime numbers
  2.  * compile as: gcc -O3 prime.c -lgmp -o prime
  3.  * run as: ./prime */
  4.  
  5. #include <stdio.h>
  6. #include <stdlib.h>
  7. #include <stdint.h>
  8. #include <string.h>
  9. #include <time.h>
  10. #include <gmp.h>
  11.  
  12. /* bit arrays */
  13.  
  14. #define ISBITSET(x, i) ((x[i>>3] &  (1<<(i&7))) != 0)
  15. #define SETBIT  (x, i)   x[i>>3] |= (1<<(i&7))       ;
  16. #define CLEARBIT(x, i)   x[i>>3] &= (1<<(i&7)) ^ 0xFF;
  17.  
  18. /* basic list library */
  19.  
  20. typedef struct list {
  21.     void *data;
  22.     struct list *next;
  23. } List;
  24.  
  25. List *insert(void *data, List *next)
  26. {
  27.     List *new;
  28.  
  29.     new = malloc(sizeof(List));
  30.     new->data = data;
  31.     new->next = next;
  32.     return new;
  33. }
  34.  
  35. List *reverse(List *list) {
  36.     List *new = NULL;
  37.     List *next;
  38.  
  39.     while (list != NULL)
  40.     {
  41.         next = list->next;
  42.         list->next = new;
  43.         new = list;
  44.         list = next;
  45.     }
  46.  
  47.     return new;
  48. }
  49.  
  50. int length(List *xs)
  51. {
  52.     int len = 0;
  53.     while (xs != NULL)
  54.     {
  55.         len += 1;
  56.         xs = xs->next;
  57.     }
  58.     return len;
  59. }
  60.  
  61. /* sieve of eratosthenes */
  62.  
  63. List *primes(long n)
  64. {
  65.     int m = (n-1) / 2;
  66.     unsigned char b[(m+7)/8];
  67.     int i = 0;
  68.     int p = 3;
  69.     List *ps = NULL;
  70.     int j;
  71.  
  72.     ps = insert((void *) 2, ps);
  73.  
  74.     memset(b, 0xFF, sizeof(b));
  75.  
  76.     while (p*p <= n)
  77.     {
  78.         if (ISBITSET(b,i))
  79.         {
  80.             ps = insert((void *) p, ps);
  81.             j = (p*p - 3) / 2;
  82.             while (j < m)
  83.             {
  84.                 CLEARBIT(b, j);
  85.                 j += p;
  86.             }
  87.         }
  88.         i += 1; p += 2;
  89.     }
  90.  
  91.     while (i < m)
  92.     {
  93.         if (ISBITSET(b,i))
  94.         {
  95.             ps = insert((void *) p, ps);
  96.         }
  97.         i += 1; p += 2;
  98.     }
  99.  
  100.     return reverse(ps);
  101. }
  102.  
  103. /* test primality by trial division */
  104.  
  105. int td_prime(unsigned long long int n)
  106. {
  107.     if (n % 2 == 0)
  108.     {
  109.         return n == 2;
  110.     }
  111.  
  112.     unsigned long long int d = 3;
  113.  
  114.     while (d*d <= n)
  115.     {
  116.         if (n % d == 0)
  117.         {
  118.             return 0;
  119.         }
  120.         d += 2;
  121.     }
  122.  
  123.     return 1;
  124. }
  125.  
  126. /* factorization by trial division */
  127.  
  128. List *td_factors(unsigned long long int n)
  129. {
  130.     List *fs = NULL;
  131.  
  132.     while (n % 2 == 0)
  133.     {
  134.         fs = insert((void *) 2, fs);
  135.         n /= 2;
  136.     }
  137.  
  138.     if (n == 1)
  139.     {
  140.         return reverse(fs);
  141.     }
  142.  
  143.     unsigned long long int f = 3;
  144.  
  145.     while (f*f <= n)
  146.     {
  147.         if (n % f == 0)
  148.         {
  149.             fs = insert((void *) f, fs);
  150.             n /= f;
  151.         }
  152.         else
  153.         {
  154.             f += 2;
  155.         }
  156.     }
  157.  
  158.     fs = insert((void *) n, fs);
  159.     return reverse(fs);
  160. }
  161.  
  162. /* miller-rabin primality checker */
  163.  
  164. int is_spsp(mpz_t n, mpz_t a)
  165. {
  166.     mpz_t d, n1, t;
  167.     mpz_inits(d, n1, t, NULL);
  168.     mpz_sub_ui(n1, n, 1);
  169.     mpz_set(d, n1);
  170.     int s = 0;
  171.  
  172.     while (mpz_even_p(d))
  173.     {
  174.         mpz_divexact_ui(d, d, 2);
  175.         s += 1;
  176.     }
  177.  
  178.     mpz_powm(t, a, d, n);
  179.     if (mpz_cmp_ui(t, 1) == 0 || mpz_cmp(t, n1) == 0)
  180.     {
  181.         mpz_clears(d, n1, t, NULL);
  182.         return 1;
  183.     }
  184.  
  185.     while (--s > 0)
  186.     {
  187.         mpz_mul(t, t, t);
  188.         mpz_mod(t, t, n);
  189.         if (mpz_cmp(t, n1) == 0)
  190.         {
  191.             mpz_clears(d, n1, t, NULL);
  192.             return 1;
  193.         }
  194.     }
  195.  
  196.     mpz_clears(d, n1, t, NULL);
  197.     return 0;
  198. }
  199.  
  200. int is_prime(mpz_t n)
  201. {
  202.     static gmp_randstate_t gmpRandState;
  203.     static int is_seeded = 0;
  204.  
  205.     if (! is_seeded)
  206.     {
  207.         gmp_randinit_default(gmpRandState);
  208.         gmp_randseed_ui(gmpRandState, time(NULL));
  209.         is_seeded = 1;
  210.     }
  211.  
  212.     mpz_t a, n3, t;
  213.     mpz_inits(a, n3, t, NULL);
  214.     mpz_sub_ui(n3, n, 3);
  215.     int i;
  216.     int k = 25;
  217.     long unsigned ps[] = { 2, 3, 5, 7 };
  218.  
  219.     if (mpz_cmp_ui(n, 2) < 0)
  220.     {
  221.         mpz_clears(a, n3, t, NULL);
  222.         return 0;
  223.     }
  224.  
  225.     for (i = 0; i < sizeof(ps) / sizeof(long unsigned); i++)
  226.     {
  227.         mpz_mod_ui(t, n, ps[i]);
  228.         if (mpz_cmp_ui(t, 0) == 0)
  229.         {
  230.             mpz_clears(a, n3, t, NULL);
  231.             return mpz_cmp_ui(n, ps[i]) == 0;
  232.         }
  233.     }
  234.  
  235.     while (k > 0)
  236.     {
  237.         mpz_urandomm(a, gmpRandState, n3);
  238.         mpz_add_ui(a, a, 2);
  239.         if (! is_spsp(n, a))
  240.         {
  241.             mpz_clears(a, n3, t, NULL);
  242.             return 0;
  243.         }
  244.         k -= 1;
  245.     }
  246.  
  247.     mpz_clears(a, n3, t, NULL);
  248.     return 1;
  249. }
  250.  
  251. /* integer factorization by pollard rho */
  252.  
  253. List *insert_in_order(void *x, List *xs)
  254. {
  255.     if (xs == NULL || mpz_cmp(x, xs->data) < 0)
  256.     {
  257.         return insert(x, xs);
  258.     }
  259.     else
  260.     {
  261.         List *head = xs;
  262.         while (xs->next != NULL && mpz_cmp(x, xs->next->data) > 0)
  263.         {
  264.             xs = xs->next;
  265.         }
  266.         xs->next = insert(x, xs->next);
  267.         return head;
  268.     }
  269. }
  270.  
  271. void rho_factor(mpz_t f, mpz_t n, long unsigned c)
  272. {
  273.     mpz_t t, h, d, r;
  274.  
  275.     mpz_init_set_ui(t, 2);
  276.     mpz_init_set_ui(h, 2);
  277.     mpz_init_set_ui(d, 1);
  278.     mpz_init(r);
  279.  
  280.     while (mpz_cmp_si(d, 1) == 0)
  281.     {
  282.         mpz_mul(t, t, t);
  283.         mpz_add_ui(t, t, c);
  284.         mpz_mod(t, t, n);
  285.  
  286.         mpz_mul(h, h, h);
  287.         mpz_add_ui(h, h, c);
  288.         mpz_mod(h, h, n);
  289.  
  290.         mpz_mul(h, h, h);
  291.         mpz_add_ui(h, h, c);
  292.         mpz_mod(h, h, n);
  293.  
  294.         mpz_sub(r, t, h);
  295.         mpz_gcd(d, r, n);
  296.     }
  297.  
  298.     if (mpz_cmp(d, n) == 0) /* cycle */
  299.     {
  300.         rho_factor(f, n, c+1);
  301.     }
  302.     else if (mpz_probab_prime_p(d, 25)) /* success */
  303.     {
  304.         mpz_set(f, d);
  305.     }
  306.     else /* found composite factor */
  307.     {
  308.         rho_factor(f, d, c+1);
  309.     }
  310.  
  311.     mpz_clears(t, h, d, r, NULL);
  312. }
  313.  
  314. void rho_factors(List **fs, mpz_t n)
  315. {
  316.  
  317.     while (mpz_even_p(n))
  318.     {
  319.         mpz_t *f = malloc(sizeof(*f));
  320.         mpz_init_set_ui(*f, 2);
  321.         *fs = insert(*f, *fs);
  322.         mpz_divexact_ui(n, n, 2);
  323.     }
  324.  
  325.     if (mpz_cmp_ui(n, 1) == 0) return;
  326.  
  327.     while (! (mpz_probab_prime_p(n, 25)))
  328.     {
  329.         mpz_t *f = malloc(sizeof(*f));
  330.         mpz_init_set_ui(*f, 0);
  331.  
  332.         rho_factor(*f, n, 1);
  333.         *fs = insert_in_order(*f, *fs);
  334.         mpz_divexact(n, n, *f);
  335.     }
  336.  
  337.     *fs = insert_in_order(n, *fs);
  338. }
  339.  
  340. /* demonstration */
  341.  
  342. int main(int argc, char *argv[])
  343. {
  344.     mpz_t n;
  345.     mpz_init(n);
  346.     List *ps = NULL;
  347.     List *fs = NULL;
  348.  
  349.     ps = primes(100); /* 2 3 5 7 11 13 17 19 23 29 31 37 41 */
  350.     while (ps != NULL) /* 43 47 53 59 61 67 71 73 79 83 89 97 */
  351.     {
  352.         printf("%ld%s", (long) ps->data,
  353.             (ps->next == NULL) ? "\n" : " ");
  354.         ps = ps->next;
  355.     }
  356.  
  357.     printf("%d\n", length(primes(1000000))); /* 78498 */
  358.  
  359.     printf("%d\n", td_prime(600851475143LL)); /* composite */
  360.  
  361.     fs = td_factors(600851475143LL); /* 71 839 1471 6857 */
  362.     while (fs != NULL)
  363.     {
  364.         printf("%llu%s", (unsigned long long int) fs->data,
  365.             (fs->next == NULL) ? "\n" : " ");
  366.         fs = fs->next;
  367.     }
  368.  
  369.     mpz_t a;
  370.     mpz_init(a);
  371.     mpz_set_str(n, "2047", 10);
  372.     mpz_set_str(a, "2", 10);
  373.     printf("%d\n", is_spsp(n, a)); /* pseudo-prime */
  374.  
  375.     mpz_set_str(n, "600851475143", 10); /* composite */
  376.     printf("%d\n", is_prime(n));
  377.  
  378.     mpz_set_str(n, "2305843009213693951", 10); /* prime */
  379.     printf("%d\n", is_prime(n));
  380.  
  381.     mpz_set_str(n, "600851475143", 10);
  382.     rho_factors(&fs, n); /* 71 839 1471 6857 */
  383.     while (fs != NULL) {
  384.         printf("%s%s", mpz_get_str(NULL, 10, fs->data),
  385.             (fs->next == NULL) ? "\n" : " ");
  386.         fs = fs->next;
  387.     }
  388. }
Not running #stdin #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/VHDJL
language:
C (gcc 8.3)
created:
12 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!