/*
 * Verify Cert
 * version 0.95
 *
 * Copyright (c) 2013 Dana Jacobsen (dana@acm.org).
 * This is free software; you can redistribute it and/or modify it under
 * the same terms as the Perl 5 programming language system itself.
 *
 * Verifies Primo and MPU certificates.
 *
 * Return values:
 *   0  all numbers are verified prime.
 *   1  at least one number was verified composite.
 *   2  the certificate does not provide a complete proof.
 *   3  there is an error in the certificate.
 *
 * TODO: Allow multiple proofs per input file
 * TODO: Projective EC for ~4x faster operation
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <gmp.h>

/*****************************************************************************/
/* Preliminary definitions                                                   */
/*****************************************************************************/

/* Projective doesn't work yet */
#define USE_AFFINE_EC 1

#define RET_PRIME 0
#define RET_COMPOSITE 1
#define RET_INVALID 2
#define RET_ERROR 3

#define CERT_UNKNOWN 0
#define CERT_PRIMO   1
#define CERT_MPU     2

#define MAX_LINE_LEN 60000
#define MAX_STEPS    20000
#define MAX_QARRAY   100
#define BAD_LINES_ALLOWED  5    /* Similar to WraithX's verifier */

typedef unsigned long UV;
typedef   signed long IV;
#define croak(fmt,...)  { gmp_printf(fmt,##__VA_ARGS__); exit(RET_ERROR); }
#define MPUassert(c,text) if (!(c)) { croak("Internal error: " text); }

#define BGCD_PRIMES      168
#define BGCD_LASTPRIME   997
#define BGCD_NEXTPRIME  1009

void GMP_pn_primorial(mpz_t prim, UV n);
UV trial_factor(mpz_t n, UV from_n, UV to_n);
int miller_rabin_ui(mpz_t n, UV base);
int miller_rabin(mpz_t n, mpz_t a);
void lucas_seq(mpz_t U, mpz_t V, mpz_t n, IV P, IV Q, mpz_t k, mpz_t Qk, mpz_t t);
#define mpz_mulmod(r, a, b, n, t)  \
  do { mpz_mul(t, a, b); mpz_mod(r, t, n); } while (0)

/*****************************************************************************/
/* Some global variables and functions we'll use                             */
/*****************************************************************************/

int _verbose = 0;
int _quiet = 0;
int _testcount = 0;
int _base = 10;
int _step = 0;
int _format = CERT_UNKNOWN;
char _line[MAX_LINE_LEN+1];
char _vstr[MAX_LINE_LEN+1];
const char* _filename;
FILE* _fh;
mpz_t PROOFN, N, A, B, M, Q, X, Y, LQ, LP, S, R, T, J, T1, T2;
mpz_t QARRAY[MAX_QARRAY];
mpz_t AARRAY[MAX_QARRAY];
mpz_t _bgcd;

int   _num_chains = 0;
mpz_t _chain_n[MAX_STEPS];
mpz_t _chain_q[MAX_STEPS];

static void var_init(void) {
  int i;
  mpz_init(PROOFN);
  mpz_init(N);
  mpz_init(A);
  mpz_init(B);
  /* MPU: */
  mpz_init(M);
  mpz_init(Q);
  mpz_init(X);
  mpz_init(Y);
  mpz_init(LQ);
  mpz_init(LP);
  /* Primo: */
  mpz_init(S);
  mpz_init(R);
  mpz_init(T);
  mpz_init(J);

  mpz_init(_bgcd);
  GMP_pn_primorial(_bgcd, BGCD_PRIMES);

  mpz_init(T1);
  mpz_init(T2);

  for (i = 0; i < MAX_QARRAY; i++) {
    mpz_init(QARRAY[i]);
    mpz_init(AARRAY[i]);
  }
}
static void var_free(void) {
  int i;
  mpz_clear(PROOFN);
  mpz_clear(N);
  mpz_clear(A);
  mpz_clear(B);
  mpz_clear(M);
  mpz_clear(Q);
  mpz_clear(X);
  mpz_clear(Y);
  mpz_clear(LQ);
  mpz_clear(LP);
  mpz_clear(S);
  mpz_clear(R);
  mpz_clear(T);
  mpz_clear(J);

  mpz_clear(_bgcd);

  mpz_clear(T1);
  mpz_clear(T2);

  for (i = 0; i < MAX_QARRAY; i++) {
    mpz_clear(QARRAY[i]);
    mpz_clear(AARRAY[i]);
  }
}

static void quit_prime(void) {
  if (!_quiet) printf("                                             \r");
  if (!_quiet) printf("PRIME\n");
  var_free();
  exit(RET_PRIME);
}
static void quit_composite(void) {
  if (!_quiet) printf("                                             \r");
  if (!_quiet) printf("COMPOSITE\n");
  var_free();
  exit(RET_COMPOSITE);
}
static void quit_invalid(const char* type, const char* msg) {
  if (!_quiet) printf("\n");
  if (!_quiet) gmp_printf("%s: step %d, %Zd failed condition %s\n", type, _step, N, msg);
  var_free();
  exit(RET_INVALID);
}
static void quit_error(const char* msg1, const char* msg2) {
  if (!_quiet) printf("\n");
  if (!_quiet) gmp_printf("ERROR: step %d, %s%s\n", _step, msg1, msg2);
  var_free();
  exit(RET_ERROR);
}

#if USE_AFFINE_EC
/*****************************************************************************/
/* EC: affine with point (x,y,1)                                             */
/*****************************************************************************/

struct ec_affine_point  { mpz_t x, y; };
/* P3 = P1 + P2 */
static void _ec_add_AB(mpz_t n,
                    struct ec_affine_point P1,
                    struct ec_affine_point P2,
                    struct ec_affine_point *P3,
                    mpz_t m,
                    mpz_t t1,
                    mpz_t t2)
{
  if (!mpz_cmp(P1.x, P2.x)) {
    mpz_add(t2, P1.y, P2.y);
    mpz_mod(t1, t2, n);
    if (!mpz_cmp_ui(t1, 0) ) {
      mpz_set_ui(P3->x, 0);
      mpz_set_ui(P3->y, 1);
      return;
    }
  }

  mpz_sub(t1, P2.x, P1.x);
  mpz_mod(t2, t1, n);

  /* t1 = 1/deltay mod n */
  if (!mpz_invert(t1, t2, n)) {
    /* We've found a factor!  In multiply, gcd(mult,n) will be a factor. */
    mpz_set_ui(P3->x, 0);
    mpz_set_ui(P3->y, 1);
    return;
  }

  mpz_sub(m, P2.y, P1.y);
  mpz_mod(t2, m, n);        /* t2 = deltay   mod n */
  mpz_mul(m, t1, t2);
  mpz_mod(m, m, n);         /* m = deltay / deltax   mod n */

  /* x3 = m^2 - x1 - x2 mod n */
  mpz_mul(t1, m, m);
  mpz_sub(t2, t1, P1.x);
  mpz_sub(t1, t2, P2.x);
  mpz_mod(P3->x, t1, n);
  /* y3 = m(x1 - x3) - y1 mod n */
  mpz_sub(t1, P1.x, P3->x);
  mpz_mul(t2, m, t1);
  mpz_sub(t1, t2, P1.y);
  mpz_mod(P3->y, t1, n);
}

/* P3 = 2*P1 */
static void _ec_add_2A(mpz_t a,
                    mpz_t n,
                    struct ec_affine_point P1,
                    struct ec_affine_point *P3,
                    mpz_t m,
                    mpz_t t1,
                    mpz_t t2)
{
  /* m = (3x1^2 + a) * (2y1)^-1 mod n */
  mpz_mul_ui(t1, P1.y, 2);
  if (!mpz_invert(m, t1, n)) {
    mpz_set_ui(P3->x, 0);
    mpz_set_ui(P3->y, 1);
    return;
  }
  mpz_mul_ui(t1, P1.x, 3);
  mpz_mul(t2, t1, P1.x);
  mpz_add(t1, t2, a);
  mpz_mul(t2, m, t1);
  mpz_tdiv_r(m, t2, n);

  /* x3 = m^2 - 2x1 mod n */
  mpz_mul(t1, m, m);
  mpz_mul_ui(t2, P1.x, 2);
  mpz_sub(t1, t1, t2);
  mpz_tdiv_r(P3->x, t1, n);

  /* y3 = m(x1 - x3) - y1 mod n */
  mpz_sub(t1, P1.x, P3->x);
  mpz_mul(t2, t1, m);
  mpz_sub(t1, t2, P1.y);
  mpz_tdiv_r(P3->y, t1, n);
}

static int ec_affine_multiply(mpz_t a, mpz_t k, mpz_t n, struct ec_affine_point P, struct ec_affine_point *R, mpz_t d)
{
  int found = 0;
  struct ec_affine_point A, B, C;
  mpz_t t, t2, t3, mult;

  mpz_init(A.x); mpz_init(A.y);
  mpz_init(B.x); mpz_init(B.y);
  mpz_init(C.x); mpz_init(C.y);
  mpz_init(t);   mpz_init(t2);   mpz_init(t3);
  mpz_init_set_ui(mult, 1);  /* holds intermediates, gcd at end */

  mpz_set(A.x, P.x);  mpz_set(A.y, P.y);
  mpz_set_ui(B.x, 0); mpz_set_ui(B.y, 1);

  /* Binary ladder multiply. */
  while (mpz_cmp_ui(k, 0) > 0) {
    if (mpz_odd_p(k)) {
      mpz_sub(t, B.x, A.x);
      mpz_mul(t2, mult, t);
      mpz_mod(mult, t2, n);
      if ( !mpz_cmp_ui(A.x, 0) && !mpz_cmp_ui(A.y, 1) ) {
        /* nothing */
      } else if ( !mpz_cmp_ui(B.x, 0) && !mpz_cmp_ui(B.y, 1) ) {
        mpz_set(B.x, A.x);  mpz_set(B.y, A.y);
      } else {
        _ec_add_AB(n, A, B, &C, t, t2, t3);
        /* If the add failed to invert, then we have a factor. */
        mpz_set(B.x, C.x);  mpz_set(B.y, C.y);
      }
      mpz_sub_ui(k, k, 1);
    } else {
      mpz_mul_ui(t, A.y, 2);
      mpz_mul(t2, mult, t);
      mpz_mod(mult, t2, n);

      _ec_add_2A(a, n, A, &C, t, t2, t3);
      mpz_set(A.x, C.x);  mpz_set(A.y, C.y);
      mpz_tdiv_q_2exp(k, k, 1);
    }
  }
  mpz_gcd(d, mult, n);
  found = (mpz_cmp_ui(d, 1) && mpz_cmp(d, n));

  mpz_tdiv_r(R->x, B.x, n);
  mpz_tdiv_r(R->y, B.y, n);

  mpz_clear(mult);
  mpz_clear(t);   mpz_clear(t2);   mpz_clear(t3);
  mpz_clear(A.x); mpz_clear(A.y);
  mpz_clear(B.x); mpz_clear(B.y);
  mpz_clear(C.x); mpz_clear(C.y);

  return found;
}

#else

/*****************************************************************************/
/* EC: projective with point (X,1,Z) (Montgomery)                            */
/*****************************************************************************/

/* (xout:zout) = (x1:z1) + (x2:z2) */
static void pec_add3(mpz_t xout, mpz_t zout,
                     mpz_t x1, mpz_t z1,
                     mpz_t x2, mpz_t z2,
                     mpz_t xin, mpz_t zin,
                     mpz_t n, mpz_t u, mpz_t v, mpz_t w)
{
  mpz_sub(u, x2, z2);
  mpz_add(v, x1, z1);
  mpz_mulmod(u, u, v, n, w);     /* u = (x2 - z2) * (x1 + z1) % n */

  mpz_add(v, x2, z2);
  mpz_sub(w, x1, z1);
  mpz_mulmod(v, v, w, n, v);     /* v = (x2 + z2) * (x1 - z1) % n */

  mpz_add(w, u, v);              /* w = u+v */
  mpz_sub(v, u, v);              /* v = u-v */

  mpz_mulmod(w, w, w, n, u);     /* w = (u+v)^2 % n */
  mpz_mulmod(v, v, v, n, u);     /* v = (u-v)^2 % n */

  mpz_set(u, xin);
  mpz_mulmod(xout, w, zin, n, w);
  mpz_mulmod(zout, v, u,   n, w);
  /* 6 mulmods, 6 adds */
}

/* (x2:z2) = 2(x1:z1) */
static void pec_double(mpz_t x2, mpz_t z2, mpz_t x1, mpz_t z1,
                       mpz_t b, mpz_t n, mpz_t u, mpz_t v, mpz_t w)
{
  mpz_add(u, x1, z1);
  mpz_mulmod(u, u, u, n, w);     /* u = (x1+z1)^2 % n */

  mpz_sub(v, x1, z1);
  mpz_mulmod(v, v, v, n, w);     /* v = (x1-z1)^2 % n */

  mpz_mulmod(x2, u, v, n, w);    /* x2 = uv % n */

  mpz_sub(w, u, v);              /* w = u-v = 4(x1 * z1) */
  mpz_mulmod(u, b, w, n, z2);
  mpz_add(u, u, v);              /* u = (v+b*w) mod n */
  mpz_mulmod(z2, w, u, n, v);    /* z2 = (w*u) mod n */
  /* 5 mulmods, 4 adds */
}

#define NORMALIZE(f, u, v, x, z, n) \
    mpz_gcdext(f, u, NULL, z, n); \
    mpz_mulmod(x, x, u, n, v); \
    mpz_set_ui(z, 1);

static void pec_mult(mpz_t a, mpz_t b, mpz_t k, mpz_t n, mpz_t x, mpz_t z)
{
  mpz_t u, v, w, x1, x2, z1, z2, r;
  int l = -1;

  mpz_init(u); mpz_init(v); mpz_init(w);
  mpz_init(x1);  mpz_init(x2);  mpz_init(z1);  mpz_init(z2);

  mpz_sub_ui(k, k, 1);
  mpz_init_set(r, k);
  while (mpz_cmp_ui(r, 1) > 0) {
    mpz_tdiv_q_2exp(r, r, 1);
    l++;
  }
  mpz_clear(r);

  //gmp_printf("x is %Zd  z is %Zd  k is %Zd\n", x, z, k);
  if (mpz_tstbit(k, l)) {
    pec_double(x2, z2, x, z, b, n, u, v, w);
    pec_add3(x1, z1, x2, z2, x, z, x, z, n, u, v, w);
    pec_double(x2, z2, x2, z2, b, n, u, v, w);
  } else {
    pec_double(x1, z1, x, z, b, n, u, v, w);
    pec_add3(x2, z2, x, z, x1, z1, x, z, n, u, v, w);
  }
  l--;
  while (l >= 1) {
    if (mpz_tstbit(k, l)) {
      pec_add3(x1, z1, x1, z1, x2, z2, x, z, n, u, v, w);
      pec_double(x2, z2, x2, z2, b, n, u, v, w);
    } else {
      pec_add3(x2, z2, x2, z2, x1, z1, x, z, n, u, v, w);
      pec_double(x1, z1, x1, z1, b, n, u, v, w);
    }
    l--;
  }
  if (mpz_tstbit(k, 0)) {
    pec_double(x, z, x2, z2, b, n, u, v, w);
  } else {
    pec_add3(x, z, x2, z2, x1, z1, x, z, n, u, v, w);
  }
  mpz_clear(u);  mpz_clear(v);  mpz_clear(w);
  mpz_clear(x1);  mpz_clear(x2);  mpz_clear(z1);  mpz_clear(z2);
}

#endif

/*****************************************************************************/
/* M-R, Lucas, BPSW                                                          */
/*****************************************************************************/
int miller_rabin_ui(mpz_t n, UV base)
{
  int rval;
  mpz_t a;
  mpz_init_set_ui(a, base);
  rval = miller_rabin(n, a);
  mpz_clear(a);
  return rval;
}

int miller_rabin(mpz_t n, mpz_t a)
{
  mpz_t nminus1, d, x;
  UV s, r;
  int rval;

  {
    int cmpr = mpz_cmp_ui(n, 2);
    if (cmpr == 0)     return 1;  /* 2 is prime */
    if (cmpr < 0)      return 0;  /* below 2 is composite */
    if (mpz_even_p(n)) return 0;  /* multiple of 2 is composite */
  }
  if (mpz_cmp_ui(a, 1) <= 0)
    croak("Base %ld is invalid", mpz_get_si(a));
  mpz_init_set(nminus1, n);
  mpz_sub_ui(nminus1, nminus1, 1);
  mpz_init_set(x, a);

  /* Handle large and small bases.  Use x so we don't modify their input a. */
  if (mpz_cmp(x, n) >= 0)
    mpz_mod(x, x, n);
  if ( (mpz_cmp_ui(x, 1) <= 0) || (mpz_cmp(x, nminus1) >= 0) ) {
    mpz_clear(nminus1);
    mpz_clear(x);
    return 1;
  }

  mpz_init_set(d, nminus1);
  s = mpz_scan1(d, 0);
  mpz_tdiv_q_2exp(d, d, s);

  mpz_powm(x, x, d, n);
  mpz_clear(d); /* done with a and d */
  rval = 0;
  if (!mpz_cmp_ui(x, 1) || !mpz_cmp(x, nminus1)) {
    rval = 1;
  } else {
    for (r = 1; r < s; r++) {
      mpz_powm_ui(x, x, 2, n);
      if (!mpz_cmp_ui(x, 1)) {
        break;
      }
      if (!mpz_cmp(x, nminus1)) {
        rval = 1;
        break;
      }
    }
  }
  mpz_clear(nminus1); mpz_clear(x);
  return rval;
}

/* Returns Lucas sequence  U_k mod n and V_k mod n  defined by P,Q */
void lucas_seq(mpz_t U, mpz_t V, mpz_t n, IV P, IV Q, mpz_t k,
               mpz_t Qk, mpz_t t)
{
  UV b = mpz_sizeinbase(k, 2);
  IV D = P*P - 4*Q;

  MPUassert( mpz_cmp_ui(n, 2) >= 0, "lucas_seq: n is less than 2" );
  MPUassert( mpz_cmp_ui(k, 0) >= 0, "lucas_seq: k is negative" );
  MPUassert( P >= 0 && mpz_cmp_si(n, P) >= 0, "lucas_seq: P is out of range" );
  MPUassert( mpz_cmp_si(n, Q) >= 0, "lucas_seq: Q is out of range" );
  MPUassert( D != 0, "lucas_seq: D is zero" );

  if (mpz_cmp_ui(k, 0) <= 0) {
    mpz_set_ui(U, 0);
    mpz_set_ui(V, 2);
    return;
  }
  mpz_set_ui(U, 1);
  mpz_set_si(V, P);
  mpz_set_si(Qk, Q);

  if (Q == 1) {
    /* Use the fast V method if possible.  Much faster with small n. */
    mpz_set_si(t, P*P-4);
    if (P > 2 && mpz_invert(t, t, n)) {
      /* Compute V_k and V_{k+1}, then computer U_k from them. */
      mpz_set_si(V, P);
      mpz_init_set_si(U, P*P-2);
      while (b > 1) {
        b--;
        if (mpz_tstbit(k, b-1)) {
          mpz_mul(V, V, U);  mpz_sub_ui(V, V, P);  mpz_mod(V, V, n);
          mpz_mul(U, U, U);  mpz_sub_ui(U, U, 2);  mpz_mod(U, U, n);
        } else {
          mpz_mul(U, V, U);  mpz_sub_ui(U, U, P);  mpz_mod(U, U, n);
          mpz_mul(V, V, V);  mpz_sub_ui(V, V, 2);  mpz_mod(V, V, n);
        }
      }
      mpz_mul_ui(U, U, 2);
      mpz_submul_ui(U, V, P);
      mpz_mul(U, U, t);
    } else {
      /* Fast computation of U_k and V_k, specific to Q = 1 */
      while (b > 1) {
        mpz_mulmod(U, U, V, n, t);     /* U2k = Uk * Vk */
        mpz_mul(V, V, V);
        mpz_sub_ui(V, V, 2);
        mpz_mod(V, V, n);               /* V2k = Vk^2 - 2 Q^k */
        b--;
        if (mpz_tstbit(k, b-1)) {
          mpz_mul_si(t, U, D);
                                      /* U:  U2k+1 = (P*U2k + V2k)/2 */
          mpz_mul_si(U, U, P);
          mpz_add(U, U, V);
          if (mpz_odd_p(U)) mpz_add(U, U, n);
          mpz_fdiv_q_2exp(U, U, 1);
                                      /* V:  V2k+1 = (D*U2k + P*V2k)/2 */
          mpz_mul_si(V, V, P);
          mpz_add(V, V, t);
          if (mpz_odd_p(V)) mpz_add(V, V, n);
          mpz_fdiv_q_2exp(V, V, 1);
        }
      }
    }
  } else {
    while (b > 1) {
      mpz_mulmod(U, U, V, n, t);     /* U2k = Uk * Vk */
      mpz_mul(V, V, V);
      mpz_submul_ui(V, Qk, 2);
      mpz_mod(V, V, n);               /* V2k = Vk^2 - 2 Q^k */
      mpz_mul(Qk, Qk, Qk);            /* Q2k = Qk^2 */
      b--;
      if (mpz_tstbit(k, b-1)) {
        mpz_mul_si(t, U, D);
                                    /* U:  U2k+1 = (P*U2k + V2k)/2 */
        mpz_mul_si(U, U, P);
        mpz_add(U, U, V);
        if (mpz_odd_p(U)) mpz_add(U, U, n);
        mpz_fdiv_q_2exp(U, U, 1);
                                    /* V:  V2k+1 = (D*U2k + P*V2k)/2 */
        mpz_mul_si(V, V, P);
        mpz_add(V, V, t);
        if (mpz_odd_p(V)) mpz_add(V, V, n);
        mpz_fdiv_q_2exp(V, V, 1);

        mpz_mul_si(Qk, Qk, Q);
      }
      mpz_mod(Qk, Qk, n);
    }
  }
  mpz_mod(U, U, n);
  mpz_mod(V, V, n);
}

static int lucas_selfridge_params(IV* P, IV* Q, mpz_t n, mpz_t t)
{
  IV D = 5;
  UV Dui = (UV) D;
  while (1) {
    UV gcd = mpz_gcd_ui(NULL, n, Dui);
    if ((gcd > 1) && mpz_cmp_ui(n, gcd) != 0)
      return 0;
    mpz_set_si(t, D);
    if (mpz_jacobi(t, n) == -1)
      break;
    if (Dui == 21 && mpz_perfect_square_p(n))
      return 0;
    Dui += 2;
    D = (D > 0)  ?  -Dui  :  Dui;
    if (Dui > 1000000)
      croak("lucas_selfridge_params: D exceeded 1e6");
  }
  if (P) *P = 1;
  if (Q) *Q = (1 - D) / 4;
  return 1;
}

static int lucas_extrastrong_params(IV* P, IV* Q, mpz_t n, mpz_t t, UV inc)
{
  UV tP = 3;
  if (inc < 1 || inc > 256)
    croak("Invalid lucas parameter increment: %lu\n", (unsigned long)inc);
  while (1) {
    UV D = tP*tP - 4;
    UV gcd = mpz_gcd_ui(NULL, n, D);
    if (gcd > 1 && mpz_cmp_ui(n, gcd) != 0)
      return 0;
    mpz_set_ui(t, D);
    if (mpz_jacobi(t, n) == -1)
      break;
    if (tP == (3+20*inc) && mpz_perfect_square_p(n))
      return 0;
    tP += inc;
    if (tP > 65535)
      croak("lucas_extrastrong_params: P exceeded 65535");
  }
  if (P) *P = (IV)tP;
  if (Q) *Q = 1;
  return 1;
}
int is_lucas_pseudoprime(mpz_t n, int strength)
{
  mpz_t d, U, V, Qk, t;
  IV P, Q;
  UV s = 0;
  int rval;

  {
    int cmpr = mpz_cmp_ui(n, 2);
    if (cmpr == 0)     return 1;  /* 2 is prime */
    if (cmpr < 0)      return 0;  /* below 2 is composite */
    if (mpz_even_p(n)) return 0;  /* multiple of 2 is composite */
  }

  mpz_init(t);
  rval = (strength < 2) ? lucas_selfridge_params(&P, &Q, n, t)
                        : lucas_extrastrong_params(&P, &Q, n, t, 1);
  if (!rval) {
    mpz_clear(t);
    return 0;
  }

  mpz_init(U);  mpz_init(V);  mpz_init(Qk);
  mpz_init_set(d, n);
  mpz_add_ui(d, d, 1);

  if (strength > 0) {
    s = mpz_scan1(d, 0);
    mpz_tdiv_q_2exp(d, d, s);
  }

  lucas_seq(U, V, n, P, Q, d, Qk, t);
  mpz_clear(d);

  rval = 0;
  if (strength == 0) {
    /* Standard checks U_{n+1} = 0 mod n. */
    rval = (mpz_sgn(U) == 0);
  } else if (strength == 1) {
    if (mpz_sgn(U) == 0) {
      rval = 1;
    } else {
      while (s--) {
        if (mpz_sgn(V) == 0) {
          rval = 1;
          break;
        }
        if (s) {
          mpz_mul(V, V, V);
          mpz_submul_ui(V, Qk, 2);
          mpz_mod(V, V, n);
          mpz_mulmod(Qk, Qk, Qk, n, t);
        }
      }
    }
  } else {
    mpz_sub_ui(t, n, 2);
    if ( mpz_sgn(U) == 0 && (mpz_cmp_ui(V, 2) == 0 || mpz_cmp(V, t) == 0) ) {
      rval = 1;
    } else {
      s--;  /* The extra strong test tests r < s-1 instead of r < s */
      while (s--) {
        if (mpz_sgn(V) == 0) {
          rval = 1;
          break;
        }
        if (s) {
          mpz_mul(V, V, V);
          mpz_sub_ui(V, V, 2);
          mpz_mod(V, V, n);
        }
      }
    }
  }
  mpz_clear(Qk); mpz_clear(V); mpz_clear(U); mpz_clear(t);
  return rval;
}

int is_prob_prime(mpz_t n)
{
  /*  Step 1: Look for small divisors.  This is done purely for performance.
   *          It is *not* a requirement for the BPSW test. */

  /* If less than 1009, make trial factor handle it. */
  if (mpz_cmp_ui(n, BGCD_NEXTPRIME) < 0)
    return trial_factor(n, 2, BGCD_LASTPRIME) ? 0 : 2;

  /* Check for tiny divisors (GMP can do these really fast) */
  if ( mpz_even_p(n)
    || mpz_divisible_ui_p(n, 3)
    || mpz_divisible_ui_p(n, 5) ) return 0;
  /* Do a big GCD with all primes < 1009 */
  {
    mpz_t t;
    mpz_init(t);
    mpz_gcd(t, n, _bgcd);
    if (mpz_cmp_ui(t, 1) != 0) { mpz_clear(t); return 0; }
    mpz_clear(t);
  }
  /* No divisors under 1009 */
  if (mpz_cmp_ui(n, BGCD_NEXTPRIME*BGCD_NEXTPRIME) < 0)
    return 2;

  /*  Step 2: The BPSW test.  psp base 2 and slpsp. */

  /* Miller Rabin with base 2 */
  if (miller_rabin_ui(n, 2) == 0)
    return 0;

  /* Extra-Strong Lucas test */
  if (is_lucas_pseudoprime(n, 2 /*extra strong*/) == 0)
    return 0;
  /* BPSW is deterministic below 2^64 */
  if (mpz_sizeinbase(n, 2) <= 64)
    return 2;

  return 1;
}

/* These primorial and trial factor functions are really slow for numerous
 * reasons, but most of all because mpz_nextprime is dog slow.  We don't
 * really use them, so don't worry about it too much. */
void GMP_pn_primorial(mpz_t prim, UV n)
{
  mpz_t p;
  mpz_init_set_ui(p, 2);
  mpz_set_ui(prim, 1);
  while (n--) {
    mpz_mul(prim, prim, p);
    mpz_nextprime(p, p);
  }
  mpz_clear(p);
}
UV trial_factor(mpz_t n, UV from_n, UV to_n)
{
  mpz_t p;
  UV f = 0;

  if (mpz_cmp_ui(n, 4) < 0)
    return (mpz_cmp_ui(n, 1) <= 0) ? 1 : 0;   /* 0,1 => 1   2,3 => 0 */
  
  if      (from_n <= 2 && to_n >= 2 && mpz_even_p(n)           ) return 2;
  else if (from_n <= 3 && to_n >= 3 && mpz_divisible_ui_p(n, 3)) return 2;
  if (from_n < 5)
    from_n = 5;
  if (from_n > to_n)
    return 0;

  mpz_init(p);
  mpz_sqrt(p, n);
  if (mpz_cmp_ui(p, to_n) < 0)
    to_n = mpz_get_ui(p);      /* limit to_n to sqrtn */
  mpz_set_ui(p, from_n-1);
  mpz_nextprime(p, p);         /* Set p to the first prime >= from_n */

  while (mpz_cmp_ui(p, to_n) <= 0) {
    if (mpz_divisible_p(n, p)) {
      f = mpz_get_ui(p);
      break;
    }
    mpz_nextprime(p, p);
  }
  mpz_clear(p);
  return f;
}

/*****************************************************************************/
/* Proof verification                                                        */
/*****************************************************************************/

/* What each of these does is verify:
 *      Assume Q is prime.
 *      Then N is prime based on the proof given.
 * We verify any necessary conditions on Q (e.g. it must be odd, or > 0, etc.
 * but do not verify Q prime.  That is done in another proof step.
 */

/* ECPP using N, A, B, M, Q, X, Y
 *
 * A.O.L. Atkin and F. Morain, "Elliptic Curves and primality proving"
 * Mathematics of Computation, v61, 1993, pages 29-68.
 * http://www.ams.org/journals/mcom/1993-61-203/S0025-5718-1993-1199989-X/
 *
 * Page 10, Theorem 5.2:
 *  "Let N be an integer prime to 6, E an elliptic curve over Z/NZ, together
 *   with a point P on E and m and s two integers with s | m.  For each prime
 *   divisor q of s, we put (m/q)P = (x_q : y_q : z_q).  We assume that
 *   mP = O_E and gcd(z_q,N) = 1 for all q.  Then, if p is a prime divisor
 *   of N, one has #E(Z/pZ) = 0 mod s."
 * Page 10, Corollary 5.1:
 *  "With the same conditions, if s > (N^(1/4) + 1)^2, then N is prime."
 *
 * Basically this same result is repeated in Crandall and Pomerance 2005,
 * Theorem 7.6.1 "Goldwasser-Kilian ECPP theorem".
 *
 * Wikipedia, "Elliptic curve primality testing":
 *  "Let N be a positive integer, and E be the set which is defined by the
 *   equation y^2 = x^3 + ax + b (mod N).  Consider E over Z/NZ, use the
 *   usual addition law on E, and write O for the neutral element on E.
 *   Let m be an integer.  If there is a prime q which divides m, and is
 *   greater than (N^(1/4) + 1)^2 and there exists a point P on E such that
 *   (1) mP = O, (2) (m/q)P is defined and not equal to O, then N is prime."
 *
 * We use the restricted form as stated by Wikipedia and used in the
 * Atkin/Morain ECPP algorithm, where s is a prime (hence the "for each prime
 * divisor q of s" of the general theorem is just s).
 */
void verify_ecpp(void) {
  mpz_mod(A, A, N);
  mpz_mod(B, B, N);

  if (mpz_cmp_ui(N, 0) <= 0)   quit_invalid("ECPP", "N > 0");
  if (mpz_gcd_ui(NULL, N, 6) != 1)  quit_invalid("ECPP", "gcd(N, 6) = 1");
  mpz_mul(T1, A, A);
  mpz_mul(T1, T1, A);
  mpz_mul_ui(T1, T1, 4);
  mpz_mul(T2, B, B);
  mpz_mul_ui(T2, T2, 27);
  mpz_add(T1, T1, T2);
  mpz_gcd(T1, T1, N);
  if (mpz_cmp_ui(T1, 1) != 0)  quit_invalid("ECPP", "gcd(4*a^3 + 27*b^2, N) = 1");
  mpz_mul(T1, X, X);
  mpz_add(T1, T1, A);
  mpz_mul(T1, T1, X);
  mpz_add(T1, T1, B);
  mpz_mod(T1, T1, N);
  mpz_mul(T2, Y, Y);
  mpz_mod(T2, T2, N);
  if (mpz_cmp(T1, T2) != 0)   quit_invalid("ECPP", "Y^2 = X^3 + A*X + B mod N");
  mpz_mul_ui(T2, N, 4);
  mpz_sqrt(T2, T2);
  mpz_add_ui(T1, N, 1);
  mpz_sub(T1, T1, T2);
  if (mpz_cmp(M, T1) < 0)     quit_invalid("ECPP", "M >= N + 1 - 2*sqrt(N)");
  mpz_add_ui(T1, N, 1);
  mpz_add(T1, T1, T2);
  if (mpz_cmp(M, T1) > 0)     quit_invalid("ECPP", "M <= N + 1 + 2*sqrt(N)");
  mpz_root(T1, N, 4);
  mpz_add_ui(T1, T1, 1);
  mpz_mul(T1, T1, T1);
  if (mpz_cmp(Q, T1) <= 0)    quit_invalid("ECPP", "Q > (N^(1/4)+1)^2");
  if (mpz_cmp(Q, N) >= 0)     quit_invalid("ECPP", "Q < N");
  /* While M = Q is odd to compute in a proof, it is allowed.
   * In Primo terms, this means S=1 is allowed.
   *   if (mpz_cmp(M, Q) == 0)     quit_invalid("ECPP", "M != Q");
   */
  if (!mpz_divisible_p(M, Q)) quit_invalid("ECPP", "Q divides M");

  {
#if USE_AFFINE_EC
    struct ec_affine_point P0, P1, P2;
    mpz_init_set(P0.x, X);  mpz_init_set(P0.y, Y);
    mpz_init(P1.x); mpz_init(P1.y);
    mpz_init(P2.x); mpz_init(P2.y);
    mpz_divexact(T1, M, Q);
    if (ec_affine_multiply(A, T1, N, P0, &P2, T2))
      quit_invalid("ECPP", "Factor found for N");
    /* Check that P2 is not (0,1) */
    if (mpz_cmp_ui(P2.x, 0) == 0 && mpz_cmp_ui(P2.y, 1) == 0)
      quit_invalid("ECPP", "(M/Q) * EC(A,B,N,X,Y) is not identity");
    mpz_set(T1, Q);
    if (ec_affine_multiply(A, T1, N, P2, &P1, T2))
      quit_invalid("ECPP", "Factor found for N");
    /* Check that P1 is (0, 1) */
    if (! (mpz_cmp_ui(P1.x, 0) == 0 && mpz_cmp_ui(P1.y, 1) == 0) )
      quit_invalid("ECPP", "M * EC(A,B,N,X,Y) is identity");
    mpz_clear(P0.x); mpz_clear(P0.y);
    mpz_clear(P1.x); mpz_clear(P1.y);
    mpz_clear(P2.x); mpz_clear(P2.y);
#else
    mpz_t PX, PY, PA, PB;
    mpz_init(PX);  mpz_init(PY);  mpz_init(PA);  mpz_init(PB);

    /* We have A,B,X,Y in affine coordinates, for the curve:
     *   Y^2 = X^3 + AX + B
     * and want to turn this into points on a Montgomery curve:
     *   by^2 = x^3 + ax^2 + x
     * so we can use the much faster (~4x) multiplication routines.
     * The inverse of this operation is:
     *    X = (3x+a)/3b
     *    Y = y/b
     *    A = (3-a^2)/(3b^2)
     *    B = (2a^3-9a)/27b^3
     * In our case we need to do the harder job of going the other direction.
     */

    /* Make Montgomery variables from affine (TODO: make this work) */
    mpz_add(PB, X, A);
    mpz_mul(PB, PB, X);
    mpz_add_ui(PB, PB, 1);
    mpz_mul(PB, PB, X);
    mpz_mod(PB, PB, N);

    mpz_mul_ui(T2, PB, 3);
    mpz_mul(T2, T2, PB);
    mpz_mod(T2, T2, N);
    mpz_gcdext(T2, T1, NULL, T2, N); /* T1 = 1/3g^2 */
    if (mpz_cmp_ui(T2,1) != 0) quit_invalid("ECPP", "Factor found during gcd");

    mpz_mul_ui(PX, X, 3);
    mpz_add(PX, PX, A);
    mpz_mul(PX, PX, PB);
    mpz_mul(PX, PX, T1);
    mpz_mod(PX, PX, N);

    mpz_set(PY, Y);
    mpz_mul_ui(PY, PY, 3);
    mpz_mul(PY, PY, PB);
    mpz_mul(PY, PY, T1);
    mpz_mod(PY, PY, N);  /* y = (3gY)/(3g^2) = Y/g */

    mpz_mul(PA, A, A);
    mpz_sub_ui(PA, PA, 3);
    mpz_neg(PA, PA);
    mpz_mul(PA, PA, T1);
    mpz_mod(PA, PA, N);

    //gmp_printf("PX: %Zd  PY: %Zd\n", PX, PY);

    mpz_divexact(T1, M, Q);
    pec_mult(PA, PB, T1, N, PX, PY);
    //gmp_printf("PX: %Zd  PY: %Zd\n", PX, PY);
    /* Check that point is not (0,0) */
    if (mpz_cmp_ui(PX, 0) == 0 && mpz_cmp_ui(PY, 0) == 0)
      quit_invalid("ECPP", "(M/Q) * EC(A,B,N,X,Y) is not identity");
    mpz_set(T1, Q);
    pec_mult(PA, PB, T1, N, PX, PY);
    //gmp_printf("PX: %Zd  PY: %Zd\n", PX, PY);
    /* Check that point is (0, 0) */
    if (! (mpz_cmp_ui(PX, 0) == 0 && mpz_cmp_ui(PY, 0) == 0) )
      quit_invalid("ECPP", "M * EC(A,B,N,X,Y) is identity");
    mpz_clear(PX);  mpz_clear(PY);  mpz_clear(PA);  mpz_clear(PB);
#endif
  }
}

/* Basic N+1 using N, Q, LP, LQ
 *
 * John Brillhart, D.H. Lehmer, J.L. Selfridge,
 *   "New Primality Criteria and Factorizations of 2^m +/- 1"
 * Mathematics of Computation, v29, n130, April 1975, pp 620-647.
 * http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
 *
 * Page 631, Theorem 15:
 *  "Let N+1 = mq, where q is an odd prime such that 2q-1 > sqrt(N).
 *   If there exists a Lucas sequence {V_k} of discriminant D with
 *   (D|N) = -1 for which N|V_{(N+1)/2}, but N∤V_{m/2}, then N is prime."
 */
void verify_bls15(void) {
  if (mpz_even_p(Q))            quit_invalid("BLS15", "Q odd");
  if (mpz_cmp_ui(Q, 2) <= 0)    quit_invalid("BLS15", "Q > 2");
  mpz_add_ui(T2, N, 1);
  if (!mpz_divisible_p(T2, Q))  quit_invalid("BLS15", "Q divides N+1");
  mpz_divexact(M, T2, Q);
  mpz_mul(T1, M, Q);
  mpz_sub_ui(T1, T1, 1);
  if (mpz_cmp(T1, N) != 0)      quit_invalid("BLS15", "MQ-1 = N");
  if (mpz_cmp_ui(M, 0) <= 0)    quit_invalid("BLS15", "M > 0");
  mpz_mul_ui(T1, Q, 2);
  mpz_sub_ui(T1, T1, 1);
  mpz_sqrt(T2, N);
  if (mpz_cmp(T1, T2) <= 0)     quit_invalid("BLS15", "2Q-1 > sqrt(N)");
  mpz_mul(T1, LP, LP);
  mpz_mul_ui(T2, LQ, 4);
  mpz_sub(T1, T1, T2);
  if (mpz_sgn(T1) == 0)         quit_invalid("BLS15", "D != 0");
  if (mpz_jacobi(T1, N) != -1)  quit_invalid("BLS15", "jacobi(D,N) = -1");
  {
    mpz_t U, V, k;
    IV iLP, iLQ;
    mpz_init(U);  mpz_init(V);  mpz_init(k);

    iLP = mpz_get_si(LP);
    iLQ = mpz_get_si(LQ);
    if (mpz_cmp_si(LP, iLP) != 0)  quit_error("BLS15 LP out of range", "");
    if (mpz_cmp_si(LQ, iLQ) != 0)  quit_error("BLS15 LQ out of range", "");

    mpz_tdiv_q_2exp(k, M, 1);
    lucas_seq(U, V, N, iLP, iLQ, k, T1, T2);
    if (mpz_sgn(V) == 0)        quit_invalid("BLS15", "V_{m/2} mod N != 0");
    mpz_add_ui(k, N, 1);
    mpz_tdiv_q_2exp(k, k, 1);
    lucas_seq(U, V, N, iLP, iLQ, k, T1, T2);
    if (mpz_sgn(V) != 0)        quit_invalid("BLS15", "V_{(N+1)/2} mod N == 0");

    mpz_clear(U);  mpz_clear(V);  mpz_clear(k);
  }
}

/* Simplistic N-1 using N, Q, A
 *
 * Hans Riesel, "Prime Numbers and Computer Methods for Factorization"
 * page 103-104, Theorem 4.6:
 *  "Suppose N-1 = R*F = R prod(j=1,n,q_j^{B_j}), with all q_j's distinct
 *   primes, with GCD(R,F) = 1 and R < F.  If an integer a can be found, s.t.
 *     GCD(A^((N-1)/q_j)-1,N) = 1 for all j=1..n
 *   and satisfying
 *     a^(N-1) = 1 mod N
 *   then N is a prime."
 *
 * Now make the severe restriction that F must be a single prime q.  This then
 * reduces to the Wikipedia "Pocklington primality test":
 *  "Let N > 1 be an integer, and suppose there exist numbers a and q such that
 *   (1) q is prime, q|N-1 and q > sqrt(N)-1
 *   (2) a^(N-1) = 1 mod N
 *   (3) gcd(a^((N-1)/q)-1,N) = 1
 *   Then N is prime."
 *
 * Note that BLS75 theorem 3 is similar, but also allows a smaller q.  Also,
 * BLS75 theorem 5 is a much better method than generalized Pocklington.
 */
void verify_pocklington(void)
{
  mpz_sub_ui(T2, N, 1);
  if (!mpz_divisible_p(T2, Q))  quit_invalid("Pocklington", "Q divides N-1");
  mpz_divexact(M, T2, Q);
  if (mpz_odd_p(M))             quit_invalid("Pocklington", "M is even");
  if (mpz_cmp_ui(M, 0) <= 0)    quit_invalid("Pocklington", "M > 0");
  if (mpz_cmp(M, Q) >= 0)       quit_invalid("Pocklington", "M < Q");
  mpz_mul(T1, M, Q);
  mpz_add_ui(T1, T1, 1);
  if (mpz_cmp(T1, N) != 0)      quit_invalid("Pocklington", "MQ+1 = N");
  if (mpz_cmp_ui(A, 1) <= 0)    quit_invalid("Pocklington", "A > 1");
  mpz_powm(T1, A, T2, N);
  if (mpz_cmp_ui(T1, 1) != 0)   quit_invalid("Pocklington", "A^(N-1) mod N = 1");
  mpz_powm(T1, A, M, N);
  if (mpz_sgn(T1)) mpz_sub_ui(T1, T1, 1);
  else             mpz_set(T1, T2);
  mpz_gcd(T1, T1, N);
  if (mpz_cmp_ui(T1, 1) != 0)   quit_invalid("Pocklington", "gcd(A^M - 1, N) = 1");
}

/* Basic N-1 using N, Q, A
 *
 * John Brillhart, D.H. Lehmer, J.L. Selfridge,
 *   "New Primality Criteria and Factorizations of 2^m +/- 1"
 * Mathematics of Computation, v29, n130, April 1975, pp 620-647.
 * http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
 *
 * Page 622-623, Theorem 3:
 *  "Let N-1 = mp, where p is an odd prime such that 2p+1 > sqrt(N).
 *   If there exists an a for which a^((N-1)/2) = -1 mod N,
 *   but a^(m/2) != -1 mod N, then N is prime."
 */
void verify_bls3(void) {
  if (mpz_even_p(Q))            quit_invalid("BLS3", "Q odd");
  if (mpz_cmp_ui(Q, 2) <= 0)    quit_invalid("BLS3", "Q > 2");
  mpz_sub_ui(T2, N, 1);
  if (!mpz_divisible_p(T2, Q))  quit_invalid("BLS3", "Q divides N-1");
  mpz_divexact(M, T2, Q);
  mpz_mul(T1, M, Q);
  mpz_add_ui(T1, T1, 1);
  if (mpz_cmp(T1, N) != 0)      quit_invalid("BLS3", "MQ+1 = N");
  if (mpz_cmp_ui(M, 0) <= 0)    quit_invalid("BLS3", "M > 0");
  mpz_mul_ui(T1, Q, 2);
  mpz_add_ui(T1, T1, 1);
  mpz_sqrt(T2, N);
  if (mpz_cmp(T1, T2) <= 0)     quit_invalid("BLS3", "2Q+1 > sqrt(N)");
  mpz_sub_ui(T2, N, 1);
  mpz_divexact_ui(T1, T2, 2);
  mpz_powm(T1, A, T1, N);
  if (mpz_cmp(T1, T2) != 0)     quit_invalid("BLS3", "A^((N-1)/2) = N-1 mod N");
  mpz_divexact_ui(T1, M, 2);
  mpz_powm(T1, A, T1, N);
  if (mpz_cmp(T1, T2) == 0)     quit_invalid("BLS3", "A^(M/2) != N-1 mod N");
}

/* Sophisticated N-1 using N, QARRAY, AARRAY
 *
 * John Brillhart, D.H. Lehmer, J.L. Selfridge,
 *   "New Primality Criteria and Factorizations of 2^m +/- 1"
 * Mathematics of Computation, v29, n130, April 1975, pp 620-647.
 * http://www.ams.org/journals/mcom/1975-29-130/S0025-5718-1975-0384673-1/S0025-5718-1975-0384673-1.pdf
 *
 * Page 621-622: "The expression 'N is a psp base a' will be used for a
 * number N which satisfies the congruence a^(N-1) = 1 mod N, 1 < a < N-1,
 * i.e., N is a "pseudoprime" base a."
 * Page 623: "Throughout the rest of this paper the notation N-1 = F_1 R_1
 * will be used, where F_1 is the even factored portion of N-1, R_1 is > 1,
 * and (F_1,R_1) = 1."
 * Page 623: "(I) For each prime p_i dividing F_1 there exists an a_i such
 * that N is a psp base a_i and (a_i^((N-1)/p_i),N) = 1."
 * Page 624, Theorem 5.
 *  "Assume (I) and let m be >= 1.
 *   When m > 1, assume further that λF_1 + 1 ∤ N for 1 <= λ < m.
 *   If N < (mF_1 + 1) [2(F_1)^2 + (r-m)F_1 + 1],
 *   where r and s are defined by R_1 = (N-1)/F_1 = 2(F_1)s + r, 1 <= r < 2F_1,
 *   then N is prime if and only if s = 0 or r^2 - 8s is not a perfect square.
 *   r != 0 since R_1 is odd."
 *
 * Note that we are using m=1, which is simple for the verifier.
 *
 */
void verify_bls5(int num_qs) {
  int i;
  mpz_t F, R, s, r;

  if (mpz_cmp_ui(N, 2) <= 0)   quit_invalid("BLS5", "N > 2");
  if (mpz_even_p(N))           quit_invalid("BLS5", "N odd");
  mpz_sub_ui(T2, N, 1);
  mpz_init_set_ui(F, 1);
  mpz_init_set(R, T2);
  mpz_init(s);  mpz_init(r);
  for (i = 0; i < num_qs; i++) {
    if (mpz_cmp_ui(QARRAY[i], 1 ) <= 0)  quit_invalid("BLS5", "Q > 1");
    if (mpz_cmp(   QARRAY[i], T2) >= 0)  quit_invalid("BLS5", "Q < N-1");
    if (mpz_cmp_ui(AARRAY[i], 1 ) <= 0)  quit_invalid("BLS5", "A > 1");
    if (mpz_cmp(   AARRAY[i], T2) >= 0)  quit_invalid("BLS5", "A < N-1");
    if (!mpz_divisible_p(T2, QARRAY[i])) quit_invalid("BLS5", "Q divides N-1");
    while (mpz_divisible_p(R, QARRAY[i])) {
      mpz_mul(F, F, QARRAY[i]);
      mpz_divexact(R, R, QARRAY[i]);
    }
  }
  mpz_mul(T1, R, F);
  if (mpz_cmp(T1, T2) != 0)    quit_invalid("BLS5", "R == (N-1)/F");
  if (mpz_odd_p(F))            quit_invalid("BLS5", "F is even");
  mpz_gcd(T1, F, R);
  if (mpz_cmp_ui(T1, 1) != 0)  quit_invalid("BLS5", "gcd(F, R) = 1");
  mpz_mul_ui(T1, F, 2);
  mpz_tdiv_qr(s, r, R, T1);

  mpz_mul_ui(T1, F, 2);   /* T1 = 2*F */
  mpz_sub_ui(T2, r, 1);
  mpz_add(T1, T1, T2);    /* T1 = 2*F + (r-1) */
  mpz_mul(T1, T1, F);     /* T1 = 2*F*F + (r-1)*F */
  mpz_add_ui(T1, T1, 1);  /* T1 = 2*F*F + (r-1)*F + 1 */
  mpz_add_ui(T2, F, 1);
  mpz_mul(T1, T1, T2);    /* T1 = (F+1) * (2*F*F + (r-1)*F + 1) */
  if (mpz_cmp(N, T1) >= 0) quit_invalid("BLS5", "N < (F+1)(2F^2+(r-1)F+1)");
  if (mpz_sgn(s) != 0) {
    mpz_mul(T2, r, r);
    mpz_submul_ui(T2, s, 8);  /* T2 = r*r - 8*s */
    if (mpz_perfect_square_p(T2))  quit_invalid("BLS5", "S=0 OR R^2-8S not a perfect square");
  }
  mpz_clear(F); mpz_clear(R); mpz_clear(s); mpz_clear(r);
  mpz_sub_ui(T2, N, 1);
  for (i = 0; i < num_qs; i++) {
    mpz_powm(T1, AARRAY[i], T2, N);
    if (mpz_cmp_ui(T1, 1) != 0)   quit_invalid("BLS5", "A[i]^(N-1) mod N = 1");
    mpz_divexact(T1, T2, QARRAY[i]);
    mpz_powm(T1, AARRAY[i], T1, N);
    if (mpz_sgn(T1)) mpz_sub_ui(T1, T1, 1);
    else             mpz_set(T1, T2);
    mpz_gcd(T1, T1, N);
    if (mpz_cmp_ui(T1, 1) != 0)   quit_invalid("BLS5", "gcd(A[i]^((N-1)/Q[i]) - 1, N) = 1");
  }
}

/* Most basic N-1 using N, QARRAY, A
 *
 * D.H. Lehmer, "Tests for Primality by the Converse of Fermat's Theorem"
 * Bull. AMS, v33, n3, 1927, pp 327-340.
 * http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.bams/1183492108
 *
 * Page 330, Theorem 2:
 *  "If a^x = 1 mod N for x = N-1, but not for x a quotient of N-1 on
 *   division by any of its prime factors, then N is a prime."
 */
void verify_lucas(int num_qs) {
  int i;
  mpz_sub_ui(T2, N, 1);
  mpz_set(R, T2);
  if (mpz_cmp_ui(A, 1) <= 0)  quit_invalid("Lucas", "A > 1");
  if (mpz_cmp(   A, N) >= 0)  quit_invalid("Lucas", "A < N");

  mpz_powm(T1, A, T2, N);
  if (mpz_cmp_ui(T1, 1) != 0)   quit_invalid("Lucas", "A^(N-1) mod N = 1");

  for (i = 1; i < num_qs; i++) {
    if (mpz_cmp_ui(QARRAY[i], 1 ) <= 0)  quit_invalid("Lucas", "Q > 1");
    if (mpz_cmp(   QARRAY[i], T2) >= 0)  quit_invalid("Lucas", "Q < N-1");
    if (!mpz_divisible_p(T2, QARRAY[i])) quit_invalid("Lucas", "Q divides N-1");

    mpz_divexact(T1, T2, QARRAY[i]);
    mpz_powm(T1, A, T1, N);
    if (mpz_cmp_ui(T1, 1) == 0)   quit_invalid("Lucas", "A^((N-1)/Q[i]) mod N != 1");
    while (mpz_divisible_p(R, QARRAY[i]))
      mpz_divexact(R, R, QARRAY[i]);
  }
  if (mpz_cmp_ui(R, 1) != 0)
    quit_invalid("Lucas", "N-1 has only factors Q[i]");
}

void verify_primo_ecpp(void) {
  /* Calculate a,b,x,y from A, B, T, then verify */
  mpz_mul(T1, T, T);
  mpz_add(T1, T1, A);
  mpz_mul(T1, T1, T);
  mpz_add(T1, T1, B);
  mpz_mod(T1, T1, N);   /* L = T1 = T^3 + AT + B mod N */

  mpz_mul(T2, T1, T1);
  mpz_mul(A, A, T2);
  mpz_mod(A, A, N);     /* a = AL^2 mod N */

  mpz_mul(T2, T2, T1);
  mpz_mul(B, B, T2);
  mpz_mod(B, B, N);     /* b = BL^3 mod N */

  mpz_mul(X, T, T1);
  mpz_mod(X, X, N);     /* x = TL mod N */
  mpz_mul(Y, T1, T1);
  mpz_mod(Y, Y, N);     /* y = L^2 mod N */

  mpz_mul(M, R, S);     /* M = R*S */
  mpz_set(Q, R);        /* Q = R */

  verify_ecpp();  /* N, A, B, M, Q, X, Y */
}

void verify_ecpp4(void) {
  mpz_mul_ui(T1, J, 2);
  if (mpz_cmpabs(T1, N) > 0) quit_invalid("Primo Type 4", "|J| <= N/2");
  if (mpz_cmp_ui(T, 0) < 0)  quit_invalid("Primo Type 4", "T >= 0");
  if (mpz_cmp(T, N) >= 0)    quit_invalid("Primo Type 4", "T < N");

  mpz_set_ui(T2, 1728);
  mpz_sub(T2, T2, J);
  mpz_mul(A, T2, J);
  mpz_mul_ui(A, A, 3);    /* A = 3 J (1728-J)   */

  mpz_mul(T2, T2, T2);
  mpz_mul(B, T2, J);
  mpz_mul_ui(B, B, 2);    /* B = 2 J (1728-J)^2 */

  verify_primo_ecpp();  /* (A,B,T)->(a,b,x,y)  (R,S)->(M,Q)  Verify */
}

void verify_ecpp3(void) {
  /* TODO: The latest primo.html doesn't require A,B to meet these conditions */
  mpz_mul_ui(T1, A, 2);
  mpz_mul_ui(T2, B, 2);
  if (mpz_cmpabs(T1, N) > 0) quit_invalid("Primo Type 3", "|A| <= N/2");
  if (mpz_cmpabs(T2, N) > 0) quit_invalid("Primo Type 3", "|B| <= N/2");
  if (mpz_cmp_ui(T, 0) < 0)  quit_invalid("Primo Type 3", "T >= 0");
  if (mpz_cmp(T, N) >= 0)    quit_invalid("Primo Type 3", "T < N");

  verify_primo_ecpp();  /* (A,B,T)->(a,b,x,y)  (R,S)->(M,Q)  Verify */
}

void verify_primo2(void) {
  mpz_set(LQ, Q);
  mpz_set_ui(LP, mpz_odd_p(LQ) ? 2 : 1);
  mpz_set(Q, R);
  verify_bls15();  /* N, Q, LP, LQ */
}

void verify_primo1(void) {
  mpz_set(Q, R);
  mpz_set(A, B);
  mpz_mul(T1, S, R);
  mpz_add_ui(T1, T1, 1);
  if (mpz_cmp(T1, N) != 0)  quit_invalid("Primo Type 1", "SR+1 = N");
  verify_pocklington();  /* N, Q, A */
}

void verify_small(void) {
  if (mpz_sizeinbase(N, 2) > 64)  quit_invalid("Small", "N <= 2^64");
  if (is_prob_prime(N) != 2)      quit_invalid("Small", "N does not pass BPSW");
}

void add_chain(mpz_t n, mpz_t q) {
  mpz_init_set(_chain_n[_num_chains], n);
  mpz_init_set(_chain_q[_num_chains], q);
  _num_chains++;
}
void free_chains(void) {
  while (_num_chains > 0) {
    _num_chains--;
    mpz_clear(_chain_n[_num_chains]);
    mpz_clear(_chain_q[_num_chains]);
  }
}
void verify_chain(mpz_t n) {
  int i, found;
  if (mpz_sizeinbase(n, 2) <= 64) {
    mpz_set(N, n);
    verify_small();
    return;
  }
  found = 0;
  for (i = 0; i < _num_chains; i++) {
    if (mpz_cmp(n, _chain_n[i]) == 0) {
      found = 1;
      verify_chain(_chain_q[i]);
    }
  }
  mpz_set(N, n);
  if (!found)  quit_invalid("Final", "q value has no proof");
}


void verify_final(void) {
  if (_format == CERT_PRIMO) {
    if (mpz_cmp_ui(N, 18) <= 0)  quit_invalid("Primo Type 0", "N > 18");
    if (mpz_cmp_ui(N, 340000000000000UL) >= 0)  quit_invalid("Primo Type 0", "N < 34 * 10^13");
    if (!miller_rabin_ui(N,  2) || !miller_rabin_ui(N,  3) ||
        !miller_rabin_ui(N,  5) || !miller_rabin_ui(N,  7) ||
        !miller_rabin_ui(N, 11) || !miller_rabin_ui(N, 13) ||
        !miller_rabin_ui(N, 17))
      quit_invalid("Primo Type 0", "N is SPSP(2,3,5,7,11,13,17)");
  } else {
    verify_chain(PROOFN);
    free_chains();
  }
}

/*****************************************************************************/
/* File parsing                                                              */
/*****************************************************************************/

static int get_line(int signaleof) {
  size_t i;
  /* Read in a line */
  if (fgets(_line, MAX_LINE_LEN, _fh) != _line) {
    if (signaleof) return 1;
    else           quit_error("Error reading from file: ", _filename);
  }
  _line[MAX_LINE_LEN] = '\0';
  /* Remove trailing newlines and spaces */
  i = strlen(_line);
  while ( i > 0 && isspace(_line[i-1]) )
    _line[--i] = '\0';
  return 0;
}

#define PROCESS_VAR(v) \
  do { \
    mpz_set_str(v, _vstr, _base); \
    for (i = 0; i < nargs; i++) { \
      if (vlist[i] != 0 && strcmp(vlist[i], #v) == 0) { \
        vfound++; \
        vlist[i] = 0; \
        break; \
      } \
    } \
    if (i >= nargs) \
      quit_error("Unknown variable: ", #v); \
  } while (0)

void read_vars(const char* vars) {
  char* varstring = strdup(vars);
  char* vlist[10];
  char  varname;
  int i;
  int nargs = 0;
  int vfound = 0;
  int bad_lines = 0;

  vlist[0] = strtok(varstring, " ");
  while (vlist[nargs] != 0)
    vlist[++nargs] = strtok(NULL, " ");
  while (vfound < nargs) {
    get_line(0);
    if (strlen(_line) == 0)    /* Skip extrenuous blank lines */
      continue;
    if (_format == CERT_PRIMO) {
      int varnum = 0;  /* Did we read a variable properly this line? */
      if (sscanf(_line, "%c$=%s", &varname, _vstr) == 2) {
        for (i = 0; i < nargs && varnum == 0; i++) {
          if (vlist[i] != 0 && varname == vlist[i][0]) {
            switch (varname) {
              case 'S':  mpz_set_str(S, _vstr, 16);  break;
              case 'R':  mpz_set_str(R, _vstr, 16);  break;
              case 'A':  mpz_set_str(A, _vstr, 16);  break;
              case 'B':  mpz_set_str(B, _vstr, 16);  break;
              case 'Q':  mpz_set_str(Q, _vstr, 16);  break;
              case 'T':  mpz_set_str(T, _vstr, 16);  break;
              case 'J':  mpz_set_str(J, _vstr, 16);  break;
              default: quit_error("Internal error: bad Primo variable type","");
                       break;
            }
            varnum = i+1;
          }
        }
      } else if (sscanf(_line, "[%d]", &i) == 1) {
        quit_error("Variables missing from proof step", "");
      }
      if (varnum != 0) {
        vfound++;             /* We found a variable on the list */
        vlist[varnum-1] = 0;  /* It should only appear once */
        bad_lines = 0;
      } else {
        if (_verbose) { printf("%60s\r", ""); printf("skipping bad line: %s\n", _line); }
        if (bad_lines++ >= BAD_LINES_ALLOWED)
          quit_error("Too many bad lines reading variables", "");
      }
    } else {
      if      (sscanf(_line, "N %s", _vstr) == 1) PROCESS_VAR(N);
      else if (sscanf(_line, "A %s", _vstr) == 1) PROCESS_VAR(A);
      else if (sscanf(_line, "B %s", _vstr) == 1) PROCESS_VAR(B);
      else if (sscanf(_line, "M %s", _vstr) == 1) PROCESS_VAR(M);
      else if (sscanf(_line, "Q %s", _vstr) == 1) PROCESS_VAR(Q);
      else if (sscanf(_line, "X %s", _vstr) == 1) PROCESS_VAR(X);
      else if (sscanf(_line, "Y %s", _vstr) == 1) PROCESS_VAR(Y);
      else if (sscanf(_line, "LQ %s", _vstr) == 1) PROCESS_VAR(LQ);
      else if (sscanf(_line, "LP %s", _vstr) == 1) PROCESS_VAR(LP);
      /* ECPP3 and ECPP4 */
      else if (sscanf(_line, "T %s", _vstr) == 1) PROCESS_VAR(T);
      else if (sscanf(_line, "J %s", _vstr) == 1) PROCESS_VAR(J);
      else if (sscanf(_line, "S %s", _vstr) == 1) PROCESS_VAR(S);
      else if (sscanf(_line, "R %s", _vstr) == 1) PROCESS_VAR(R);
      else
        quit_error("Internal error: bad MPU variable type", "");
    }
  }
  free(varstring);
}

/* TODO:
 * rearrange so we (1) open and read everything up to proof for / candidate.
 * then (2) func that does the proof
 * this should let us call #2 if we hit another proof, so we can do multiple
 * proof in a file.
 */

void parse_top(void)
{
  do {
    if (get_line(1))
      quit_error("Count not find primality certificate indicator", "");
  } while (strstr(_line, "Primality Certificate") == 0);
  if (strcmp(_line, "[PRIMO - Primality Certificate]") == 0)
    _format = CERT_PRIMO;
  else if (strcmp(_line, "[MPU - Primality Certificate]") == 0)
    _format = CERT_MPU;
  else
    quit_error("First line in file is not primality certificate indicator","");

  if (_format == CERT_PRIMO) {
    int items_found = 0;
    while (items_found < 3) {
      get_line(0);
      if (sscanf(_line, "TestCount=%d", &_testcount) == 1) items_found++;
      if (strcmp(_line, "[Candidate]") == 0)  items_found++;
      if (sscanf(_line, "N$=%s", _vstr) == 1) items_found++;
    }
    mpz_set_str(PROOFN, _vstr, 16);
  } else {
    while (1) {
      get_line(0);
      if (_line[0] == '#')  continue;
      if (sscanf(_line, "Base %d", &_base) == 1) continue;
      if (strcmp(_line, "Proof for:") == 0) {
        read_vars("N");
        mpz_set(PROOFN, N);
        break;
      }
    }
  }
  if (!_quiet) {
    printf("%60s\r", "");
    printf("N: ");
    if (_verbose) gmp_printf("%Zd ", PROOFN);
    printf("(%d digits)\n", (int)mpz_sizeinbase(PROOFN, 10));
    printf("Verifying probable prime status.\r");
    fflush(stdout);
  }
  if (is_prob_prime(PROOFN) == 0)
    quit_composite();
  mpz_set(N, PROOFN);
}

void process_file(const char* filename)
{
  _step = 0;
  if (strcmp(filename, "-") == 0)
    _fh = stdin;
  else if ((_fh = fopen(filename, "r")) == NULL)
    quit_error("Unable to open file: ", filename);
  _filename = filename;

  parse_top();

  if (_format == CERT_PRIMO) {
    int type;
    while (1) {
      int rstep;
      get_line(0);
      if (sscanf(_line, "[%d]", &rstep) == 1) {
        if (rstep != _step+1)
          quit_error("Wrong step number found", "");
        _step++;
      }
      if (sscanf(_line, "Type=%d", &type) == 1) {
        if (!_quiet) { printf("%60s\r", ""); printf("Step %3d/%-3d %5d digits  Type %d\r", _step, _testcount, (int)mpz_sizeinbase(N,10), type); fflush(stdout); }
        switch (type) {
          case 4:  read_vars("S R J T");   verify_ecpp4();  break;
          case 3:  read_vars("S R A B T"); verify_ecpp3();  break;
          case 2:  read_vars("S R Q");     verify_primo2(); break;
          case 1:  read_vars("S R B");     verify_primo1(); break;
          case 0:  /* verify_small */ break;
          default: quit_error("Parsing", "Unknown type");   break;
        }
        if (type == 0) break;
        mpz_set(N, R);
      }
    }
  } else {
    char type[MAX_LINE_LEN+1];
    while (1) {
      if (get_line(1)) break;
      if (sscanf(_line, "Type %s", type) == 1) {
        { /* Convert type to upper case */
          char* s = type;
          while (*s != '\0') {
            if (islower(*s))  *s = toupper(*s);
            s++;
          }
        }
        _step++;
        /* TODO: Quick parse of the file to count steps? */
        if (!_quiet) { printf("%60s\r", ""); printf("Step %-4d %5d digits  Type %s\r", _step, (int)mpz_sizeinbase(N,10), type); fflush(stdout); }
        if        (strcmp(type, "ECPP" ) == 0) { read_vars("N A B M Q X Y");
                                                 verify_ecpp();
                                                 add_chain(N, Q);
        } else if (strcmp(type, "ECPP3") == 0) { read_vars("N S R A B T");
                                                 verify_ecpp3();
                                                 add_chain(N, Q);
        } else if (strcmp(type, "ECPP4") == 0) { read_vars("N S R J T");
                                                 verify_ecpp4();
                                                 add_chain(N, Q);
        } else if (strcmp(type, "BLS15") == 0) { read_vars("N Q LP LQ");
                                                 verify_bls15();
                                                 add_chain(N, Q);
        } else if (strcmp(type, "BLS3" ) == 0) { read_vars("N Q A");
                                                 verify_bls3();
                                                 add_chain(N, Q);
        } else if (strcmp(type, "POCKLINGTON") == 0) { read_vars("N Q A");
                                                 verify_pocklington();
                                                 add_chain(N, Q);
        } else if (strcmp(type, "SMALL") == 0) { read_vars("N");
                                                 verify_small();
        } else if (strcmp(type, "BLS5") == 0) {
          int i, index;
          for (index = 0; index < MAX_QARRAY; index++) {
            mpz_set_ui(QARRAY[index], 0);
            mpz_set_ui(AARRAY[index], 2);
          }
          mpz_set_ui(QARRAY[0], 2);
          index = 1;
          while (1) {
            get_line(0);
            if (_line[0] == '-') {
              break;
            } else if (sscanf(_line, "N %s", _vstr) == 1) {
              mpz_set_str(N, _vstr, _base);
            } else if (sscanf(_line, "Q[%d] %s", &i, _vstr) == 2) {
              if (i != index) quit_error("BLS5", "Invalid Q index");
              mpz_set_str(QARRAY[i], _vstr, _base);
              index++;
            } else if (sscanf(_line, "A[%d] %s", &i, _vstr) == 2) {
              if (i < 0 || i > index) quit_error("BLS5", "Invalid A index");
              mpz_set_str(AARRAY[i], _vstr, _base);
            }
          }
          verify_bls5(index);
          for (i = 0; i < index; i++)
            add_chain(N, QARRAY[i]);
        } else if (strcmp(type, "LUCAS") == 0) {
          int i, index;
          for (index = 0; index < MAX_QARRAY; index++)
            mpz_set_ui(QARRAY[index], 0);
          index = 1;
          while (1) {
            get_line(0);
            if        (sscanf(_line, "N %s", _vstr) == 1) {
              mpz_set_str(N, _vstr, _base);
            } else if (sscanf(_line, "Q[%d] %s", &i, _vstr) == 2) {
              if (i != index) quit_error("Lucas", "Invalid Q index");
              mpz_set_str(QARRAY[i], _vstr, _base);
              index++;
            } else if (sscanf(_line, "A %s", _vstr) == 1) {
              mpz_set_str(A, _vstr, _base);
              break;
            }
          }
          verify_lucas(index);
          for (i = 1; i < index; i++)
            add_chain(N, QARRAY[i]);
        } else {
          quit_error("Parsing", "Unknown type");
        }
      }
    }
  }
  verify_final();
}

static void dieusage(const char* prog) {
  printf("Verify Cert version 0.9.  Dana Jacobsen\n\n");
  printf("Usage: %s [options] <file>\n\n", prog);
  printf("Options:\n");
  printf("   -v     set verbose\n");
  printf("   -q     set quiet (no output, only exit code)\n");
  printf("   -help  this message\n");
  var_free();
  exit(RET_INVALID);
}


int main(int argc, char *argv[])
{
  int i;
  int optdone = 0;
  if (argc < 2) dieusage(argv[0]);

  var_init();
  mpz_set_ui(N, 0);

  for (i = 1; i < argc; i++) {
    if (!optdone && argv[i][0] == '-') {
      if (strcmp(argv[i], "--") == 0) {
        optdone = 1;
      } else if (argv[i][1] == '\0') {
        process_file("-");
      } else if (strcmp(argv[i], "-v") == 0) {
        _verbose++;
      } else if (strcmp(argv[i], "-q") == 0) {
        _quiet++;
      } else if (strcmp(argv[i], "-help") == 0 || strcmp(argv[i], "--help") == 0) {
        dieusage(argv[0]);
      } else {
        printf("Unknown option: %s\n\n", argv[i]);
        dieusage(argv[0]);
      }
      continue;
    }
    /* process_file will exit if not verified prime */
    process_file(argv[i]);
  }
  if (mpz_sgn(N) == 0)
    dieusage(argv[0]);
  quit_prime();
  exit(RET_PRIME);
}