123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577 |
-
- // This value is tuned for generating large keys, or small keys many times. Higher values see
- // diminishing returns in efficiency (and though I cannot prove it, I believe that there is no
- // upper bound on the efficiency attainable by setting it higher, but it just isn't worth it.)
- // The larger this number, the longer build_sieve() will take (and thus make_difflist()). Also,
- // larger values of this number will use more memory. I wouldn't recommend setting it any higher,
- // but for generating small numbers of small keys, setting it lower is a good idea.
- // Another thing I forgot: Larger values will tend to make multithreaded searches where the
- // number of threads is not a multiple of many small primes whose primorial is less then this
- // number take much longer and use more memory.
- #define PRIMORIAL_CAP_CAP 1000000
-
- #include <stdlib.h>
- #include <stdio.h>
- #include <pthread.h>
- #include <errno.h>
-
- #include <gmp.h>
-
- struct mpz_list
- {
- unsigned long int numel;
- unsigned long int capacity;
- mpz_t *data;
- };
-
- void mpz_list_init(struct mpz_list *list, unsigned long int capacity)
- {
- list->numel = 0;
- list->capacity = capacity;
- list->data = malloc(capacity * sizeof(mpz_t));
- }
-
- // Push an element to the list
- mpz_t *mpz_list_push(struct mpz_list *list)
- {
- if (list->numel == list->capacity)
- {
- list->capacity <<= 1;
- list->data = realloc(list->data, list->capacity * sizeof(mpz_t));
- }
- mpz_init(list->data[list->numel]);
- list->numel++;
- return &(list->data[list->numel - 1]);
- }
-
- // Free the memory in a list
- void destroy_mpz_list(struct mpz_list *list)
- {
- int ii;
- for (ii = 0; ii < list->numel; ii++)
- mpz_clear(list->data[ii]);
- free(list->data);
- }
-
- struct ui_list
- {
- unsigned long int numel;
- unsigned long int capacity;
- unsigned long int *data;
- };
-
- void ui_list_init(struct ui_list *list, unsigned long int capacity)
- {
- list->numel = 0;
- list->capacity = capacity;
- list->data = malloc(capacity * sizeof(unsigned long int));
- }
-
- // Push an element to the list
- unsigned long int *ui_list_push(struct ui_list *list)
- {
- if (list->numel == list->capacity)
- {
- list->capacity <<= 1;
- list->data = realloc(list->data, list->capacity * sizeof(unsigned long int));
- }
- list->numel++;
- return &(list->data[list->numel - 1]);
- }
-
- // Free the memory in a list
- void destroy_ui_list(struct ui_list *list)
- {
- free(list->data);
- }
-
- // Reuse all the variables to prevent from having to allocate memory
- struct sieve_state
- {
- mpz_t root;
- };
-
- // Use a sieve to test primality of n
- int test_with_sieve(struct mpz_list *sieve, struct sieve_state state, mpz_t n)
- {
- int ii;
- mpz_sqrt(state.root, n);
- int iters = mpz_get_ui(state.root);
- if (sieve->numel < iters)
- iters = sieve->numel;
- for (ii = 0; ii < iters; ii++)
- {
- if (mpz_divisible_p(n, sieve->data[ii]))
- return 0;
- }
- return 1;
- }
-
- // Build a sieve until the primorial of the value being tested is above a certain value
- unsigned long int build_sieve(unsigned long int max_primorial, struct mpz_list **sieve_out)
- {
- struct mpz_list *sieve = malloc(sizeof(struct mpz_list));
- mpz_list_init(sieve, 16);
-
- struct sieve_state state;
- mpz_init(state.root);
-
- unsigned long int primorial = 1;
-
- mpz_t c;
- mpz_init(c);
- mpz_set_ui(c, 2);
-
- while (1)
- {
- unsigned long int ic = mpz_get_ui(c);
- if (test_with_sieve(sieve, state, c))
- {
- if (primorial * ic > max_primorial)
- break;
- primorial *= ic;
- mpz_set(*mpz_list_push(sieve), c);
- }
-
- mpz_add_ui(c, c, 1);
- }
-
- mpz_clear(state.root);
- mpz_clear(c);
- *sieve_out = sieve;
- return primorial;
- }
-
- //unsigned long int build_totant_list(unsigned long int n, struct mpz_list **list_out)
- //{
- // struct mpz_list *list = malloc(sizeof(struct mpz_list));
- // mpz_list_init(list, 16);
- //
- // mpz_t c, g;
- // mpz_init(c);
- // mpz_init(g);
- // mpz_set_ui(c, 1);
- //
- // while (mpz_cmp_ui(c, n) != 0)
- // {
- // mpz_gcd_ui(g, c, n);
- // if (mpz_cmp_ui(g, 1) == 0)
- // {
- // mpz_set(*mpz_list_push(list), c);
- // }
- //
- // mpz_add_ui(c, c, 1);
- // }
- //
- // mpz_clear(c);
- // mpz_clear(g);
- //
- // *list_out = list;
- // return list->numel;
- //}
-
- unsigned long int build_sgprime_candidate_list(unsigned long int n, struct mpz_list **list_out)
- {
- struct mpz_list *list = malloc(sizeof(struct mpz_list));
- mpz_list_init(list, 16);
-
- mpz_t c, g, u;
- mpz_init(c);
- mpz_init(g);
- mpz_init(u);
- mpz_set_ui(c, 1);
-
- while (mpz_cmp_ui(c, n) != 0)
- {
- mpz_gcd_ui(g, c, n);
- if (mpz_cmp_ui(g, 1) == 0)
- {
- mpz_add(u, c, c);
- mpz_add_ui(u, u, 1);
- if (mpz_cmp_ui(u, n) > 0)
- mpz_sub_ui(u, u, n);
- mpz_gcd_ui(g, u, n);
- if (mpz_cmp_ui(g, 1) == 0)
- {
- mpz_set(*mpz_list_push(list), c);
- }
- }
-
- mpz_add_ui(c, c, 1);
- }
-
- mpz_clear(c);
- mpz_clear(g);
- mpz_clear(u);
-
- *list_out = list;
- return list->numel;
- }
-
- //unsigned long int build_sgnprime_candidate_list(int niter, unsigned long int n, struct mpz_list **list_out)
- //{
- // struct mpz_list *list = malloc(sizeof(struct mpz_list));
- // mpz_list_init(list, 16);
- //
- // mpz_t c, g, u;
- // mpz_init(c);
- // mpz_init(g);
- // mpz_init(u);
- // mpz_set_ui(c, 1);
- //
- // while (mpz_cmp_ui(c, n) != 0)
- // {
- // mpz_set(u, c);
- // mpz_gcd_ui(g, u, n);
- // int ii;
- // for (ii = 0; ii < niter ; ii++)
- // {
- // if (mpz_cmp_ui(g, 1) != 0)
- // goto cont;
- // mpz_add(u, u, u);
- // mpz_add_ui(u, u, 1);
- // if (mpz_cmp_ui(u, n) > 0)
- // mpz_sub_ui(u, u, n);
- // mpz_gcd_ui(g, u, n);
- // }
- // mpz_set(*mpz_list_push(list), c);
- //cont:
- // mpz_add_ui(c, c, 1);
- // }
- //
- // mpz_clear(c);
- // mpz_clear(g);
- // mpz_clear(u);
- //
- // *list_out = list;
- // return list->numel;
- //}
-
- // This function makes the liberal assumption that each element has another one within the range of an unsigned long ahead of it.
- unsigned long int differentiate_list(unsigned long int max, struct mpz_list *list, struct ui_list **out_list)
- {
- if (list->numel <= 0) return 0;
-
- struct ui_list *difflist = malloc(sizeof(struct ui_list));
- ui_list_init(difflist, list->numel);
-
- mpz_t diff;
- mpz_init(diff);
-
- unsigned long int offset = mpz_get_ui(list->data[0]);
-
- unsigned long int head;
- for (head = 0; head < list->numel - 1; head++)
- {
- mpz_sub(diff, list->data[head + 1], list->data[head]);
- *ui_list_push(difflist) = mpz_get_ui(diff);
- }
-
- mpz_set(diff, list->data[0]);
- mpz_add_ui(diff, diff, max);
- mpz_sub(diff, diff, list->data[list->numel - 1]);
- *ui_list_push(difflist) = mpz_get_ui(diff);
-
- mpz_clear(diff);
-
- *out_list = difflist;
- return offset;
- }
-
- unsigned long int gcd_ui(unsigned long int a, unsigned long int b)
- {
- unsigned long int r;
- r = a % b;
- while (r)
- {
- a = b;
- b = r;
- r = a % b;
- }
- return b;
- }
-
- unsigned long int thread_separate_difflist(int nthreads, struct ui_list *difflist, unsigned int tid, struct ui_list **out_list)
- {
- unsigned long int llength = difflist->numel;
- unsigned long int plength = llength / gcd_ui(llength, nthreads);
-
- struct ui_list *list = malloc(sizeof(struct ui_list));
- ui_list_init(list, plength);
-
- unsigned int head = tid % llength;
- unsigned int ii, ij;
-
- for (ii = 0; ii < plength; ii++)
- {
- unsigned long int cum = 0;
- for (ij = 0; ij < nthreads; ij++)
- {
- cum += difflist->data[head];
- head = (head + 1) % llength;
- }
- *ui_list_push(list) = cum;
- }
-
- *out_list = list;
-
- unsigned long int offset = 0;
- for (head = 0; head < tid; head++)
- {
- offset += difflist->data[head % llength];
- }
-
- return offset;
- }
-
- // q is passed to prevent unnecessary memory allocation
- int test_sgprime(mpz_t p, mpz_t q, int reps)
- {
- if (mpz_probab_prime_p(p, reps))
- {
- mpz_add(q, p, p);
- mpz_add_ui(q, q, 1);
- if (mpz_probab_prime_p(q, reps))
- {
- return 1;
- }
- }
- return 0;
- }
-
- // q is passed to prevent unnecessary memory allocation
- int search_sgprime_iter(mpz_t a, mpz_t q, int reps, struct ui_list *difflist, int *counter, int chunksize)
- {
- int ii;
- for (ii = 0; ii < chunksize; ii++)
- {
- if (test_sgprime(a, q, reps))
- return 1;
- mpz_add_ui(a, a, difflist->data[(*counter)++]);
- if (*counter == difflist->numel)
- *counter = 0;
- }
- return 0;
- }
-
- void st_sgprime_search(mpz_t r, int reps, unsigned long int offset, struct ui_list *difflist, unsigned long int maximum)
- {
- mpz_t q;
- mpz_init(q);
-
- mpz_mod_ui(q, r, maximum);
- mpz_sub(r, r, q);
- mpz_add_ui(r, r, offset);
-
- int c = 0;
-
- while(!search_sgprime_iter(r, q, reps, difflist, &c, difflist->numel));
-
- mpz_clear(q);
- }
-
- 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)
- {
- mpz_t q, r;
- mpz_init(r);
- mpz_init(q);
-
- mpz_set(r, p);
-
- struct ui_list *tl_difflist;
- int tl_offset = thread_separate_difflist(nthreads, difflist, tid, &tl_difflist);
-
- mpz_mod_ui(q, r, maximum);
- mpz_sub(r, r, q);
- mpz_add_ui(r, r, offset + tl_offset);
-
- int c = 0;
-
- while(!search_sgprime_iter(r, q, reps, tl_difflist, &c, chunksize))
- {
- if (pthread_mutex_trylock(success_mutex) == 0)
- {
- if (*done)
- {
- pthread_mutex_unlock(success_mutex);
- goto cleanup;
- }
- pthread_mutex_unlock(success_mutex);
- }
- }
-
- pthread_mutex_lock(success_mutex);
- if (!*done) // We're the first ones here
- mpz_set(p, r);
- *done = 1;
- pthread_mutex_unlock(success_mutex);
-
- cleanup:
- destroy_ui_list(tl_difflist);
- free(tl_difflist);
- mpz_clear(q);
- mpz_clear(r);
- }
-
- struct mt_sgprime_search_kernel_params
- {
- 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;
- };
-
- void *mt_sgprime_search_kernel_entry(void *v)
- {
- struct mt_sgprime_search_kernel_params *p = (struct mt_sgprime_search_kernel_params*) v;
- 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);
- free(p);
- return NULL;
- }
-
- 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)
- {
- pthread_mutex_t *success_mutex = malloc(sizeof(pthread_mutex_t));
- pthread_mutex_init(success_mutex, NULL);
-
- pthread_t *threads = malloc(nthreads * sizeof(pthread_t));
-
- int done = 0;
-
- int ii;
- for (ii = 0; ii < nthreads; ii++)
- {
- struct mt_sgprime_search_kernel_params *p = malloc(sizeof(struct mt_sgprime_search_kernel_params));
- p->p[0] = r[0];
- p->reps = reps;
- p->offset = offset;
- p->difflist = difflist;
- p->maximum = maximum;
- p->nthreads = nthreads;
- p->tid = ii;
- p->success_mutex = success_mutex;
- p->done = &done;
- p->chunksize = chunksize;
- pthread_create(&threads[ii], NULL, &mt_sgprime_search_kernel_entry, p);
- }
-
- void *status;
-
- for (ii = 0; ii < nthreads; ii++)
- {
- pthread_join(threads[ii], &status);
- }
-
- free(threads);
- pthread_mutex_destroy(success_mutex);
- free(success_mutex);
- }
-
- struct ui_list *make_difflist(mpz_t lower_bound, unsigned long int *offset, unsigned long int *primorial)
- {
- struct mpz_list *list;
- unsigned long int primorial_cap;
- if (mpz_cmp_ui(lower_bound, PRIMORIAL_CAP_CAP) >= 0)
- {
- primorial_cap = PRIMORIAL_CAP_CAP;
- } else
- {
- primorial_cap = mpz_get_ui(lower_bound);
- }
- *primorial = build_sieve(primorial_cap, &list);
- destroy_mpz_list(list);
- free(list);
-
- unsigned long int totient = build_sgprime_candidate_list(*primorial, &list);
- struct ui_list *difflist;
- *offset = differentiate_list(*primorial, list, &difflist);
-
- destroy_mpz_list(list);
- free(list);
-
- return difflist;
- }
-
- void get_sg_prime(struct ui_list *difflist, unsigned long int offset, unsigned long int primorial, mpz_t search, unsigned int threads, unsigned int reps, unsigned int chunksize)
- {
- if (threads <= 1)
- {
- st_sgprime_search(search, reps, offset, difflist, primorial);
- } else
- {
- mt_sgprime_search(search, reps, offset, difflist, primorial, threads, chunksize);
- }
- }
-
- void destroy_difflist(struct ui_list *difflist)
- {
- destroy_ui_list(difflist);
- free(difflist);
- }
-
- //int main()
- //{
- // struct mpz_list *list;
- //
- // unsigned long int primorial = build_sieve(100000000, &list);
- //
- //// int ii;
- //// for (ii = 0; ii < list->numel; ii++)
- //// {
- //// printf("%lu\n", mpz_get_ui(list->data[ii]));
- //// }
- // // printf("%lu\n", primorial);
- //
- // destroy_mpz_list(list);
- // free(list);
- //
- // unsigned long int totient = build_sgprime_candidate_list(primorial, &list);
- //
- //// printf("finished building candidate list\n");
- //// fflush(stdout);
- //// printf("totient: %lu\n", totient);
- //// printf("ratio: %lu\n", primorial / totient);
- //
- // struct ui_list *difflist;
- //
- // unsigned long int offset = differentiate_list(primorial, list, &difflist);
- //
- //// printf("offset: %lu\n", offset);
- //
- // mpz_t p;
- // mpz_init(p);
- //
- // gmp_randstate_t rs;
- // gmp_randinit_default(rs);
- //
- // int bits = 32;
- //
- // int ii;
- // for (ii = 0; ii < 1000; ii++)
- // {
- // mpz_urandomb(p, rs, bits);
- // mpz_setbit(p, bits);
- //
- // mt_sgprime_search(p, 31, offset, difflist, primorial, 8, 10);
- //
- // gmp_printf("%i-bit safe prime: %Zd\n", bits, p);
- // }
- //
- // gmp_randclear(rs);
- //
- // mpz_clear(p);
- //
- // destroy_ui_list(difflist);
- // free(difflist);
- // destroy_mpz_list(list);
- // free(list);
- //
- // return 0;
- //}
|