random_prime.c

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gmp.h>
#include "ptypes.h"
#include "random_prime.h"
#include "utility.h"
#include "primality.h"
#include "gmp_main.h"
#include "isaac.h"
#include "prime_iterator.h"

static char pr[31] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127};

void mpz_random_nbit_prime(mpz_t p, UV n)
{
  switch (n) {
    case 0:
    case 1:   mpz_set_ui(p, 0);  return;
    case 2:   mpz_set_ui(p, pr[ 0+isaac_rand(2)]);  return;
    case 3:   mpz_set_ui(p, pr[ 2+isaac_rand(2)]);  return;
    case 4:   mpz_set_ui(p, pr[ 4+isaac_rand(2)]);  return;
    case 5:   mpz_set_ui(p, pr[ 6+isaac_rand(5)]);  return;
    case 6:   mpz_set_ui(p, pr[11+isaac_rand(7)]);  return;
    case 7:   mpz_set_ui(p, pr[18+isaac_rand(13)]); return;
    default:  break;
  }
  /* For 32-bit inputs, use fast trivial method */
  if (n <= 32) {
    uint32_t mask = (0xFFFFFFFFU >> (34-n)) << 1,  base = mask+3;
    do {
      mpz_set_ui(p, base | (isaac_rand32() & mask));
    } while (!_GMP_is_prob_prime(p));
  } else {
#if 0
    do {                        /* Trivial method. */
      mpz_isaac_urandomb(p, n);
      mpz_setbit(p, n-1);
      mpz_setbit(p, 0);
    } while (!_GMP_is_prob_prime(p));
#else
    mpz_t base;                 /* Fouque+Tibouchi Alg 1, without modulo checks */
    mpz_init(base);
    if (n > 33) { mpz_isaac_urandomb(base, n-33); mpz_mul_2exp(base,base,1); }
    mpz_setbit(base, n-1);
    mpz_setbit(base, 0);
    do {
      mpz_set_ui(p, isaac_rand32());
      mpz_mul_2exp(p, p, n-32);
      mpz_ior(p, p, base);
    } while (!_GMP_is_prob_prime(p));
    mpz_clear(base);
#endif
  }
}

/* PRIMEINC: pick random value, select next prime. */
/* Fast but bad distribution. */
static int _random_prime_primeinc(mpz_t p, const mpz_t lo, const mpz_t hi)
{
  mpz_t r, t;
  mpz_init(t);
  mpz_init(r);
  mpz_sub(r, hi, lo);
  mpz_isaac_urandomm(t, r);
  mpz_clear(r);
  mpz_add(t, t, lo);
  mpz_sub_ui(t, t, 1);
  _GMP_next_prime(t);
  if (mpz_cmp(t,hi) > 0) {
     mpz_sub_ui(t, lo, 1);
    _GMP_next_prime(t);
    if (mpz_cmp(t,hi) > 0) {
      mpz_clear(t);
      return 0;
    }
  }
  mpz_set(p, t);
  mpz_clear(t);
  return 1;
}

/* TRIVIAL: pick random values until one is prime */
/* Perfect distribution. */
static int _random_prime_trivial(mpz_t p, const mpz_t lo_in, const mpz_t hi_in)
{
  mpz_t r, t, lo, hi;
  int res = 0, tries = 10000;

  if (mpz_cmp_ui(hi_in,2) < 0 || mpz_cmp(lo_in,hi_in) > 0)
    return 0;

  mpz_init_set(lo, lo_in);
  mpz_init_set(hi, hi_in);
  if (mpz_cmp_ui(lo,2) <= 0) {
    mpz_set_ui(lo,1);
  } else if (mpz_even_p(lo)) {
    mpz_add_ui(lo,lo,1);
  }
  if (mpz_cmp_ui(hi,2) <= 0) {
    mpz_set_ui(hi,1);
  } else if (mpz_even_p(hi)) {
    mpz_sub_ui(hi,hi,1);
  }
  /* lo and hi are now odd */
  if (mpz_cmp(lo,hi) >= 0) {
    if (mpz_cmp(lo,hi) > 0) {
      /* null range */
    } else if (mpz_cmp_ui(lo,1) == 0) {
      mpz_set_ui(p,2);
      res = 1;
    } else if (_GMP_is_prob_prime(lo)) {
      mpz_set(p,lo);
      res = 1;
    }
    mpz_clear(hi); mpz_clear(lo);
    return res;
  }
  /* lo and hi are now odd and at least one odd between them */

  mpz_init(t);
  mpz_init(r);
  mpz_sub(r, hi, lo);
  mpz_tdiv_q_2exp(r, r, 1);
  mpz_add_ui(r,r,1);
  do {
    mpz_isaac_urandomm(t, r);
    mpz_mul_2exp(t, t, 1);
    mpz_add(t, t, lo);
    if (mpz_cmp_ui(t,1) == 0) mpz_set_ui(t,2);  /* map 1 back to 2 */
  } while (!_GMP_is_prob_prime(t) && --tries > 0);

  if (tries > 0) {
    mpz_set(p, t);
    res = 1;
  } else {
    /* We couldn't find anything.  Perhaps no primes in range. */
    res = _random_prime_primeinc(p, lo, hi);
  }
  mpz_clear(r); mpz_clear(t); mpz_clear(hi); mpz_clear(lo);
  return res;
}

/* Set p to a random prime between lo and hi inclusive */
int mpz_random_prime(mpz_t p, const mpz_t lo, const mpz_t hi)
{
  return _random_prime_trivial(p,lo,hi);
}

void mpz_random_ndigit_prime(mpz_t p, UV n)
{
  mpz_t lo, hi;
  switch (n) {
    case 0:   mpz_set_ui(p,0); return;
    case 1:   mpz_set_ui(p, pr[isaac_rand(4)]);  return;
    case 2:   mpz_set_ui(p, pr[4+isaac_rand(21)]);  return;
    default:  break;
  }
  mpz_init_set_ui(lo,10);
  mpz_pow_ui(lo, lo, n-1);
  mpz_init(hi);
  mpz_mul_ui(hi, lo, 10);

  if (!mpz_random_prime(p, lo, hi))
    croak("Failed to find %"UVuf" digit prime\n", n);

  mpz_clear(lo);
  mpz_clear(hi);
}

/* Random number rop such that 2*mult*rop+1 has nbits bits. */
static void _rand_in_bit_interval(mpz_t rop, UV nbits, mpz_t mult)
{
  mpz_t t, lo, hi;
  mpz_init(t); mpz_init(lo); mpz_init(hi);

  mpz_mul_ui(t, mult, 2);

  mpz_setbit(lo, nbits-1);
  mpz_sub_ui(lo, lo, 1);
  mpz_cdiv_q(lo, lo, t);   /* lo = ceil(2^(nbits-1)-1 / (2*mult)) */

  mpz_setbit(hi, nbits);
  mpz_sub_ui(hi, hi, 2);
  mpz_fdiv_q(hi, hi, t);   /* hi = floor(2^nbits-2 / (2*mult)) */

  mpz_sub(t, hi, lo);
  mpz_isaac_urandomm(rop, t);
  mpz_add(rop, rop, lo);

  mpz_clear(t); mpz_clear(lo); mpz_clear(hi);
}

#define _SAFE_REJECT(q, p1, p2, p3, p4, p5) \
  { uint32_t qm = mpz_fdiv_ui(q, p1*p2*p3*p4*p5); \
    if ((qm % p1) == 0 || (qm % p1) == (p1>>1)) continue; \
    if ((qm % p2) == 0 || (qm % p2) == (p2>>1)) continue; \
    if ((qm % p3) == 0 || (qm % p3) == (p3>>1)) continue; \
    if ((qm % p4) == 0 || (qm % p4) == (p4>>1)) continue; \
    if ((qm % p5) == 0 || (qm % p5) == (p5>>1)) continue; \
  }

void mpz_random_safe_prime(mpz_t p, UV nbits)
{
  static const unsigned char small_safe[] = {5,7,11,23,47,59,83,107};
  mpz_t q, base;
  uint32_t qmod, PR, tlimit;
  int verbose;
  PRIME_ITERATOR(iter);

  switch (nbits) {
    case 0: case 1: case 2:   mpz_set_ui(p, 0);  return;
    case 3:   mpz_set_ui(p, small_safe[ 0+isaac_rand(2)]);  return;
    case 4:   mpz_set_ui(p, 11);  return;
    case 5:   mpz_set_ui(p, 23);  return;
    case 6:   mpz_set_ui(p, small_safe[ 4+isaac_rand(2)]);  return;
    case 7:   mpz_set_ui(p, small_safe[ 6+isaac_rand(2)]);  return;
    default:  break;
  }

  mpz_init(q); mpz_init(base);

  if (nbits > 35) {
    mpz_isaac_urandomb(base, nbits-35);
    mpz_mul_2exp(base,base,2);
  }
  mpz_setbit(base, nbits-1);
  mpz_setbit(base, 1);
  mpz_setbit(base, 0);
  verbose = get_verbose_level();
  tlimit = (nbits <= 512000) ?  (nbits*(nbits/64.0) + 0.5)  :  4000000000U;
  while (1) {
    /* 1. Generate random nbit p */
    if (nbits > 35) {
      mpz_set_ui(p, isaac_rand32());
      mpz_mul_2exp(p, p, nbits-33);
    } else {
      mpz_isaac_urandomb(p, nbits-3);
      mpz_mul_2exp(p, p, 2);
    }
    mpz_ior(p, p, base);
    /* 2. p = 2q+1  =>  q = p >> 1 */
    mpz_div_2exp(q, p, 1);
    /* 3. Fast compositeness pretests for both q and p at the same time. */
    qmod = mpz_fdiv_ui(q, 1155UL);
    if ( (qmod % 3) != 2 ||
         (qmod % 5) == 0 || (qmod % 7) == 0 || (qmod % 11) == 0 ||
         (qmod % 5) == 2 || (qmod % 7) == 3 || (qmod % 11) == 5) continue;
    if (nbits < 16) {
      /* 4. Pretest that p isn't easily composite */
      if (!primality_pretest(p)) continue;
      /* 5. Pretest that q isn't easily composite */
      if (!primality_pretest(q)) continue;
    } else {
      _SAFE_REJECT(q,  13U,  17U,  19U,  23U,  29U);
      _SAFE_REJECT(q,  31U,  37U,  41U,  43U,  47U);
      _SAFE_REJECT(q,  53U,  59U,  61U,  67U,  71U);
      _SAFE_REJECT(q,  73U,  79U,  83U,  89U,  97U);
      if (tlimit >= 101) {
        prime_iterator_setprime(&iter, PR = 101);
        while (PR <= tlimit) {
          uint32_t qm = mpz_fdiv_ui(q, PR);
          if (qm == 0 || qm == (PR>>1)) break;
          PR = prime_iterator_next(&iter);
        }
        if (PR <= tlimit) continue;
      }
    }
    if (verbose > 2) { printf("."); fflush(stdout); }
    /* 6. BPSW on q and M-R base 2 on p. */
    if (!is_euler_plumb_pseudoprime(q)) continue;  /* Start faster */
    if (verbose > 2) { printf("+"); fflush(stdout); }
    if (!miller_rabin_ui(p, 2)) continue;
    if (verbose > 2) { printf("*"); fflush(stdout); }
    if (!_GMP_is_lucas_pseudoprime(q, 2)) continue;
    if (nbits > 64 && !miller_rabin_ui(q, 2)) continue; /* Verify fast test */
    break;
  }

  mpz_clear(base); mpz_clear(q);
  prime_iterator_destroy(&iter);
}

/* Gordon's algorithm */
void mpz_random_strong_prime(mpz_t p, UV nbits)
{
  mpz_t S, T, R, P0, t, i, j;
  UV rbits, sbits, tbits;

  if (nbits < 128)  croak("random_strong_prime, bits must be >= 128");

  if (nbits < 256) {
    rbits = ((nbits+1) >> 1) - 2;
    sbits = (nbits >> 1) - 20;
    tbits = rbits - 20;
  } else {
    UV N1, N2;
    { /* Calculate FIPS 186-4 C.10 recommended parameter */
      UV t_, l2_;
      for (l2_ = 1, t_ = nbits; t_ >>= 1; ) l2_++;
      N1 =  (nbits/2)-l2_-7;
      N2 = N1/2;
    }
    if (N1 > 200) N1 = 201;
    if (N2 > 100) N2 = 101;
    if (N2 < 100) N2 += N1/4;
    rbits = sbits = N1;
    tbits = N2;
  }

  mpz_init(S);  mpz_init(T);  mpz_init(R);  mpz_init(P0);
  mpz_init(t);  mpz_init(i);  mpz_init(j);

  while (1) {
    mpz_random_nbit_prime(S, sbits);
    mpz_random_nbit_prime(T, tbits);

    _rand_in_bit_interval(i, rbits, T);
    while (1) {
      mpz_mul(t, i, T);
      mpz_mul_ui(t, t, 2);
      mpz_add_ui(R, t, 1);                 /* R = 2*i*T+1 */
      if (_GMP_is_prob_prime(R)) break;
      mpz_add_ui(i,i,1);
    }

    mpz_sub_ui(t, R, 2);
    mpz_powm(P0, S, t, R);
    mpz_mul_ui(P0, P0, 2);
    mpz_mul(P0, P0, S);
    mpz_sub_ui(P0, P0, 1);

    mpz_mul(i, R, S);
    mpz_mul_ui(t, i, 2);
    _rand_in_bit_interval(j, nbits, i);
    while (1) {
      mpz_mul(p, j, t);
      mpz_add(p, p, P0);                 /* p = 2*j*R*S+p0 */
      if (mpz_sizeinbase(p,2) > nbits) break;
      if (_GMP_is_prob_prime(p)) {
        mpz_clear(t);  mpz_clear(i);  mpz_clear(j);
        mpz_clear(S);  mpz_clear(T);  mpz_clear(R);  mpz_clear(P0);
        /* p-1 has factor R.  p+1 has factor S.  r-1 has factor T. */
        return;
      }
      mpz_add_ui(j,j,1);
    }
  }
}

/*===========================================================================*/
/* Proven primes (Maurer and Shawe-Taylor */
/*===========================================================================*/

#define MAKE_PROOF_START(proofptr, n, nums) \
  if (proofptr) { \
    char* thisproof, *thisptr; \
    int prevlen = (*proofptr == 0) ? 0 : strlen(*proofptr); \
    int thislen = (5 + mpz_sizeinbase(n,10)) * nums + 200; \
    New(0, thisproof, thislen + prevlen + 1, char); \
    thisptr = thisproof; \
    thisptr += gmp_sprintf(thisptr,
#define MAKE_PROOF_END(proofptr) \
               ); \
    if (*proofptr) { \
      thisptr += gmp_sprintf(thisptr,"\n"); \
      strcat(thisptr, *proofptr); \
      Safefree(*proofptr); \
    } \
    *proofptr = thisproof; \
  }

#define USE_THEOREM5 0
void mpz_random_maurer_prime(mpz_t n, UV k, char** proofptr)
{
  mpz_t t, a, q, I, R;
  double m, r, minr = USE_THEOREM5 ? 0.334 : 0.5;
  int i, verbose = get_verbose_level();

  /* We could use safely use k <= 64. */
  if (k <= 32)
    return mpz_random_nbit_prime(n, k);

  r = minr;        /* size of q relative to size of n */
  m = 20;          /* always use at least this many bits of randomness */
  if (k > 2*m) {
    do {
      double s = ((double)isaac_rand32()) / ((double)4294967295.0);  /* [0,1] */
      r = pow(2,s-1);  /* exp2 is C99 */
#if USE_THEOREM5
      r = 0.334 + 1.332 * (r-0.5);  /* Stretch r to cover 0.334 - 1 */
#endif
    } while ((k-r*k) <= m);
  }
#if 0  /* Improve efficiency for less than ideal distribution */
  r -= 0.25;  if (r < minr) r = minr;
#endif
  mpz_init(t); mpz_init(a); mpz_init(q); mpz_init(I); mpz_init(R);
  mpz_random_maurer_prime(q, (UV)(r*k)+1, proofptr);
  mpz_setbit(I, k-1);
  mpz_mul_ui(t, q, 2);
  mpz_fdiv_q(I, I, t);   /* I = floor(2^(k-1) / 2q) */
  if (verbose && verbose != 3)
    { gmp_printf("r = %lf  k = %lu  q = %Zd  I = %Zd\n",r,k,q,I); fflush(stdout); }

  while (1) {
    if (verbose > 2) { printf("."); fflush(stdout); }
    mpz_isaac_urandomm(R, I);  /* [0, I-1]  */
    mpz_add(R, R, I);          /* [I, 2I-I] */
    mpz_add_ui(R, R, 1);       /* [I+1,2I]  */
#if USE_THEOREM5
    mpz_setbit(R, 0);          /* We need R to be odd */
#endif
    mpz_mul(n, R, q);
    mpz_mul_ui(n,n,2);
    mpz_add_ui(n,n,1);         /* n = 2Rq+1 */

    if (!primality_pretest(n)) continue;
    if (verbose > 2) { printf("+"); fflush(stdout); }
    /* if (!is_euler_plumb_pseudoprime(n)) continue; */
    if (!miller_rabin_ui(n,2)) continue;
    /* n is a base-2 psp and probably prime */

    if (verbose > 2) { printf("*"); fflush(stdout); }

    /* See if we can use BLS75 theorem 3 */
    mpz_mul_ui(t, q, 2);
    mpz_add_ui(t, t, 1);
    mpz_mul(t, t, t);
    if (mpz_cmp(t, n) > 0) {
      for (i = 0; i < 20; i++) {
        mpz_set_ui(a, pr[i]);
        /* Check A^R mod N != N-1 */
        mpz_powm(a, a, R, n);
        mpz_add_ui(t,a,1);
        if (mpz_cmp(t, n) == 0) continue;
        /* Check A^{Rq} mod N == N-1 */
        mpz_powm(a, a, q, n);
        mpz_add_ui(t,a,1);
        if (mpz_cmp(t, n) != 0) continue;
        if (verbose > 2) { printf("(%"UVuf")",k); fflush(stdout); }
        /* Ensure all results passed BPSW.  ~20% speed penalty. */
        if (!_GMP_is_lucas_pseudoprime(n,2)) croak("Maurer internal failure");
        MAKE_PROOF_START(proofptr, n, 3)
           "Type BLS3\nN %Zd\nQ %Zd\nA %u\n", n, q, pr[i]
        MAKE_PROOF_END(proofptr)
        mpz_clear(t); mpz_clear(a); mpz_clear(q); mpz_clear(I); mpz_clear(R);
        return;
      }
      /* Blast, we couldn't find the right 'a' value fast enough.  Try a new n. */
      continue;
    }

    /* Our q is smaller than sqrt(n)/2-1, so use BLS75 theorem 5. */
#if !USE_THEOREM5
    croak("random_maurer_prime: internal bit size error");
#else
    /* Check for obvious generation problems. */
    if (mpz_even_p(R)) continue;
    if (mpz_cmp_ui(R, 1) <= 0) continue;
    mpz_gcd(t, q, R);  if (mpz_cmp_ui(t, 1) != 0) continue;

    /* Theorem 5 with m = 2, assuming (I) which we'll check after this. */
    {
      mpz_t ts, tr, F;
      mpz_init(ts); mpz_init(tr); mpz_init(F);
      mpz_mul_ui(F, q, 2);
      /* Calculate r,s from page 624 of BLS75 */
      mpz_mul_ui(t, F, 2);
      mpz_tdiv_qr(ts, tr, R, t);
      /* Verify the r,s condition */
      mpz_mul(t, tr, tr);
      mpz_submul_ui(t, ts, 8);   /* t = r^2-8s */
      if (mpz_sgn(ts) != 0 && mpz_perfect_square_p(t)) {
        /* printf("fail r/s check\n"); */
        mpz_clear(ts);  mpz_clear(tr);  mpz_clear(F);
        continue;
      }
      /* Verify size of N with m=2. a,t are temps.  Should not fail. */
      mpz_mul(t, F, tr);
      mpz_add_ui(a, t, 1);         /* a = rF + 1 */
      mpz_sub_ui(tr, tr, 1);
      mpz_mul(t, F, F);
      mpz_mul_ui(t, t, 2);
      mpz_mul(t, t, tr);
      mpz_add(a, a, t);           /* a = (r-1)2FF + rF + 1 */
      mpz_mul(t, F, F);
      mpz_mul(t, t, F);
      mpz_mul_ui(t, t, 4);
      mpz_add(a, a, t);           /* a = 4FFF + (r-1)2FF + rF + 1 */
      mpz_add_ui(t, F, 1);
      mpz_clear(tr);  mpz_clear(ts);  mpz_clear(F);
      if (mpz_cmp(n,a) >= 0) {
        /* printf("fail N size check\n"); */
        continue;
      }
      /* Check divisibility required to use m=2 */
      if (mpz_divisible_p(n,t)) {
        /* printf("fail N divisiblity check\n"); */
        continue;
      }
    }

#define SET_A_CHECK_PSP(i) \
   mpz_set_ui(a, pr[i]); \
   if (apsp[i] == -1) \
     { mpz_sub_ui(t,n,1); mpz_powm(t,a,t,n); apsp[i] = (mpz_cmp_ui(t,1) == 0); } \
   if (apsp[i] == 0) continue;
#define CHECK_GCD(t) \
    mpz_powm(t, a, t, n); mpz_sub_ui(t, t, 1); mpz_gcd(t, t, n); \
    if (mpz_cmp_ui(t,1) != 0) continue;

    {
      int j, apsp[20];  /* apsp caches psp check.  Init all to -1. */
      for (i = 0; i < 20; i++) apsp[i] = -1;
      apsp[0] = 1;  /* We passed a base 2 psp test to get here */
      /* Find an a that works for p=2 */
      for (i = 0; i < 20; i++) {
        SET_A_CHECK_PSP(i);
        mpz_mul(t, q, R);
        CHECK_GCD(t);
        /* We are good for p=2.  Find an a for p=q. */
        for (j = 0; j < 20; j++) {
          SET_A_CHECK_PSP(j);
          mpz_mul_ui(t, R, 2);
          CHECK_GCD(t);
          /* Success */
          if (verbose > 2) { printf("(%lu)",k); fflush(stdout); }
          if (i == 0 && j == 0) {
            MAKE_PROOF_START(proofptr, n, 2)
              "Type BLS5\nN %Zd\nQ[1] %Zd\n----\n", n, q
            MAKE_PROOF_END(proofptr)
          } else {
            MAKE_PROOF_START(proofptr, n, 2)
              "Type BLS5\nN %Zd\nQ[1] %Zd\nA[0]  %lu\nA[1]  %lu\n----\n", n, q, pr[i], pr[j]
            MAKE_PROOF_END(proofptr)
          }
          /* Ensure all results passed BPSW.  ~20% speed penalty. */
          if (!_GMP_is_lucas_pseudoprime(n,2)) croak("Maurer internal failure");
          mpz_clear(t); mpz_clear(a); mpz_clear(q); mpz_clear(I); mpz_clear(R);
          return;
        }
        break; /* Failed for p=q */
      }
    }
    /* Blast, we couldn't find the right 'a' value fast enough.  Try a new n. */
#endif
  }
}

/* FIPS 186-4 algorithm but using our CSPRNG (ISAAC) instead of SHA-256 */
void mpz_random_shawe_taylor_prime(mpz_t c, UV k, char** proofptr)
{
  mpz_t c0, t, u, a, z;

  if (k <= 32)
    return mpz_random_nbit_prime(c, k);

  mpz_init(c0); mpz_init(t); mpz_init(u); mpz_init(a); mpz_init(z);

  mpz_random_shawe_taylor_prime(c0, 1 + (k+1)/2, proofptr);

  mpz_isaac_urandomb(t, k-1);
  mpz_setbit(t,k-1);       /* Steps 18-21: t a random k-bit integer */
  mpz_mul_ui(u, c0, 2);
  mpz_cdiv_q(t, t, u);     /* Step 22: set t based on random integer */

  while (1) {
    /* Steps 23-24 */
    mpz_mul_ui(u, c0, 2);  /* u = 2c0 */
    mpz_mul(c, u, t);
    mpz_add_ui(c, c, 1);   /* c = 2tc0+1 */
    if (mpz_sizeinbase(c,2) > k) {
      mpz_set_ui(t,0);
      mpz_setbit(t,k-1);
      mpz_cdiv_q(t, t, u);
      mpz_mul(c, u, t);
      mpz_add_ui(c, c, 1);
    }
    /* Don't bother with Steps 26-31 for obvious composites */
    if (primality_pretest(c) && miller_rabin_ui(c,2)) {
      /* Steps 26-29 */
      mpz_sub_ui(u, c, 3);
      mpz_isaac_urandomm(a, u);
      mpz_add_ui(a, a, 2);

      /* Step 30 */
      mpz_mul_ui(u, t, 2);
      mpz_powm(z, a, u, c);

      /* Step 31 */
      mpz_sub_ui(u, z, 1);
      mpz_gcd(u, u, c);
      if (mpz_cmp_ui(u, 1) == 0) {
        mpz_powm(u, z, c0, c);
        if (mpz_cmp_ui(u, 1) == 0) {
          /* Ensure all results passed BPSW.  ~20% speed penalty. */
          if (!_GMP_is_lucas_pseudoprime(c,2)) croak("ST internal failure");
          MAKE_PROOF_START(proofptr, c, 3)
             "Type Pocklington\nN %Zd\nQ %Zd\nA %Zd\n", c, c0, a
          MAKE_PROOF_END(proofptr)
          mpz_clear(c0); mpz_clear(t); mpz_clear(u); mpz_clear(a); mpz_clear(z);
          return;
        }
      }
    }
    mpz_add_ui(t,t,1);
  }
}

/*===========================================================================*/