proth.c

   1 #include <stdio.h>
   2 #include <string.h>
   3 #include <stdbool.h>
   4
   5 #include <time.h>
   6
   7 #include <gmp.h>
   8
   9 /* the real test */
  10 //#define K_CONST 21181
  11 //#define N_CONST 44200000
  12
  13 /* recent find (k=46157, n=698207)*/
  14 //#define K_CONST 46157
  15 //#define N_CONST 698206
  16
  17 /* searching for k=301, n=7360 */
  18 #define K_CONST 301
  19 #define N_CONST 6500
  20
  21 mpz_t t_exp;
  22 mpz_t t_p1;
  23 mpz_t t_p2;
  24 mpz_t t_p3;
  25 mpz_t t;
  26
  27 time_t start;
  28 time_t end;
  29
  30 int main(void)
  31 {
  32     time(&start);
  33
  34     mpz_t exp;
  35     mpz_t p;
  36     mpz_t k;
  37     mpz_t base;
  38
  39     mpz_init(p);
  40
  41     mpz_init(k);
  42     mpz_init(exp);
  43     mpz_init(base);
  44
  45     mpz_init(t_exp);
  46     mpz_init(t_p1);
  47     mpz_init(t_p2);
  48     mpz_init(t_p3);
  49     mpz_init(t);
  50
  51     //first 3 primes
  52     mpz_set_ui(t_p1, 2);
  53     mpz_set_ui(t_p2, 3);
  54     mpz_set_ui(t_p3, 5);
  55
  56     //proth number is k*2^n+1
  57     mpz_set_ui(k, K_CONST);
  58
  59     //two is our base
  60     mpz_set_ui(base, 2);
  61
  62     unsigned long int n = N_CONST;
  63     while (1)
  64     {
  65         printf("(n=%llu)\n", n);
  66
  67         //calc 2^n
  68         mpz_pow_ui(exp, base, n);
  69
  70         //p=k*2^n
  71         mpz_mul(p, k, exp);
  72         //p+=1
  73         mpz_add_ui(p, p, 1);
  74
  75         //prime check exponent is (p-1)/2
  76         mpz_sub_ui(t_exp, p, 1);
  77         mpz_divexact_ui(t_exp, t_exp, 2);
  78
  79         //printf("prime check a=2\n");
  80         mpz_powm(t, t_p1, t_exp, p);
  81         mpz_add_ui(t, t, 1);
  82         //printf("cmp1\n");
  83         if(mpz_cmp(t, p) == 0)
  84         {
  85             printf("prime with a=2!\n");
  86             break;
  87         }
  88
  89         //printf("prime check a=3\n");
  90         mpz_powm(t, t_p2, t_exp, p);
  91         mpz_add_ui(t, t, 1);
  92         //printf("cmp2\n");
  93         if(mpz_cmp(t, p) == 0)
  94         {
  95             printf("prime with a=3!\n");
  96             break;
  97         }
  98
  99         //printf("prime check a=5\n");
 100         mpz_powm(t, t_p3, t_exp, p);
 101         mpz_add_ui(t, t, 1);
 102         //printf("cmp3\n");
 103         if(mpz_cmp(t, p) == 0)
 104         {
 105             printf("prime with a=5!\n");
 106             break;
 107         }
 108
 109         n++;
 110     }
 111
 112     //proth number is k*2^n+1
 113     gmp_printf("%Zd*2^%llu+1 is prime!\n", k, n);
 114
 115     //free
 116     mpz_clear(p);
 117
 118     mpz_clear(k);
 119     mpz_clear(exp);
 120     mpz_clear(base);
 121
 122     mpz_clear(t_exp);
 123     mpz_clear(t_p1);
 124     mpz_clear(t_p2);
 125     mpz_clear(t_p3);
 126     mpz_clear(t);
 127
 128     time(&end);
 129     printf("total time: %f seconds\n", difftime(end, start));
 130 }