Browse Source

Wrote a bunch of code to find safe primes

thajohns 11 months ago
commit
635a2ab349
1 changed files with 483 additions and 0 deletions
  1. 483
    0
      sieve.c

+ 483
- 0
sieve.c View File

@@ -0,0 +1,483 @@
1
+
2
+#include <stdlib.h>
3
+#include <stdio.h>
4
+#include <pthread.h>
5
+#include <errno.h>
6
+
7
+#include <gmp.h>
8
+
9
+#include "sieve.h"
10
+
11
+struct mpz_list
12
+{
13
+  unsigned long int numel;
14
+  unsigned long int capacity;
15
+  mpz_t *data;
16
+};
17
+
18
+void mpz_list_init(struct mpz_list *list, unsigned long int capacity)
19
+{
20
+  list->numel = 0;
21
+  list->capacity = capacity;
22
+  list->data = malloc(capacity * sizeof(mpz_t));
23
+}
24
+
25
+// Push an element to the list
26
+mpz_t *mpz_list_push(struct mpz_list *list)
27
+{
28
+      if (list->numel == list->capacity)
29
+      {
30
+        list->capacity <<= 1;
31
+        list->data = realloc(list->data, list->capacity * sizeof(mpz_t));
32
+      }
33
+      mpz_init(list->data[list->numel]);
34
+      list->numel++;
35
+      return &(list->data[list->numel - 1]);
36
+}
37
+
38
+// Free the memory in a list
39
+void destroy_mpz_list(struct mpz_list *list)
40
+{
41
+  int ii;
42
+  for (ii = 0; ii < list->numel; ii++)
43
+    mpz_clear(list->data[ii]);
44
+  free(list->data);
45
+}
46
+
47
+struct ui_list
48
+{
49
+  unsigned long int numel;
50
+  unsigned long int capacity;
51
+  unsigned long int *data;
52
+};
53
+
54
+void ui_list_init(struct ui_list *list, unsigned long int capacity)
55
+{
56
+  list->numel = 0;
57
+  list->capacity = capacity;
58
+  list->data = malloc(capacity * sizeof(unsigned long int));
59
+}
60
+
61
+// Push an element to the list
62
+unsigned long int *ui_list_push(struct ui_list *list)
63
+{
64
+      if (list->numel == list->capacity)
65
+      {
66
+        list->capacity <<= 1;
67
+        list->data = realloc(list->data, list->capacity * sizeof(unsigned long int));
68
+      }
69
+      list->numel++;
70
+      return &(list->data[list->numel - 1]);
71
+}
72
+
73
+// Free the memory in a list
74
+void destroy_ui_list(struct ui_list *list)
75
+{
76
+  free(list->data);
77
+}
78
+
79
+// Reuse all the variables to prevent from having to allocate memory
80
+struct sieve_state
81
+{
82
+  mpz_t root;
83
+};
84
+
85
+// Use a sieve to test primality of n
86
+int test_with_sieve(struct mpz_list *sieve, struct sieve_state state, mpz_t n)
87
+{
88
+  int ii;
89
+  mpz_sqrt(state.root, n);
90
+  int iters = mpz_get_ui(state.root);
91
+  if (sieve->numel < iters)
92
+    iters = sieve->numel;
93
+  for (ii = 0; ii < iters; ii++)
94
+  {
95
+    if (mpz_divisible_p(n, sieve->data[ii]))
96
+      return 0;
97
+  }
98
+  return 1;
99
+}
100
+
101
+// Build a sieve until the primorial of the value being tested is above a certain value
102
+unsigned long int build_sieve(unsigned long int max_primorial, struct mpz_list **sieve_out)
103
+{
104
+  struct mpz_list *sieve = malloc(sizeof(struct mpz_list));
105
+  mpz_list_init(sieve, 16);
106
+
107
+  struct sieve_state state;
108
+  mpz_init(state.root);
109
+
110
+  unsigned long int primorial = 1;
111
+
112
+  mpz_t c;
113
+  mpz_init(c);
114
+  mpz_set_ui(c, 2);
115
+
116
+  while (1)
117
+  {
118
+    unsigned long int ic = mpz_get_ui(c);
119
+    if (test_with_sieve(sieve, state, c))
120
+    {
121
+      if (primorial * ic > max_primorial)
122
+        break;
123
+      primorial *= ic;
124
+      mpz_set(*mpz_list_push(sieve), c);
125
+    }
126
+
127
+    mpz_add_ui(c, c, 1);
128
+  }
129
+
130
+  mpz_clear(state.root);
131
+  mpz_clear(c);
132
+  *sieve_out = sieve;
133
+  return primorial;
134
+}
135
+
136
+unsigned long int build_totant_list(unsigned long int n, struct mpz_list **list_out)
137
+{
138
+  struct mpz_list *list = malloc(sizeof(struct mpz_list));
139
+  mpz_list_init(list, 16);
140
+
141
+  mpz_t c, g;
142
+  mpz_init(c);
143
+  mpz_init(g);
144
+  mpz_set_ui(c, 1);
145
+
146
+  while (mpz_cmp_ui(c, n) != 0)
147
+  {
148
+    mpz_gcd_ui(g, c, n);
149
+    if (mpz_cmp_ui(g, 1) == 0)
150
+    {
151
+      mpz_set(*mpz_list_push(list), c);
152
+    }
153
+
154
+    mpz_add_ui(c, c, 1);
155
+  }
156
+
157
+  mpz_clear(c);
158
+  mpz_clear(g);
159
+
160
+  *list_out = list;
161
+  return list->numel;
162
+}
163
+
164
+unsigned long int build_sgprime_candidate_list(unsigned long int n, struct mpz_list **list_out)
165
+{
166
+  struct mpz_list *list = malloc(sizeof(struct mpz_list));
167
+  mpz_list_init(list, 16);
168
+
169
+  mpz_t c, g, u;
170
+  mpz_init(c);
171
+  mpz_init(g);
172
+  mpz_init(u);
173
+  mpz_set_ui(c, 1);
174
+
175
+  while (mpz_cmp_ui(c, n) != 0)
176
+  {
177
+    mpz_gcd_ui(g, c, n);
178
+    if (mpz_cmp_ui(g, 1) == 0)
179
+    {
180
+      mpz_add(u, c, c);
181
+      mpz_add_ui(u, u, 1);
182
+      if (mpz_cmp_ui(u, n) > 0)
183
+        mpz_sub_ui(u, u, n);
184
+      mpz_gcd_ui(g, u, n);
185
+      if (mpz_cmp_ui(g, 1) == 0)
186
+      {
187
+        mpz_set(*mpz_list_push(list), c);
188
+      }
189
+    }
190
+
191
+    mpz_add_ui(c, c, 1);
192
+  }
193
+
194
+  mpz_clear(c);
195
+  mpz_clear(g);
196
+  mpz_clear(u);
197
+
198
+  *list_out = list;
199
+  return list->numel;
200
+}
201
+
202
+// This function makes the liberal assumption that each element has another one within the range of an unsigned long ahead of it.
203
+unsigned long int differentiate_list(unsigned long int max, struct mpz_list *list, struct ui_list **out_list)
204
+{
205
+  if (list->numel <= 0) return 0;
206
+
207
+  struct ui_list *difflist = malloc(sizeof(struct ui_list));
208
+  ui_list_init(difflist, list->numel);
209
+
210
+  mpz_t diff;
211
+  mpz_init(diff);
212
+
213
+  unsigned long int offset = mpz_get_ui(list->data[0]);
214
+
215
+  unsigned long int head;
216
+  for (head = 0; head < list->numel - 1; head++)
217
+  {
218
+    mpz_sub(diff, list->data[head + 1], list->data[head]);
219
+    *ui_list_push(difflist) = mpz_get_ui(diff);
220
+  }
221
+
222
+  mpz_set(diff, list->data[0]);
223
+  mpz_add_ui(diff, diff, max);
224
+  mpz_sub(diff, diff, list->data[list->numel - 1]);
225
+  *ui_list_push(difflist) = mpz_get_ui(diff);
226
+
227
+  mpz_clear(diff);
228
+
229
+  *out_list = difflist;
230
+  return offset;
231
+}
232
+
233
+unsigned long int gcd_ui(unsigned long int a, unsigned long int b)
234
+{
235
+  unsigned long int r;
236
+  r = a % b;
237
+  while (r)
238
+  {
239
+    a = b;
240
+    b = r;
241
+    r = a % b;
242
+  }
243
+  return b;
244
+}
245
+
246
+unsigned long int thread_separate_difflist(int nthreads, struct ui_list *difflist, unsigned int tid, struct ui_list **out_list)
247
+{
248
+  unsigned long int llength = difflist->numel;
249
+  unsigned long int plength = llength / gcd_ui(llength, nthreads);
250
+
251
+  struct ui_list *list = malloc(sizeof(struct ui_list));
252
+  ui_list_init(list, plength);
253
+
254
+  unsigned int head = tid % llength;
255
+  unsigned int ii, ij;
256
+
257
+  for (ii = 0; ii < plength; ii++)
258
+  {
259
+    unsigned long int cum = 0;
260
+    for (ij = 0; ij < nthreads; ij++)
261
+    {
262
+      cum += difflist->data[head];
263
+      head = (head + 1) % llength;
264
+    }
265
+    *ui_list_push(list) = cum;
266
+  }
267
+
268
+  *out_list = list;
269
+
270
+  unsigned long int offset = 0;
271
+  for (head = 0; head < tid; head++)
272
+  {
273
+    offset += difflist->data[head % llength];
274
+  }
275
+  
276
+  return offset;
277
+}
278
+
279
+// q is passed to prevent unnecessary memory allocation
280
+int test_sgprime(mpz_t p, mpz_t q, int reps)
281
+{
282
+  if (mpz_probab_prime_p(p, reps))
283
+  {
284
+    mpz_add(q, p, p);
285
+    mpz_add_ui(q, q, 1);
286
+    if (mpz_probab_prime_p(q, reps))
287
+    {
288
+      return 1;
289
+    }
290
+  }
291
+  return 0;
292
+}
293
+
294
+// Search through a difflist exactly once for an SG prime. Return nonzero if one is found.
295
+// q is passed to prevent unnecessary memory allocation
296
+int search_sgprime_iter(mpz_t a, mpz_t q, int reps, struct ui_list *difflist, int *counter, int chunksize)
297
+{
298
+  int ii;
299
+  for (ii = 0; ii < chunksize; ii++)
300
+  {
301
+    if (test_sgprime(a, q, reps))
302
+      return 1;
303
+    mpz_add_ui(a, a, difflist->data[(*counter)++]);
304
+    if (*counter == difflist->numel)
305
+      *counter = 0;
306
+  }
307
+  return 0;
308
+}
309
+
310
+void st_sgprime_search(mpz_t r, int reps, unsigned long int offset, struct ui_list *difflist, unsigned long int maximum)
311
+{
312
+  mpz_t q;
313
+  mpz_init(q);
314
+
315
+  mpz_mod_ui(q, r, maximum);
316
+  mpz_sub(r, r, q);
317
+  mpz_add_ui(r, r, offset);
318
+
319
+  int c;
320
+  while(!search_sgprime_iter(r, q, reps, difflist, &c, difflist->numel));
321
+
322
+  mpz_clear(q);
323
+}
324
+
325
+void mt_sgprime_search_kernel(mpz_t p, int reps, unsigned long int offset, struct ui_list *difflist, unsigned long int maximum, int nthreads, int tid, pthread_mutex_t *success_mutex, int *done, int chunksize)
326
+{
327
+  mpz_t q, r;
328
+  mpz_init(r);
329
+  mpz_init(q);
330
+
331
+  mpz_set(r, p);
332
+
333
+  struct ui_list *tl_difflist;
334
+  int tl_offset = thread_separate_difflist(nthreads, difflist, tid, &tl_difflist);
335
+
336
+  mpz_mod_ui(q, r, maximum);
337
+  mpz_sub(r, r, q);
338
+  mpz_add_ui(r, r, offset + tl_offset);
339
+
340
+  int c = 0;
341
+
342
+  while(!search_sgprime_iter(r, q, reps, tl_difflist, &c, chunksize))
343
+  {
344
+    if (pthread_mutex_trylock(success_mutex) == 0)
345
+    {
346
+      if (*done)
347
+      {
348
+        pthread_mutex_unlock(success_mutex);
349
+        goto cleanup;
350
+      }
351
+      pthread_mutex_unlock(success_mutex);
352
+    }
353
+  }
354
+
355
+  pthread_mutex_lock(success_mutex);
356
+  if (!*done)  // We're the first ones here
357
+    mpz_set(p, r);
358
+  *done = 1;
359
+  pthread_mutex_unlock(success_mutex);
360
+
361
+cleanup:
362
+  destroy_ui_list(tl_difflist);
363
+  free(tl_difflist);
364
+  mpz_clear(q);
365
+  mpz_clear(r);
366
+}
367
+
368
+struct mt_sgprime_search_kernel_params
369
+{
370
+  mpz_t p;
371
+  int reps;
372
+  unsigned long int offset;
373
+  struct ui_list *difflist;
374
+  unsigned long int maximum;
375
+  int nthreads;
376
+  int tid;
377
+  pthread_mutex_t *success_mutex;
378
+  int *done;
379
+  int chunksize;
380
+};
381
+
382
+void *mt_sgprime_search_kernel_entry(void *v)
383
+{
384
+  struct mt_sgprime_search_kernel_params *p = (struct mt_sgprime_search_kernel_params*) v;
385
+  mt_sgprime_search_kernel(p->p, p->reps, p->offset, p->difflist, p->maximum, p->nthreads, p->tid, p->success_mutex, p->done, p->chunksize);
386
+  free(p);
387
+  return NULL;
388
+}
389
+
390
+void mt_sgprime_search(mpz_t r, int reps, unsigned long int offset, struct ui_list *difflist, unsigned long int maximum, int nthreads, int chunksize)
391
+{
392
+  pthread_mutex_t *success_mutex = malloc(sizeof(pthread_mutex_t));
393
+  pthread_mutex_init(success_mutex, NULL);
394
+
395
+  pthread_t *threads = malloc(nthreads * sizeof(pthread_t));
396
+
397
+  int done = 0;
398
+
399
+  int ii;
400
+  for (ii = 0; ii < nthreads; ii++)
401
+  {
402
+    struct mt_sgprime_search_kernel_params *p = malloc(sizeof(struct mt_sgprime_search_kernel_params));
403
+    p->p[0] = r[0];
404
+    p->reps = reps;
405
+    p->offset = offset;
406
+    p->difflist = difflist;
407
+    p->maximum = maximum;
408
+    p->nthreads = nthreads;
409
+    p->tid = ii;
410
+    p->success_mutex = success_mutex;
411
+    p->done = &done;
412
+    p->chunksize = chunksize;
413
+    pthread_create(&threads[ii], NULL, &mt_sgprime_search_kernel_entry, p);
414
+  }
415
+
416
+  void *status;
417
+
418
+  for (ii = 0; ii < nthreads; ii++)
419
+  {
420
+    pthread_join(threads[ii], &status);
421
+  }
422
+
423
+  free(threads);
424
+  pthread_mutex_destroy(success_mutex);
425
+  free(success_mutex);
426
+}
427
+
428
+int main()
429
+{
430
+  struct mpz_list *list;
431
+
432
+  unsigned long int primorial = build_sieve(100000000, &list);
433
+
434
+//  int ii;
435
+//  for (ii = 0; ii < list->numel; ii++)
436
+//  {
437
+//    printf("%lu\n", mpz_get_ui(list->data[ii]));
438
+//  }
439
+  printf("%lu\n", primorial);
440
+
441
+  destroy_mpz_list(list);
442
+  free(list);
443
+
444
+  unsigned long int totient = build_sgprime_candidate_list(primorial, &list);
445
+
446
+  printf("totient: %lu\n", totient);
447
+  printf("ratio: %lu\n", primorial / totient);
448
+
449
+  struct ui_list *difflist;
450
+
451
+  unsigned long int offset = differentiate_list(primorial, list, &difflist);
452
+  
453
+  printf("offset: %lu\n", offset);
454
+
455
+  mpz_t p;
456
+  mpz_init(p);
457
+
458
+  gmp_randstate_t rs;
459
+  gmp_randinit_default(rs);
460
+
461
+  int ii;
462
+  for (ii = 0; ii < 5; ii++)
463
+  {
464
+    mpz_urandomb(p, rs, 4096);
465
+    mpz_setbit(p, 4096);
466
+
467
+    mt_sgprime_search(p, 31, offset, difflist, primorial, 8, 100);
468
+
469
+    gmp_printf("safe prime: %Zd\n", p);
470
+  }
471
+
472
+  gmp_randclear(rs);
473
+
474
+  mpz_clear(p);
475
+
476
+  destroy_ui_list(difflist);
477
+  free(difflist);
478
+  destroy_mpz_list(list);
479
+  free(list);
480
+
481
+  return 0;
482
+}
483
+

Loading…
Cancel
Save