proth.c

   1 #include <time.h>
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <stdint.h>
   5 #include <stdbool.h>
   6
   7 #include <gmp.h>
   8
   9
  10 bool is_prime(mpz_t p);
  11 bool is_prime2(mpz_t p);
  12
  13 clock_t start;
  14
  15
  16 //ints for is_prime test
  17 mpz_t t_exp;
  18 mpz_t t_p1;
  19 mpz_t t_p2;
  20 mpz_t t_p3;
  21 mpz_t t;
  22
  23
  24 int main(void)
  25 {
  26     mpz_t exp;
  27     mpz_t p;
  28     mpz_t k;
  29     mpz_t base;
  30
  31     mpz_init(p);
  32
  33     mpz_init(k);
  34     mpz_init(exp);
  35     mpz_init(base);
  36
  37     mpz_init(t_exp);
  38     mpz_init(t_p1);
  39     mpz_init(t_p2);
  40     mpz_init(t_p3);
  41     mpz_init(t);
  42
  43     //first 3 primes
  44     mpz_set_ui(t_p1, 2);
  45     mpz_set_ui(t_p2, 3);
  46     mpz_set_ui(t_p3, 5);
  47
  48     //start = clock();
  49
  50
  51     //proth number is k*2^n+1
  52     //check where k is 19
  53     mpz_set_ui(k, 46157);
  54
  55     //calc 2^n
  56     mpz_set_ui(base, 2);
  57
  58     for(unsigned long int n = 698200; ; n++)
  59     {
  60         mpz_pow_ui(exp, base, n);
  61
  62         //p=k*2^n
  63         mpz_mul(p, k, exp);
  64         //p+=1
  65         mpz_add_ui(p, p, 1);
  66
  67         gmp_printf("(n=%d)\n%Zd\n\n", n, p);
  68         if(is_prime2(p) == true)
  69         {
  70             printf(" <--- prime\n");
  71             return 1;
  72         }
  73         printf("\n");
  74     }
  75
  76     //free p
  77     mpz_clear(p);
  78 }
  79
  80 bool is_prime(mpz_t p)
  81 {
  82     //exponent is (p-1)/2
  83     //printf("exp\n");
  84     mpz_sub_ui(t_exp, p, 1);
  85     mpz_divexact_ui(t_exp, t_exp, 2);
  86     unsigned long int ui_exp = mpz_get_ui(t_exp);
  87
  88     //printf("pow1\n");
  89     mpz_pow_ui(t, t_p1, ui_exp);
  90     mpz_add_ui(t, t, 1);
  91     //printf("div1\n");
  92     if(mpz_divisible_p(t, p) != 0)
  93         return true;
  94
  95     //printf("pow2\n");
  96     mpz_pow_ui(t, t_p2, ui_exp);
  97     mpz_add_ui(t, t, 1);
  98     //printf("div2\n");
  99     if(mpz_divisible_p(t, p) != 0)
 100         return true;
 101
 102     //printf("pow3\n");
 103     mpz_pow_ui(t, t_p3, ui_exp);
 104     mpz_add_ui(t, t, 1);
 105     //printf("div3\n");
 106     if(mpz_divisible_p(t, p) != 0)
 107         return true;
 108
 109     return false;
 110 }
 111
 112 bool is_prime2(mpz_t p)
 113 {
 114     //exponent is (p-1)/2
 115     printf("exp\n");
 116     mpz_sub_ui(t_exp, p, 1);
 117     mpz_divexact_ui(t_exp, t_exp, 2);
 118
 119     printf("pow1\n");
 120     mpz_powm(t, t_p1, t_exp, p);
 121     mpz_add_ui(t, t, 1);
 122     printf("cmp1\n");
 123     if(mpz_cmp(t, p) == 0)
 124         return true;
 125
 126     printf("pow2\n");
 127     mpz_powm(t, t_p2, t_exp, p);
 128     mpz_add_ui(t, t, 1);
 129     printf("cmp2\n");
 130     if(mpz_cmp(t, p) == 0)
 131         return true;
 132
 133     printf("pow3\n");
 134     mpz_powm(t, t_p3, t_exp, p);
 135     mpz_add_ui(t, t, 1);
 136     printf("cmp3\n");
 137     if(mpz_cmp(t, p) == 0)
 138         return true;
 139
 140     return false;
 141 }