proth.c

   1 #include <stdio.h>
   2 #include <string.h>
   3 #include <stdbool.h>
   4
   5 #include <pthread.h>
   6 #include <time.h>
   7
   8 #include <gmp.h>
   9
  10 /* the real test */
  11 //#define K_CONST 21181
  12 //#define N_CONST 44200000
  13
  14 /* recent find (k=46157, n=698207)*/
  15 //#define K_CONST 46157
  16 //#define N_CONST 698203
  17
  18 /* searching for k=301, n=7360 */
  19 #define K_CONST 301
  20 #define N_CONST 6500
  21
  22 #define NUM_CORES 6
  23
  24 static void * proth_thread(void *arg);
  25
  26 //mutex protected, thread safe n
  27 unsigned long int n_val = N_CONST;
  28 pthread_mutex_t n_lock;
  29 pthread_mutex_t found_lock;
  30
  31 int main(void)
  32 {
  33     time_t start;
  34     time_t end;
  35     time(&start);
  36
  37     //array to keep track of the thread id's
  38     pthread_t thread_ids[NUM_CORES];
  39
  40     //create the mutexes
  41     pthread_mutex_init(&n_lock, NULL);
  42     pthread_mutex_init(&found_lock, NULL);
  43
  44     //create the threads
  45     for(int i = 0; i < NUM_CORES; i++)
  46     {
  47         pthread_create(&thread_ids[i], NULL, &proth_thread, NULL);
  48     }
  49
  50     //lock main thread until prime is found
  51     pthread_mutex_lock(&found_lock);
  52     pthread_mutex_lock(&found_lock);
  53
  54     //destroy the threads
  55     for(int i = 0; i < NUM_CORES; i++)
  56     {
  57         pthread_cancel(thread_ids[i]);
  58     }
  59
  60     //destroy the mutexes
  61     pthread_mutex_destroy(&n_lock);
  62     pthread_mutex_destroy(&found_lock);
  63
  64     time(&end);
  65     printf("total time: %f seconds\n", difftime(end, start));
  66 }
  67
  68 unsigned long get_n()
  69 {
  70     pthread_mutex_lock(&n_lock);
  71     unsigned long n = n_val;
  72     n_val++;
  73     pthread_mutex_unlock(&n_lock);
  74     return n;
  75 }
  76
  77 static void * proth_thread(void *arg)
  78 {
  79     // variables to test the primality
  80     mpz_t t;
  81     mpz_t t_exp;
  82     mpz_t t_p1;
  83     mpz_t t_p2;
  84     mpz_t t_p3;
  85
  86     mpz_init(t);
  87     mpz_init(t_exp);
  88     mpz_init(t_p1);
  89     mpz_init(t_p2);
  90     mpz_init(t_p3);
  91
  92     // variables to calculate the proth number p=k*2^n+1
  93     mpz_t proth;    //the proth number
  94     mpz_t p_minus1; // proth number minus 1
  95     mpz_t k;        // k
  96     mpz_t two;      // 2
  97     mpz_t two_n;     //2^n
  98
  99     mpz_init(proth);
 100     mpz_init(p_minus1);
 101     mpz_init(k);
 102     mpz_init(two);
 103     mpz_init(two_n);
 104
 105     //first 3 primes
 106     mpz_set_ui(t_p1, 2);
 107     mpz_set_ui(t_p2, 3);
 108     mpz_set_ui(t_p3, 5);
 109
 110     //proth number is p=k*2^n+1
 111     mpz_set_ui(k, K_CONST);
 112
 113     //two=2
 114     mpz_set_ui(two, 2);
 115
 116     unsigned long int n;
 117
 118     while (1)
 119     {
 120         //get the thread-safe n
 121         n = get_n();
 122         printf("(n=%lu)\n", n);
 123
 124         //calc exp=2^n
 125         mpz_pow_ui(two_n, two, n);
 126         //calc p=k*2^n
 127         mpz_mul(p_minus1, k, two_n);
 128         //calc p=k*2^n+1 (this is our proth number)
 129         mpz_add_ui(proth, p_minus1, 1);
 130
 131         //the prime check exponent is (p-1)/2
 132         //mpz_sub_ui(t_exp, p, 1);
 133         mpz_divexact_ui(t_exp, p_minus1, 2);
 134
 135 #if 0
 136         //printf("prime check a=2\n");
 137         mpz_powm(t, t_p1, t_exp, proth);
 138         mpz_add_ui(t, t, 1);
 139         //printf("cmp1\n");
 140         if(mpz_cmp(t, proth) == 0)
 141         {
 142             printf("prime with a=2!\n");
 143             break;
 144         }
 145 #endif
 146
 147         //printf("prime check a=3\n");
 148         mpz_powm(t, t_p2, t_exp, proth);
 149         mpz_add_ui(t, t, 1);
 150         //printf("cmp2\n");
 151         if(mpz_cmp(t, proth) == 0)
 152         {
 153             printf("prime with a=3!\n");
 154             break;
 155         }
 156
 157         //printf("prime check a=5\n");
 158         mpz_powm(t, t_p3, t_exp, proth);
 159         mpz_add_ui(t, t, 1);
 160         //printf("cmp3\n");
 161         if(mpz_cmp(t, proth) == 0)
 162         {
 163             printf("prime with a=5!\n");
 164             break;
 165         }
 166     }
 167
 168     //proth number is k*2^n+1
 169     gmp_printf("%Zd*2^%llu+1 is prime!\n", k, n);
 170
 171     //free
 172     mpz_clear(proth);
 173     mpz_init(p_minus1);
 174     mpz_clear(k);
 175     mpz_clear(two);
 176     mpz_clear(two_n);
 177
 178     mpz_clear(t_exp);
 179     mpz_clear(t_p1);
 180     mpz_clear(t_p2);
 181     mpz_clear(t_p3);
 182     mpz_clear(t);
 183
 184     pthread_mutex_unlock(&found_lock);
 185
 186     return NULL;
 187 }