text: /*
   *        Prime number generator and tests utilizing Rabin's
   *        compositeness test and Gordon's strong prime concept.
   *
   *        coded by Jack Lacy, December, 1991
   *
   *        Copyright (c) 1991 Bell Laboratories
   */
  #include "lib9.h"
  #include <libcrypt.h>
  
  static int long_log(int);
  static int fakelog(int);
  static Boolean probPrimeTest(BigInt, BigInt*);
  static Boolean first53Test(BigInt);
  static Boolean genGordonPrimeSet(int, BigInt, int, BigInt);
  static Boolean genNISTPrimeSet(int, BigInt, int, BigInt);
  
  static int first_53_primes[53] = {
          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,
          131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
          181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
          239, 241, 251,
  };
  
  #define RLEN(a, i) ((a >= i) && (a < 2*i))
  
  static int
  long_log(int n)
  {
  	
  	int i, j;
  	
  	i = 8192;
  	j = 13;
  	for(;;) {
  		if (RLEN(n, i))
  			break;
  		else {
  			j++;
  			i *= 2;
  		}
  	}
  	return j;
  }
  
  static int
  fakelog(int n)
  {
  	
  	switch(n) {
  	case 2:    return(1);
  	case 4:    return(2);
  	case 8:    return(3);
  	case 16:   return(4);
  	case 32:   return(5);
  	case 64:   return(6);
  	case 128:  return(7);
  	case 256:  return(8);
  	case 512:  return(9);
  	case 1024: return(10);
  	case 2048: return(11);
  	case 4096: return(12);
  	default:
  		if (RLEN(n, 2))
  			return fakelog(2);
  		else if (RLEN(n, 4))
  			return fakelog(4);
  		else if (RLEN(n, 8))
  			return fakelog(8);
  		else if (RLEN(n, 16))
  			return fakelog(16);
  		else if (RLEN(n, 32))
  			return fakelog(32);
  		else if (RLEN(n, 64))
  			return fakelog(64);
  		else if (RLEN(n, 128))
  			return fakelog(128);
  		else if (RLEN(n, 256))
  			return fakelog(256);
  		else if (RLEN(n, 512))
  			return fakelog(512);
  		else if (RLEN(n, 1024))
  			return fakelog(1024);
  		else if (RLEN(n, 2048))
  			return fakelog(2048);
  		else
  			return long_log(n);
  	}
  }
  
  int primeTestAttempts = 5;
  
  void
  setPrimeTestAttempts(int i)
  {
  	primeTestAttempts = i;
  }
  
  
  static Boolean
  probPrimeTest(BigInt n, BigInt *bigs)
  {
  	Boolean retval = FALSE;
  	int j, k = 0;
  	BigInt nminus1, x, y, q;
  	
  	nminus1 = bigs[0];
  	q       = bigs[1];
  	x       = bigs[2];
  	y       = bigs[3];
  
  	bigsub(n, one, nminus1);
  	bigCopy(nminus1, q);
  	
  	while (even(q)) {
  		k++;
  		bigRightShift(q, (int)1, q);
  	}
  	
  	getRandBetween(one, n, x, PSEUDO);
  	bigPow(x, q, n, y);
  	
  	j = 0;
  	if (bigCompare(y, one) == 0) {
  		retval = TRUE;
  	}
  	
  	while ((j < k) && (retval == FALSE)) {
  		if (bigCompare(y, nminus1) == 0) {
  			retval = TRUE;
  			break;
  		}
  		if (bigCompare(y, one) == 0) {
  			retval = FALSE;
  			break;
  		}
  		j++;
  		bigMultiply(y, y, q);
  		bigMod(q, n, y);
  	}
  	return retval;
  }
  
  
  static  Boolean
  first53Test(BigInt n)
  {
  	BigData np;
  	int i, j;
  	int N;
  	ulong m, divisor;
  	ushort m2;
  	
  	N = (int)LENGTH(n);
  	np = NUM(n);
  	i = 0;
  	do {
  		divisor = first_53_primes[i++];
  		m = 0;
  		for (j=N-1; j >= 0; j--) {
  			m2 = (ushort)(np[j]>>16) & 0xffff;
  			m = ((m<<16) + m2)%divisor;
  			m2 = (ushort)(np[j] & 0xffff);
  			m = ((m<<16) + m2)%divisor;
  		}
  		if (m == 0) {
  			return(FALSE);
  		}
  	} while (i < 53);
  	return(TRUE);
  }
  
  
  static Boolean
  primeTest1(BigInt n, BigInt *bigs)
  {
  	int k;
  	int accuracy;
  	BigData np = NUM(n);
  	
  	if (LENGTH(n) == 1) {
  		k = 0;
  		do {
  			if (np[0] == (NumType)first_53_primes[k++])
  				return(TRUE);
  		} while (k < 53);
  	}
  	
  	if (even(n))
  		return(FALSE);
  	
  	if (first53Test(n) == FALSE)
  		return(FALSE);
  	
  	accuracy = primeTestAttempts;
  	while ((accuracy > 0) && probPrimeTest(n, bigs) == TRUE)
  		accuracy--;
  	
  	if (accuracy == 0)
  		return(TRUE);
  	else
  		return(FALSE);
  }
  
  static BigInt*
  getprimebigs(void)
  {
  	BigInt *bigs;
  
  	bigs = crypt_malloc(4*sizeof(BigInt));
  	bigs[0] = bigInit(0);
  	bigs[1] = bigInit(0);
  	bigs[2] = bigInit(0);
  	bigs[3] = bigInit(0);
  	return bigs;
  }
  
  static void
  freeprimebigs(BigInt *bigs)
  {
  	freeBignum(bigs[0]);
  	freeBignum(bigs[1]);
  	freeBignum(bigs[2]);
  	freeBignum(bigs[3]);
  	free(bigs);
  }
  
  Boolean
  primeTest(BigInt n)
  {
  	BigInt *bigs;
  	Boolean rv;
  
  	bigs = getprimebigs();
  	rv = primeTest1(n, bigs);
  	freeprimebigs(bigs);
  
  	return rv;
  }
  
  void
  getPrime(int numbits, BigInt prime)
  {
  	int numbytes;
  	int shiftdiff;
  	BigInt *bigs;
  	
  	numbytes = (numbits/8);
  	if (numbits%8 != 0)
  		numbytes++;
  	
  	bigRand(numbytes, prime, PSEUDO);
  	
  	shiftdiff = (bigBits(prime) - numbits);
  	if (shiftdiff > 0)
  		bigRightShift(prime, (int)shiftdiff, prime);
  	else if (shiftdiff < 0)
  		bigLeftShift(prime, (int)(-shiftdiff), prime);
  	
  	if (EVEN(prime)) {
  		NUM(prime)[0] |= 1;
  	}
  
  	bigs = getprimebigs();
  	while (1) {
  		if (primeTest1(prime, bigs) == TRUE)
  			break;
  		NUM(prime)[0] += 2;
  	}
  	freeprimebigs(bigs);
  }
  
  /*
     Generate strong primes using Gordon's method and
     Rabin's probabilistic primality test.  This function
     returns the strong prime p as well as r, a prime factor
     of p-1.
     */
  static Boolean
  genGordonPrimeSet(int numbits, BigInt prime, int factorbits, BigInt factor)
  {
  	BigInt r, pzero, twors;
  	BigInt s, t, p, rs, ss, rr;
  	BigInt rminus1, sminus1, twot;
  	int pbits, tbits, rbits, sbits;
  	Boolean found;
  	BigInt *bigs;
  	
  	if (numbits < 32)
  		return genNISTPrimeSet(numbits, prime, factorbits, factor);
  	
  	/* the two return values */
  	p = prime;
  	r = factor;
  	
  	pbits = numbits;
  	rbits = (pbits - fakelog(pbits) - 1);
  	sbits = rbits/2;
  	rbits = sbits;
  	tbits = rbits - fakelog(rbits) - 1;
  	
  	s = bigInit(0);
  	t = bigInit(0);
  	getPrime(sbits, s);
  	getPrime(tbits, t);
  	
  	/* find r -- r = 2Lt + 1 and is prime */
  	twot = bigInit(0);
  	bigLeftShift(t, 1, twot);
  	
  	reset_big(r, 1);
  	bigAdd(r, twot, r);
  
  	bigs = getprimebigs();
  	while (primeTest1(r, bigs) == FALSE)
  		bigAdd(r, twot, r);
  	
  	freeBignum(t);
  	freeBignum(twot);
  	
  	/*
  	   find p -- p = p0 + 2krs where:
  	   p0 = u(r,s);        u(r,s) odd
  	   p0 = u(r,s) + rs;   u(r,s) even
  	   u(r,s) = (s^(r-1) - r^(s-1))mod rs.
  	   */
  	rs = bigInit(0);
  	rminus1 = bigInit(0);
  	sminus1 = bigInit(0);
  	
  	bigMultiply(r, s, rs);
  	bigSubtract(r, one, rminus1);
  	bigSubtract(s, one, sminus1);
  	
  	ss = bigInit(0);
  	rr = bigInit(0);
  	bigPow(s, rminus1, rs, ss);
  	bigPow(r, sminus1, rs, rr);
  	
  	pzero = bigInit(0);
  	bigSubtract(ss, rr, pzero);
  	
  	if (SIGN(pzero) == NEG) {
  		negate(pzero, rs, pzero);
  	}
  	
  	if (EVEN(pzero))
  		bigAdd(pzero, rs, pzero);
  	
  	freeBignum(s);
  	freeBignum(ss);
  	freeBignum(rr);
  	freeBignum(rminus1);
  	freeBignum(sminus1);
  	
  	twors = bigInit(0);
  	bigLeftShift(rs, 1, twors);
  	freeBignum(rs);
  	
  	reset_big(p, (NumType)0);
  	bigAdd(pzero, p, p);
  	
  	while (bigBits(p) < pbits)
  		bigAdd(p, twors, p);
  	
  	found = TRUE;
  	
  	while (primeTest1(p, bigs) == FALSE) {
  		bigAdd(p, twors, p);
  		if (bigBits(p) > pbits) {
  			found = FALSE;
  			break;
  		}
  	}
  
  	freeprimebigs(bigs);
  	
  	freeBignum(twors);
  	
  #ifndef NDEBUG
  	if (found) {
  		BigInt pminus1;
  		/* verify that p and r are ok */
  		pminus1 = bigInit(0);
  		bigSubtract(p, one, pminus1);
  		bigMod(pminus1, r, pzero);
  		freeBignum(pminus1);
  		
  		if (bigCompare(pzero, zero) != 0)
  			handle_exception(WARNING, "genGordonPrimeSet: DEBUG: pzero non-zero.");
  	}
  #endif
  	
  	freeBignum(pzero);
  	
  	return found;
  }
  
  static Boolean
  genNISTPrimeSet(int numbits, BigInt prime, int facbits, BigInt factor)
  {
  	int qlen;
  	BigInt p, q, twoq, tmp, ignore, n;
  	BigInt smallseed, bigseed;
  	Boolean primeVal;
  	BigInt *bigs;
  
  	p                = prime;
  	q                = factor;
  	twoq        = bigInit(0);
  	tmp                = bigInit(0);
  	ignore        = bigInit(0);
  	n                = bigInit(0);
  	smallseed        = bigInit(0);
  	bigseed        = bigInit(0);
  	
  	qlen = facbits;
  		/* was commented out */
  		if (numbits >= 512)
  			qlen = 160;
  		else
  			qlen = (numbits/2) - ((3*numbits)/16);
  	getPrime(qlen, q);
  	bigLeftShift(q, (int)1, twoq);
  	
  	bigLeftShift(one, numbits, tmp);
  	bigSubtract(tmp, one, tmp);
  	bigDivide(tmp, twoq, bigseed, ignore);
  	
  	bigLeftShift(one, (int)(numbits-1), tmp);
  	bigSubtract(tmp, one, tmp);
  	bigDivide(tmp, twoq, smallseed, ignore);
  	
  	freeBignum(tmp);
  	freeBignum(ignore);
  	
  	primeVal = FALSE;
  	getRandBetween(bigseed, smallseed, n, PSEUDO);
  	bigSubtract(n, one, n);
  
  	bigs = getprimebigs();
  	while (primeVal == FALSE) {
  		bigAdd(n, one, n);
  		bigMultiply(n, twoq, p);
  		bigAdd(p, one, p);
  		primeVal = primeTest1(p, bigs);
  	}
  	freeprimebigs(bigs);
  	freeBignum(smallseed);
  	freeBignum(bigseed);
  	freeBignum(twoq);
  	freeBignum(n);
  	
  	return primeVal;
  	
  }
  
  
  void
  genStrongPrimeSet(int numbits, BigInt prime, int facbits, BigInt factor, PrimeType type)
  {
  	if (type == NIST)
  		genNISTPrimeSet(numbits, prime, facbits, factor);
  	else if (type == GORDON) {
  		while (genGordonPrimeSet(numbits, prime, facbits, factor) == FALSE) {
  			reset_big(prime, (NumType)0);
  			reset_big(factor, (NumType)0);
  		}
  	}
  }
  
  void
  genStrongPrime(int numbits, BigInt prime)
  {
  	BigInt factor;
  	
  	factor = bigInit(0);
  	
  	genStrongPrimeSet(numbits, prime, 160, factor, NIST);
  	
  	freeBignum(factor);
  }
  
  /* f is a factor of p-1 */
  void
  getPrimitiveElement(BigInt a, BigInt p, BigInt f)
  {
  	BigInt pminus1, tmp, d;
  	
  	pminus1 = bigInit(0);
  	tmp = bigInit(0);
  	d = bigInit(0);
  	
  	bigSubtract(p, one, pminus1);
  	bigDivide(pminus1, f, d, tmp);
  	getRandBetween(p, one, a, PSEUDO);
  	bigPow(a, d, p, tmp);
  	while (bigCompare(tmp, one) == 0) {
  		bigAdd(a, one, a);
  		bigPow(a, d, p, tmp);
  	}
  	
  	freeBignum(pminus1);
  	freeBignum(tmp);
  	freeBignum(d);
  }
  
  