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