#include <stdio.h>
#include <stdlib.h>
#include <stdint.h> // for intptr_t
#include <string.h>
#include <math.h>
#include <gmp.h>
#include "pbc.h"
#include "pbc_utils.h"

struct cell_s {
  int ind;
  mpz_t data;
};
typedef struct cell_s *cell_ptr;

static cell_ptr newcell(void)
{
  cell_ptr res;
  res = pbc_malloc(sizeof(struct cell_s));
  //res->data = pbc_malloc(sizeof(mpz_t));
  //mpz_init(res->data);
  mpz_init(res->data);
  return res;
}

static void delcell(void *p)
{
  cell_ptr cp = p;
  mpz_clear(cp->data);
  pbc_free(p);
}

static int is_gen(mpz_t x, mpz_t q, darray_ptr fac, darray_ptr mul) {
  int result = 1;
  mpz_t z;
  mpz_t q1;
  int i;
  UNUSED_VAR(mul);

  mpz_init(z);
  mpz_init(q1);

  mpz_sub_ui(q1, q, 1);
  for (i=0; i<fac->count; i++) {
    mpz_divexact(z, q1, fac->item[i]);
    mpz_powm(z, x, z, q);
    if (!mpz_cmp_ui(z, 1)) {
      result = 0;
      break;
    }
  }

  mpz_clear(q1);
  mpz_clear(z);
  return result;
}

// Garner's Algorithm.
// See Algorithm 14.71, Handbook of Cryptography.
static void CRT(mpz_t x, mpz_ptr *v, mpz_ptr *m, int t) {
  mpz_t u;
  mpz_t C[t];
  int i, j;

  mpz_init(u);
  for (i=1; i<t; i++) {
    mpz_init(C[i]);
    mpz_set_ui(C[i], 1);
    for (j=0; j<i; j++) {
      mpz_invert(u, m[j], m[i]);
      mpz_mul(C[i], C[i], u);
      mpz_mod(C[i], C[i], m[i]);
    }
  }
  mpz_set(u, v[0]);
  mpz_set(x, u);
  for (i=1; i<t; i++) {
    mpz_sub(u, v[i], x);
    mpz_mul(u, u, C[i]);
    mpz_mod(u, u, m[i]);
    for (j=0; j<i; j++) {
      mpz_mul(u, u, m[j]);
    }
    mpz_add(x, x, u);
  }

  for (i=1; i<t; i++) mpz_clear(C[i]);
  mpz_clear(u);
}

//TODO: http://www.cecm.sfu.ca/CAG/abstracts/aaron27Jan06.pdf
//TODO: don't need to store last element of list in row[i]
//TODO: linked lists might be better than dynamic arrays (avoids memmove())
//TODO: allow holes in the table
//(if drought lasts too long)
void index_calculus_step1(mpz_t *ind, int r, mpz_t g, mpz_t q,
    darray_ptr fac, darray_ptr mul) {
  int count = 0;
  int i, j;
  mpz_t z, z0, z1;
  int relcount;
  unsigned int *prime = pbc_malloc(sizeof(unsigned int) * r);
  int bundlecount = (r - 10 + 19) / 20;
  mpz_t *bundle = pbc_malloc(sizeof(mpz_t) * bundlecount);
  int faci;
  mpz_t k, km;

  cell_ptr *rel = pbc_malloc(sizeof(cell_ptr) * r);
  cell_ptr *relm = pbc_malloc(sizeof(cell_ptr) * r);
  //''matrix'' is actually a list of matrices
  //(we solve over different moduli and combine using CRT)
  //darray_t **matrix = pbc_malloc(sizeof(darray_t *) * fac->count);
  darray_t *matrix[fac->count];
  int minfound[fac->count];

  for (i=0; i<r; i++) {
    rel[i] = newcell();
    relm[i] = newcell();
  }
  for (i=0; i<fac->count; i++) {
    //similarly ''row'' refers to a list of rows
    darray_t *row = pbc_malloc(sizeof(darray_t) * r);
    for (j=0; j<r; j++) {
      darray_init(row[j]);
    }
    matrix[i] = row;
    minfound[i] = 0;
  }

  mpz_init(k);
  mpz_init(km);
  mpz_init(z);
  mpz_init(z1);
  mpz_init(z0);

  printf("building prime table...\n");
  prime[0] = 2;
  mpz_set_ui(z, 2);
  for (i=1; i<r; i++) {
    mpz_nextprime(z, z);
    prime[i] = mpz_get_ui(z);
  }

  for (i=0; i<bundlecount; i++) {
    mpz_init(bundle[i]);
    mpz_set_ui(bundle[i], 1);
    for (j=0; j<20; j++) {
      int jj = 10 + 20 * i + j;
      if (jj >= r) break;
      mpz_mul_ui(bundle[i], bundle[i], prime[jj]);
    }
    element_printf("bundle %d: %Zd\n", i, bundle[i]);
  }
  printf("searching for r-smooth numbers\n");

  mpz_set_ui(z, 1);
  mpz_init(k);
  int try = 0;
  do {
    mpz_mul(z, z, g);
    mpz_mod(z, z, q);
    mpz_add_ui(k, k, 1);

    /*
    pbc_mpz_random(k, q);
    mpz_powm(z, g, k, q);
    */

    try++;

    mpz_set(z1, z);
    relcount = 0;
    for (i=0; i<10; i++) {
      if (i >= r) break;
      j = 0;
      while (mpz_divisible_ui_p(z1, prime[i])) {
        mpz_divexact_ui(z1, z1, prime[i]);
        j++;
      }
      if (j) {
        rel[relcount]->ind = i;
        mpz_set_ui(rel[relcount]->data, j);
        relcount++;
        if (!mpz_cmp_ui(z1, 1)) goto found;
      }
    }
    for (i=0; i<bundlecount; i++) {
      mpz_gcd(z0, bundle[i], z1);
      if (mpz_cmp_ui(z0, 1)) {
        int ii;
        for (ii = 0; ii < 20; ii++) {
          int jj = 10 + i * 20 + ii;
          if (jj >= r) break;
          j = 0;
          while (mpz_divisible_ui_p(z1, prime[jj])) {
            mpz_divexact_ui(z1, z1, prime[jj]);
            j++;
          }
          if (j) {
            rel[relcount]->ind = jj;
            mpz_set_ui(rel[relcount]->data, j);
            relcount++;
            if (!mpz_cmp_ui(z1, 1)) goto found;
          }
        }
      }
    }
    continue;
found:

/*
    printf("found r-smooth number after %d tries\n", try);

    gmp_printf("g^%Zd = %Zd:", k, z);
    for (i=0; i<relcount; i++) {
      gmp_printf(" %u:%Zd", rel[i]->ind, rel[i]->data);
    }
    printf("\n");
*/
    try = 0;

    for (faci=0; faci<fac->count; faci++) {
      darray_t *row = matrix[faci];
      mpz_ptr order = fac->item[faci];
      int relmcount = 0;
      mpz_mod(km, k, order);

      //gmp_printf("mod %Zd\n", order);
      for (i=0; i<relcount; i++) {
        mpz_mod(z0, rel[i]->data, order);
        if (mpz_sgn(z0)) {
          mpz_set(relm[relmcount]->data, z0);
          relm[relmcount]->ind = rel[i]->ind;
          relmcount++;
        }
      }

      while (relmcount) {
        //start from the sparse end
        int rind = relm[relmcount - 1]->ind;
        darray_ptr rp = row[rind];

        if (rind < minfound[faci]) break;

        mpz_set(z0, relm[relmcount - 1]->data);
        if (!rp->count) {
          mpz_invert(z0, z0, order);
          cell_ptr cnew = newcell();
          cnew->ind = -1;
          mpz_mul(z1, km, z0);
          mpz_mod(cnew->data, z1, order);
          darray_append(rp, cnew);
          for (j=0; j<relmcount; j++) {
            cnew = newcell();
            cnew->ind = relm[j]->ind;
            mpz_mul(z1, relm[j]->data, z0);
            mpz_mod(cnew->data, z1, order);
            darray_append(rp, cnew);
          }
          count++;
printf("%d / %d\n", count, r * fac->count);
/*
for (i=1; i<rp->count; i++) {
  cnew = rp->item[i];
  gmp_printf(" %u:%Zd", cnew->ind, cnew->data);
}
cnew = rp->item[0];
gmp_printf(" %Zd\n", cnew->data);
*/

          if (rind == minfound[faci]) {
            do {
              if (!minfound[faci]) {
              printf("found log p_%d\n", minfound[faci]);
              cnew = rp->item[0];
              gmp_printf("km = %Zd mod %Zd\n", cnew->data, order);
              }
              minfound[faci]++;
              if (minfound[faci] >= r) break;
              rp = row[minfound[faci]];
            } while (rp->count);
          }
          break;

        }

/*
{
//gmp_printf("mod = %Zd\n", order);
printf("before:");
for (i=0; i<relmcount; i++) {
  gmp_printf(" %u:%Zd", relm[i]->ind, relm[i]->data);
}
gmp_printf(" %Zd\n", km);
cell_ptr cp;
printf("sub %d:", rind);
for (i=1; i<rp->count; i++) {
  cp = rp->item[i];
  gmp_printf(" %u:%Zd", cp->ind, cp->data);
}
cp = rp->item[0];
gmp_printf(" %Zd\n", cp->data);
}
*/
        cell_ptr cpi, cpj;
        relmcount--;
        i=0; j=1;
        while (i<relmcount && j<rp->count - 1) {
          cpi = relm[i];
          cpj = rp->item[j];
          if (cpi->ind == cpj->ind) {
            mpz_mul(z1, z0, cpj->data);
            mpz_mod(z1, z1, order);
            int res = mpz_cmp(z1, cpi->data);
            if (!res) {
              memmove(&relm[i], &relm[i + 1], (relmcount - i - 1) * sizeof(cell_ptr));
              relm[relmcount - 1] = cpi;
              relmcount--;
              j++;
            } else if (res > 0) {
              mpz_sub(z1, order, z1);
              mpz_add(cpi->data, cpi->data, z1);
              i++;
              j++;
            } else {
              mpz_sub(cpi->data, cpi->data, z1);
              i++;
              j++;
            }
          } else if (cpi->ind > cpj->ind) {
            cpi = relm[relmcount];
            memmove(&relm[i + 1], &relm[i], (relmcount - i) * sizeof(cell_ptr));
            relm[i] = cpi;
            relmcount++;

            cpi->ind = cpj->ind;
            mpz_mul(z1, z0, cpj->data);
            mpz_mod(z1, z1, order);
            mpz_sub(cpi->data, order, z1);
            //cpi->data = order - ((u0 * cpj->data) % order);
            i++;
            j++;
          } else {
            i++;
          }
        }

        if (i == relmcount) {
          while (j < rp->count - 1) {
            cpi = relm[relmcount];
            cpj = rp->item[j];
            cpi->ind = cpj->ind;
            mpz_mul(z1, z0, cpj->data);
            mpz_mod(z1, z1, order);
            mpz_sub(cpi->data, order, z1);
            //cpi->data = order - ((u0 * cpj->data) % order);
            relmcount++;
            j++;
          }
        }

        cpj = rp->item[0];
        mpz_mul(z1, z0, cpj->data);
        mpz_mod(z1, z1, order);
        //u1 = (u0 * cpj->data) % order;
        if (mpz_cmp(km, z1) >= 0) {
          mpz_sub(km, km, z1);
        } else {
          mpz_sub(z1, order, z1);
          mpz_add(km, km, z1);
        }

/*
printf("after:");
for (i=0; i<relmcount; i++) {
  gmp_printf(" %u:%Zd", relm[i]->ind, relm[i]->data);
}
gmp_printf(" %Zd\n", km);
*/
      }
    }

  } while (count < r * fac->count);

  for (faci=0; faci<fac->count; faci++) {
    darray_t *row = matrix[faci];
    mpz_ptr order = fac->item[faci];
    for (i=1; i<r; i++) {
      darray_ptr rp = row[i];
      cell_ptr c0 = rp->item[0];
      for (j=1; j<rp->count-1; j++) {
        cell_ptr cp = rp->item[j];
        darray_ptr r2 = row[cp->ind];
        cell_ptr c2 = r2->item[0];
        mpz_mul(z0, cp->data, c2->data);
        mpz_sub(c0->data, c0->data, z0);
        mpz_mod(c0->data, c0->data, order);
      }
    }
  }

  mpz_ptr *tmp = pbc_malloc(sizeof(mpz_ptr) * fac->count);
  for (i=0; i<fac->count; i++) {
    tmp[i] = pbc_malloc(sizeof(mpz_t));
    mpz_init(tmp[i]);
    mpz_pow_ui(fac->item[i], fac->item[i], (unsigned int) mul->item[i]);
  }

  for (i=0; i<r; i++) {
    for (faci=0; faci<fac->count; faci++) {
      darray_t *row = matrix[faci];
      cell_ptr cp = row[i]->item[0];
      mpz_set(tmp[faci], cp->data);
    }
    CRT(ind[i], tmp, (mpz_ptr *) fac->item, fac->count);
  }

  for (i=0; i<fac->count; i++) {
    mpz_clear(tmp[i]);
  }
  pbc_free(tmp);

  for (faci=0; i<fac->count; faci++) {
    //similarly ''row'' refers to a list of rows
    darray_t *row = matrix[faci];
    for (j=0; j<r; j++) {
      darray_forall(row[j], delcell);
      darray_clear(row[j]);
    }
    pbc_free(matrix[faci]);
  }

  for (i=0; i<r; i++) {
    delcell(rel[i]);
    delcell(relm[i]);
  }

  pbc_free(prime);
  pbc_free(rel);
  pbc_free(relm);
  mpz_clear(k);
  mpz_clear(km);
  mpz_clear(z);
  mpz_clear(z0);
  mpz_clear(z1);
}

// Brute-force: does not use the fact that matrices are sparse.
void slow_index_calculus_step1(mpz_t *ind, int r, mpz_t g, mpz_t q,
    darray_ptr fac, darray_ptr mul) {
  int count = 0;
  int i, j;
  mpz_t z, z0, z1;
  //mpz_t rel[r + 1];
  //mpz_t relm[r + 1];
  mpz_t *rel = pbc_malloc(sizeof(mpz_t) * (r + 1));
  mpz_t *relm = pbc_malloc(sizeof(mpz_t) * (r + 1));
  unsigned int *prime = pbc_malloc(sizeof(unsigned int) * r);
  darray_t matrix;
  int faci;
  mpz_t k;
  int minfound[fac->count];

  for (i=0; i<r+1; i++) mpz_init(rel[i]);
  for (i=0; i<r+1; i++) mpz_init(relm[i]);

  mpz_init(k);
  mpz_init(z);
  mpz_init(z1);
  mpz_init(z0);

  darray_init(matrix);

  for (i=0; i<fac->count; i++) {
    darray_append(matrix, pbc_malloc(r * sizeof(mpz_t *)));
    minfound[i] = 0;
  }

  for (j=0; j<fac->count; j++) {
    mpz_t **row = matrix->item[j];
    for (i=0; i<r; i++) row[i] = NULL;
  }

  printf("building prime table...\n");
  prime[0] = 2;
  mpz_set_ui(z, 2);
  for (i=1; i<r; i++) {
    mpz_nextprime(z, z);
    prime[i] = mpz_get_ui(z);
  }
  printf("searching for r-smooth numbers\n");

  mpz_set_ui(z, 1);
  mpz_init(k);
  int try = 0;
  do {
    mpz_mul(z, z, g);
    mpz_mod(z, z, q);

    mpz_add_ui(k, k, 1);
    /*
    pbc_mpz_random(k, q);
    mpz_powm(z, g, k, q);
    */

    try++;

    mpz_set(z1, z);
    for (i=0; i<r; i++) {
      mpz_set_ui(rel[i], 0);
      while (mpz_divisible_ui_p(z1, prime[i])) {
        mpz_add_ui(rel[i], rel[i], 1);
        mpz_divexact_ui(z1, z1, prime[i]);
      }
    }
    if (mpz_cmp_ui(z1, 1)) {
      continue;
    }
    mpz_set(rel[r], k);

/*
    printf("found r-smooth number after %d tries\n", try);
    gmp_printf("g^%Zd = %Zd:", rel[r], z);
    for (i=0; i<r; i++) {
      if (mpz_sgn(rel[i])) {
        gmp_printf(" %u:%Zd", i, rel[i]);
      }
    }
    printf("\n");
*/

    try = 0;

    for (faci=0; faci<fac->count; faci++) {
      mpz_t **row = matrix->item[faci];
      mpz_ptr order = fac->item[faci];
      //gmp_printf("mod %Zd\n", order);
      for (i=0; i<r+1; i++) {
        mpz_mod(relm[i], rel[i], order);
      }

      for (;;) {
        /*
        for (i=0; i<r && !mpz_sgn(relm[i]); i++);
        if (i == r) {
          //printf("redundant relation\n");
          break;
        }
        */
        for (i=r-1; i>=0 && !mpz_sgn(relm[i]); i--);
        if (i < 0) {
          //printf("redundant relation\n");
          break;
        }
        if (i < minfound[faci]) {
          break;
        }
        mpz_set(z0, relm[i]);
        if (!row[i]) {
          row[i] = pbc_malloc(sizeof(mpz_t) * (r + 1));
          mpz_invert(z1, z0, order);
          for (j=0; j<r+1; j++) {
            mpz_init(row[i][j]);
            mpz_mul(row[i][j], z1, relm[j]);
            mpz_mod(row[i][j], row[i][j], order);
          }
          count++;
printf("%d / %d\n", count, r * fac->count);
/*
for (j=0; j<r; j++) {
  if (mpz_sgn(row[i][j])) {
    gmp_printf(" %d:%Zd", j, row[i][j]);
  }
}
gmp_printf(" %Zd\n", row[i][j]);
*/

          if (i == minfound[faci]) {
            do {
              if (!minfound[faci]) {
              printf("found log p_%d\n", minfound[faci]);
              gmp_printf("km = %Zd mod %Zd\n", row[i][r], order);
              }
              minfound[faci]++;
              if (minfound[faci] >= r) break;
            } while (row[minfound[faci]]);
          }
          break;
        }

        /*
        printf("before:");
        for (j=0; j<r; j++) {
          if (mpz_sgn(relm[j])) {
            gmp_printf(" %d:%Zd", j, relm[j]);
          }
        }
        gmp_printf(" %Zd\n", relm[j]);

        printf("sub %d:", i);
        for (j=0; j<r; j++) {
          if (mpz_sgn(row[i][j])) {
            gmp_printf(" %d:%Zd", j, row[i][j]);
          }
        }
        gmp_printf(" %Zd\n", row[i][j]);
        */

        for (j=0; j<r+1; j++) {
          mpz_mul(z1, z0, row[i][j]);
          mpz_sub(relm[j], relm[j], z1);
          mpz_mod(relm[j], relm[j], order);
        }

        /*
        printf("after:");
        for (j=0; j<r; j++) {
          if (mpz_sgn(relm[j])) {
            gmp_printf(" %d:%Zd", j, relm[j]);
          }
        }
        gmp_printf(" %Zd\n", relm[j]);
        */
      }
    }

  } while (count < r * fac->count);

  for (faci=0; faci<fac->count; faci++) {
    mpz_t **row = matrix->item[faci];
    mpz_ptr order = fac->item[faci];
    /*
    gmp_printf("mod %Zd:\n", order);
    for (i=0; i<r; i++) {
      for (j=0; j<r+1; j++) {
        gmp_printf(" %Zd", row[i][j]);
      }
      printf("\n");
    }
    printf("\n");
    */

    for (i=1; i<r; i++) {
      for (j=0; j<i; j++) {
        if (mpz_sgn(row[i][j])) {
          mpz_mul(z0, row[i][j], row[j][r]);
          mpz_sub(row[i][r], row[i][r], z0);
          mpz_mod(row[i][r], row[i][r], order);
        }
      }
    }
    /*
    for (i=r-2; i>=0; i--) {
      for (j=i+1; j<r; j++) {
        if (mpz_sgn(row[i][j])) {
          mpz_mul(z0, row[i][j], row[j][r]);
          mpz_sub(row[i][r], row[i][r], z0);
          mpz_mod(row[i][r], row[i][r], order);
        }
      }
    }
    */

    /*
    for (i=0; i<r; i++) {
      mpz_set(rel[i], row[i][r]);
      gmp_printf(" %Zd", row[i][r]);
      printf("\n");
    }
    */
  }

  mpz_ptr *tmp = pbc_malloc(sizeof(mpz_ptr) * fac->count);
  for (i=0; i<fac->count; i++) {
    tmp[i] = pbc_malloc(sizeof(mpz_t));
    mpz_init(tmp[i]);
    mpz_pow_ui(fac->item[i], fac->item[i], (unsigned int) mul->item[i]);
  }

  for (i=0; i<r; i++) {
    for (faci=0; faci<fac->count; faci++) {
      mpz_t **row = matrix->item[faci];
      mpz_set(tmp[faci], row[i][r]);
    }
    CRT(ind[i], tmp, (mpz_ptr *) fac->item, fac->count);
  }

  for (i=0; i<fac->count; i++) {
    mpz_clear(tmp[i]);
  }
  pbc_free(tmp);

  for (faci=0; faci<matrix->count; faci++) {
    mpz_t **row = matrix->item[faci];
    for (j=0; j<r; j++) {
      for (i=0; i<r+1; i++) {
        mpz_clear(row[j][i]);
      }
      pbc_free(row[j]);
    }
    pbc_free(row);
  }
  darray_clear(matrix);
  for (i=0; i<r+1; i++) mpz_clear(rel[i]);
  for (i=0; i<r+1; i++) mpz_clear(relm[i]);
  pbc_free(prime);
  pbc_free(rel);
  pbc_free(relm);
  mpz_clear(k);
  mpz_clear(z);
  mpz_clear(z0);
  mpz_clear(z1);

  printf("step 1 completed\n");
  for (i=0; i<r; i++) element_printf(" %Zd", ind[i]);
  printf("\n");
}

static void index_calculus_step2(mpz_t x, mpz_t *ind, int r,
    mpz_t g, mpz_t h, mpz_t q) {
  mpz_t prime;
  mpz_t s;
  mpz_t z, z1;
  mpz_t rel[r];
  int i;

  mpz_init(z);
  mpz_init(z1);
  mpz_init(s);
  mpz_init(prime);
  for (i=0; i<r; i++) mpz_init(rel[i]);

  mpz_set(z, h);

  for (;;) {
    mpz_mul(z, z, g);
    mpz_mod(z, z, q);
    mpz_add_ui(s, s, 1);

    mpz_set(z1, z);
    mpz_set_ui(prime, 1);
    for (i=0; i<r; i++) {
      mpz_set_ui(rel[i], 0);
      mpz_nextprime(prime, prime);
      while (mpz_divisible_p(z1, prime)) {
        mpz_add_ui(rel[i], rel[i], 1);
        mpz_divexact(z1, z1, prime);
      }
    }
    if (mpz_cmp_ui(z1, 1)) continue;
    element_printf("found r-smooth number on try #%Zd\n", s);
    mpz_set_ui(x, 0);
    for (i=0; i<r; i++) {
      mpz_mul(z, rel[i], ind[i]);
      mpz_add(x, x, z);
    }
    mpz_sub(x, x, s);
    mpz_sub_ui(z, q, 1);
    mpz_mod(x, x, z);
    break;
  }
}

static void mpzclear(void *p) {
  mpz_clear(p);
  pbc_free(p);
}

struct addfm_scope_var {
  darray_ptr fac, mul;
};

static int addfm(mpz_t f, unsigned int m, struct addfm_scope_var *v) {
  darray_append(v->fac, f);
  darray_append(v->mul, int_to_voidp(m));
  return 0;
}

void pbc_mpz_index_calculus(mpz_t x, mpz_t g, mpz_t h, mpz_t q) {
  int i, r;
  mpz_t q1, z0;

  mpz_init(q1);
  mpz_init(z0);

  mpz_sub_ui(q1, q, 1);
  mpz_setbit(z0, 6);

  darray_t fac, mul;
  darray_init(fac);
  darray_init(mul);
  struct addfm_scope_var v = {.fac = fac, .mul = mul};
  pbc_trial_divide((int(*)(mpz_t,unsigned,void*))addfm, &v, q1, z0);

  for (i=0; i<mul->count; i++) {
    unsigned int m = (unsigned int) mul->item[i];
    if (m != 1) {
      //TODO
      printf("p-adics not implemented yet\n");
      return;
    }
  }

  {
    double dq = mpz_get_d(q);
    //r = exp(sqrt(log(dq)*log(log(dq))));
    //printf("r = %d\n", r);
    r = exp(1.2 * sqrt(log(dq)));
    printf("r = %d\n", r);
  }
  mpz_t *ind = pbc_malloc(sizeof(mpz_t) * r);
  for (i=0; i<r; i++) mpz_init(ind[i]);

  if (is_gen(g, q, fac, mul)) {

    index_calculus_step1(ind, r, g, q, fac, mul);

    index_calculus_step2(x, ind, r, g, h, q);
  } else {
    mpz_t y, z;
    mpz_t d;

    mpz_init(d);
    mpz_init(y);
    mpz_init(z);
    do {
      pbc_mpz_random(z, q);
    } while (!is_gen(z, q, fac, mul));

    gmp_printf("new gen: %Zd\n", z);

    index_calculus_step1(ind, r, z, q, fac, mul);
    //slow_index_calculus_step1(ind, r, z, q, fac, mul);

    index_calculus_step2(x, ind, r, z, g, q);
    index_calculus_step2(y, ind, r, z, h, q);
    //want y / x mod q-1
    mpz_gcd(d, x, q1);
    mpz_divexact(q1, q1, d);
    mpz_divexact(x, x, d);
    //if y not divisible by d there is no solution
    mpz_divexact(y, y, d);
    mpz_invert(x, x, q1);
    mpz_mul(x, y, x);
    mpz_mod(x, x, q1);

    do {
      mpz_powm(z0, g, x, q);
      if (!mpz_cmp(z0, h)) {
        break;
      }
      mpz_add(x, x, q1);
      mpz_sub_ui(d, d, 1);
    } while (mpz_sgn(d));

    mpz_clear(d);
    mpz_clear(y);
    mpz_clear(z);
  }

  for (i=0; i<r; i++) mpz_clear(ind[i]);
  pbc_free(ind);

  darray_forall(fac, mpzclear);
  darray_clear(mul);
  darray_clear(fac);
  mpz_clear(q1);
  mpz_clear(z0);
}