cs456-a3/sgprime.c

// 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;
//}