/* Math.java             nov 2000
 *
 * Erwan Lemonnier & Eric Nordenstam
 *
 * A personal implementation for a set of operation on large 
 * numbers.
 ***********************************************************************/

import java.math.BigInteger;
import java.util.Random;

/**
 * An implementation for operations on large integers.
 * <p>
 * @author E.Lemonnier & E.Nordenstam
 */
 
public class Math {
	
	private static BigInteger ZERO = BigInteger.ZERO;
	private static BigInteger ONE  = BigInteger.ONE;
	private static BigInteger TWO  = new BigInteger("2");
	
	//list of the prime numbers <100. Used for primality testing (see polFactor())
	private static BigInteger[] primes =
		{	new BigInteger("3"),  new BigInteger("53"), 
			new BigInteger("5"),  new BigInteger("59"),
			new BigInteger("7"),  new BigInteger("47"),
			new BigInteger("11"), new BigInteger("61"), 
			new BigInteger("13"), new BigInteger("67"), 
			new BigInteger("17"), new BigInteger("71"), 
			new BigInteger("19"), new BigInteger("73"), 
			new BigInteger("23"), new BigInteger("79"), 
			new BigInteger("27"), new BigInteger("83"),  
			new BigInteger("29"), new BigInteger("89"), 
			new BigInteger("31"), new BigInteger("97"),  
			new BigInteger("37"), new BigInteger("41"), 
			new BigInteger("43")			
		};
	
	
	/**
	 * Compute the greatest common divisor (gcd) of two integers. It uses
	 * the Euclidean algorithm.
	 * <p>
	 * @param a large integer.
	 * @param b large integer.
	 * @return the gcd of a and b.
	 */
	public static BigInteger gcd(BigInteger a, BigInteger b) {
		BigInteger t;
		
		while (b.compareTo(ZERO) != 0) {
			t = a.mod(b);
			a = b;
			b = t;
		}	
		
		return a;
	}
	
	
	/**
	 * Compute the inverse of an integer <code>a</code> modulo <code>n</code>.
	 * It is assumed that <code>a</code> and <code>n</code> are relative primes,
	 * ie <code>gcd(a,n)=1</code>, otherwise an <code>ArithmeticExpression</code>
	 * is thrown. Uses the extended Euclidean algorithm.
	 * <p>
	 * @param a integer to inverse.
	 * @param n modulo.
	 * @return the inverse of a modulo n.
	 */
	//This implementation follows the pseudocode given in
	//'Cryptographie, Theorie and Practice' (Douglas R. Stinson)
	public static BigInteger modInverse(BigInteger a, BigInteger n) throws ArithmeticException {
		BigInteger n0 = n;
		BigInteger a0 = a;
		BigInteger t0 = ZERO;
		BigInteger t  = ONE;
		BigInteger q  = n0.divide(a0);
		BigInteger r  = n0.subtract(q.multiply(a0));
		BigInteger temp;
		
		while (r.compareTo(ZERO) == 1) {
			temp = (t0.subtract(q.multiply(t))).mod(n);
							
			t0 = t;
			t  = temp;
			n0 = a0;
			a0 = r;
			q  = n0.divide(a0);
			r  = n0.subtract(q.multiply(a0));
		}
		
		if (a0.compareTo(ONE) != 0)	throw new ArithmeticException();
		
		return t.mod(n);
	}
	
	/**
	 * Compute <code>a</code> power <code>e</code> modulo <code>n</code>.
	 * <p>
	 * @param a integer to exponentiate.
	 * @param e exponent.
	 * @param n modulo.
	 * @return a^e mod p.
	 */
	public static BigInteger modPow(BigInteger a, BigInteger e, BigInteger n) {
		BigInteger r = ONE;
		
		//ameliorer avec lookup
		for(int i=e.bitLength()-1; i>=0; i--) {
			r = (r.multiply(r)).mod(n);
			if (e.testBit(i))
				r = (r.multiply(a)).mod(n);
		}
		
		return r;
	}

	/**
	 * Same as <code>modPow</code>. In theory, it should be faster than modPow
	 * for large numbers, but somehow, java makes it slower...
	 * <p>
	 * @param a integer to exponentiate.
	 * @param e exponent.
	 * @param n modulo.
	 * @return a^e mod p.
	 */
	 //like modPow, but check the bits of the exponents by group of 4 at a time.
	 //Theoretically faster for 512 bits integers.
	public static BigInteger mod4Pow(BigInteger a, BigInteger e, BigInteger n) {
		BigInteger[] val = new BigInteger[16];
		val[0] = ONE;
		for(int i=1; i<16; i++) {
			val[i] = val[i-1].multiply(a);
		}
		
		BigInteger r = ONE;
		byte[] b = e.toByteArray();
		
		byte LOW = (byte)15;
		byte HIGH = (byte)240;
		byte bits;
		int p;
		for(int i=0; i<b.length; i++) {

			//catch 4 bits of higher value			
			p = (int)(b[i] & HIGH);
			if (p<0)	p+=256;
			p=p/16;
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(val[p])).mod(n);
			
			//catch 4 lower bits
			p = (int)(b[i] & LOW);
			if (p<0)	p+=256;
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(r)).mod(n);
			r = (r.multiply(val[p])).mod(n);
		}
		return r;
	}
	
	
	/**
	 * Factorises a large integer. Uses Pollard's Rho algorithm.
	 * <p>
	 * @param n an integer to factorise.
	 * @return a factor of this integer. It can be a composite number
	 * or the integer itself. 
	 */
	public static BigInteger polFactors(BigInteger n) {
		BigInteger xi = new BigInteger(100, new Random());
		BigInteger x2i = ((xi.multiply(xi)).add(ONE)).mod(n);
		BigInteger d;
		
		while(true) {
			d = gcd(x2i.subtract(xi), n);
			
			if (d.compareTo(ONE)>0)	return d;
			
			xi = ((xi.multiply(xi)).add(ONE)).mod(n);
			x2i = ((x2i.multiply(x2i)).add(ONE)).mod(n);
			x2i = ((x2i.multiply(x2i)).add(ONE)).mod(n);
		}
	}
	
	/**
	 * Create a large prime integer. This integer has a length of
	 * <code>bitLength</code> bits, and its first and last bits are set
	 * to 1. <code>isProbablePrime</code> is used to check whether the generated
	 * integer is prime. The generated integer has a probability
	 * of being prime equal to: 1-1/2^<code>certainty</code>.
	 * <p>
	 * @param bitLength the number of bits in the integer.
	 * @param certainty degree of probability of the integer of being prime.
	 * @return a prime integer.
	 */
	public static BigInteger makePrime(int bitLength, int certainty) {
		BigInteger a = new BigInteger(bitLength, new Random());
		a = a.setBit(bitLength-1);
		a = a.setBit(0);
		
		while (!isProbablePrime(a,certainty)) {
			a = a.add(TWO);
		}
		return a;
	}
	
	/**
	 * Check if an integer is prime. It first test whether the integer
	 * is 1 or even. Then it uses an array of all the primes between 3 and
	 * 100, and check whether the integer is equal to or multiple of one
	 * of them. Finally, it uses the Miller-Rabin probabilistic algorithm to
	 * check if the integer is prime. The probability of the integer of being
	 * prime is then 1-1/2^<code>certainty</code>.
	 * <p>
	 * @param n the integer to check.
	 * @param certainty degree of probability of the integer of being prime.
	 * @return true if it is prime, false otherwise.
	 */
	public static boolean isProbablePrime(BigInteger n, int certainty) {
		//test if n=1
		if (n.compareTo(TWO)<1)
			return true;
		
		//test if n is even
		if (!n.testBit(0))
			return false;
		
		//try to divide by the first primes between 3 and 100
		for (int i=0; i<primes.length; i++) {
			if (n.compareTo(primes[i])==0)
				return true;
			if ( (n.mod(primes[i])).compareTo(ZERO)==0 )
				return false;
		}		
		
		// the Miller-Rabin algorithm itself
		BigInteger n1 = n.subtract(ONE);
		BigInteger t = n1;
		
		int s = 0;
		while(!t.testBit(0)) {			
			t = t.shiftRight(1);
			s++;
		}		
		
		BigInteger a = new BigInteger("3");		//start the sequence at 3
		
		for(int j=0; j<certainty; j++) {
			BigInteger u0 = modPow(a,t,n);
			if (!testPrime(u0, a, n1, s, n))
				return false;
			a = a.add(ONE);
		}
			
		return true;
	}
	
	//part of the Miller-Rabin algorithm. See 'isProbablePrime()'
	private static boolean testPrime(BigInteger u, BigInteger a, BigInteger n1, int s, BigInteger n) {
		if ( u.compareTo(ONE)==0)
			return true;
		
		for (int i=0; i<s; i++) {
			
			if ( u.compareTo(n1)==0 )
				return true;
			else
				u = (u.multiply(u)).mod(n);
		}
		
		return false;
	}
	
	
}