/*
 * Based on the public domain C reference code for New Hope.
 * This Java version is also placed into the public domain.
 * 
 * Original authors: Erdem Alkim, Léo Ducas, Thomas Pöppelmann, Peter Schwabe
 * Java port: Rhys Weatherley
 */

package com.southernstorm.noise.crypto;

import java.security.SecureRandom;
import java.util.Arrays;

/**
 * NewHope key exchange algorithm.
 * 
 * This class implements the standard "ref" version of the New Hope
 * algorithm.
 * 
 * @see NewHopeTor
 */
public class NewHope {

	// -------------- params.h --------------
	
	static final int PARAM_N = 1024;
	static final int PARAM_Q = 12289;
	static final int POLY_BYTES = 1792;
	static final int SEEDBYTES = 32;
	static final int RECBYTES = 256;

	/**
	 * Number of bytes in the public key value sent by Alice.
	 */
	public static final int SENDABYTES = POLY_BYTES + SEEDBYTES;
	
	/**
	 * Number of bytes in the public key value sent by Bob.
	 */
	public static final int SENDBBYTES = POLY_BYTES + RECBYTES;
	
	/**
	 * Number of bytes in shared secret values computed by shareda() and sharedb().
	 */
	public static final int SHAREDBYTES = 32;

	// -------------- newhope.c --------------

	private Poly sk;

	/**
	 * Constructs a NewHope object.
	 */
	public NewHope()
	{
		sk = null;
	}

	@Override
	protected void finalize()
	{
		destroy();
	}

	/**
	 * Destroys sensitive material in this object.
	 * 
	 * This function should be called once the application has finished
	 * with the private key contained in this object.  This function
	 * will also be called when the object is finalized, but the point
	 * of finalization is unpredictable.  This function provides a more
	 * predictable place where the sensitive data is destroyed.
	 */
	public void destroy()
	{
		if (sk != null) {
			sk.destroy();
			sk = null;
		}
	}

	/**
	 * Generates the keypair for Alice.
	 * 
	 * @param send Buffer to place the public key for Alice in, to be sent to Bob.
	 * @param sendOffset Offset of the first byte in the send buffer to populate.
	 * 
	 * The send buffer must have space for at least NewHope.SENDABYTES bytes
	 * starting at sendOffset.
	 * 
	 * @see #sharedb(byte[], int, byte[], int, byte[], int)
   * @see #shareda(byte[], int, byte[], int)
	 */
	public void keygen(byte[] send, int sendOffset)
	{
	  Poly a = new Poly();
	  Poly e = new Poly();
	  Poly r = new Poly();
	  Poly pk = new Poly();
	  byte[] seed = new byte [SEEDBYTES + 32];
	  byte[] noiseseed = new byte [32];

	  try {
		  randombytes(seed);
		  sha3256(seed, 0, seed, 0, SEEDBYTES); /* Don't send output of system RNG */
		  System.arraycopy(seed, SEEDBYTES, noiseseed, 0, 32);
	
		  uniform(a.coeffs, seed);
	
		  if (sk == null)
			  sk = new Poly();
		  sk.getnoise(noiseseed,(byte)0);
		  sk.ntt();
	
		  e.getnoise(noiseseed,(byte)1);
		  e.ntt();
	
		  r.pointwise(sk,a);
		  pk.add(e,r);
	
		  encode_a(send, sendOffset, pk, seed);
	  } finally {
		  a.destroy();
		  e.destroy();
		  r.destroy();
		  pk.destroy();
		  Arrays.fill(seed, (byte)0);
		  Arrays.fill(noiseseed, (byte)0);
	  }
	}

	/**
	 * Generates the public key and shared secret for Bob.
	 * 
	 * @param sharedkey Buffer to place the shared secret for Bob in.
	 * @param sharedkeyOffset Offset of the first byte in the sharedkey buffer to populate.
	 * @param send Buffer to place the public key for Bob in to be sent to Alice.
	 * @param sendOffset Offset of the first byte in the send buffer to populate.
	 * @param received Buffer containing the public key value received from Alice.
	 * @param receivedOffset Offset of the first byte of the value received from Alice.
	 * 
	 * The sharedkey buffer must have space for at least NewHope.SHAREDBYTES
	 * bytes starting at sharedkeyOffset.
	 * 
	 * The send buffer must have space for at least NewHope.SENDBBYTES bytes
	 * starting at sendOffset.
	 * 
	 * The received buffer must have space for at least NewHope.SENDABYTES
	 * bytes starting at receivedOffset.
	 * 
	 * @see #shareda(byte[], int, byte[], int)
   * @see #keygen(byte[], int)
	 */
	public void sharedb(byte[] sharedkey, int sharedkeyOffset,
						byte[] send, int sendOffset,
						byte[] received, int receivedOffset)
	{
	  Poly sp = new Poly();
	  Poly ep = new Poly();
	  Poly v = new Poly();
	  Poly a = new Poly();
	  Poly pka = new Poly();
	  Poly c = new Poly();
	  Poly epp = new Poly();
	  Poly bp = new Poly();
	  byte[] seed = new byte [SEEDBYTES];
	  byte[] noiseseed = new byte [32];
	  byte[] skey = new byte [32];

	  try {
		  randombytes(noiseseed);
	
		  decode_a(pka, seed, received, receivedOffset);
		  uniform(a.coeffs, seed);
	
		  sp.getnoise(noiseseed,(byte)0);
		  sp.ntt();
		  ep.getnoise(noiseseed,(byte)1);
		  ep.ntt();
	
		  bp.pointwise(a, sp);
		  bp.add(bp, ep);
	
		  v.pointwise(pka, sp);
		  v.invntt();
	
		  epp.getnoise(noiseseed,(byte)2);
		  v.add(v, epp);
	
		  helprec(c, v, noiseseed, (byte)3);
	
		  encode_b(send, sendOffset, bp, c);
	
		  rec(skey, v, c);
	
		  sha3256(sharedkey, sharedkeyOffset, skey, 0, 32);
	  } finally {
		  sp.destroy();
		  ep.destroy();
		  v.destroy();
		  a.destroy();
		  pka.destroy();
		  c.destroy();
		  epp.destroy();
		  bp.destroy();
		  Arrays.fill(seed, (byte)0);
		  Arrays.fill(noiseseed, (byte)0);
		  Arrays.fill(skey, (byte)0);
	  }
	}

	/**
	 * Generates the shared secret for Alice.
	 * 
	 * @param sharedkey Buffer to place the shared secret for Alice in.
	 * @param sharedkeyOffset Offset of the first byte in the sharedkey buffer to populate.
	 * @param received Buffer containing the public key value received from Bob.
	 * @param receivedOffset Offset of the first byte of the value received from Bob.
	 * 
	 * The sharedkey buffer must have space for at least NewHope.SHAREDBYTES
	 * bytes starting at sharedkeyOffset.
	 * 
	 * The received buffer must have space for at least NewHope.SENDBBYTES bytes
	 * starting at receivedOffset.
	 * 
	 * @see #shareda(byte[], int, byte[], int)
   * @see #keygen(byte[], int)
	 */
	public void shareda(byte[] sharedkey, int sharedkeyOffset,
						byte[] received, int receivedOffset)
	{
	  Poly v = new Poly();
	  Poly bp = new Poly();
	  Poly c = new Poly();
	  byte[] skey = new byte [32];

	  try {
		  decode_b(bp, c, received, receivedOffset);
	
		  v.pointwise(sk,bp);
		  v.invntt();
	
		  rec(skey, v, c);
		  
		  sha3256(sharedkey, sharedkeyOffset, skey, 0, 32);
	  } finally {
		  v.destroy();
		  bp.destroy();
		  c.destroy();
		  Arrays.fill(skey, (byte)0);
	  }
	}

	/**
	 * Generates random bytes for use in the NewHope implementation.
	 * 
	 * @param buffer The buffer to fill with random bytes.
	 * 
	 * This function may be overridden in subclasses to provide a better
	 * random number generator or to provide static data for test vectors.
	 */
	protected void randombytes(byte[] buffer)
	{
		SecureRandom random = new SecureRandom();
		random.nextBytes(buffer);
	}

	private static void encode_a(byte[] r, int roffset, Poly pk, byte[] seed)
	{
	  int i;
	  pk.tobytes(r, roffset);
	  for(i=0;i<SEEDBYTES;i++)
	    r[POLY_BYTES+roffset+i] = seed[i];
	}

	private static void decode_a(Poly pk, byte[] seed, byte[] r, int roffset)
	{
	  int i;
	  pk.frombytes(r, roffset);
	  for(i=0;i<SEEDBYTES;i++)
	    seed[i] = r[POLY_BYTES+roffset+i];
	}

	private static void encode_b(byte[] r, int roffset, Poly b, Poly c)
	{
	  int i;
	  b.tobytes(r,roffset);
	  for(i=0;i<PARAM_N/4;i++)
	    r[POLY_BYTES+roffset+i] = (byte)(c.coeffs[4*i] | (c.coeffs[4*i+1] << 2) | (c.coeffs[4*i+2] << 4) | (c.coeffs[4*i+3] << 6));
	}

	private static void decode_b(Poly b, Poly c, byte[] r, int roffset)
	{
	  int i;
	  b.frombytes(r, roffset);
	  for(i=0;i<PARAM_N/4;i++)
	  {
	    c.coeffs[4*i+0] = (char)( r[POLY_BYTES+roffset+i]       & 0x03);
	    c.coeffs[4*i+1] = (char)((r[POLY_BYTES+roffset+i] >> 2) & 0x03);
	    c.coeffs[4*i+2] = (char)((r[POLY_BYTES+roffset+i] >> 4) & 0x03);
	    c.coeffs[4*i+3] = (char)(((r[POLY_BYTES+roffset+i] & 0xff) >> 6));
	  }
	}

	// -------------- poly.c --------------
	
	private class Poly
	{
		public char[] coeffs;

		public Poly()
		{
			coeffs = new char [PARAM_N];
		}
		
		protected void finalize()
		{
			destroy();
		}
		
		public void destroy()
		{
			Arrays.fill(coeffs, (char)0);
		}
		
		public void frombytes(byte[] a, int offset)
		{
			int i;
			for (i = 0; i < PARAM_N/4; i++)
			{
			    coeffs[4*i+0] = (char)(                                   (a[offset+7*i+0] & 0xff)       | ((a[offset+7*i+1] & 0x3f) << 8));
			    coeffs[4*i+1] = (char)(((a[offset+7*i+1] & 0xc0) >> 6) | ((a[offset+7*i+2] & 0xff) << 2) | ((a[offset+7*i+3] & 0x0f) << 10));
			    coeffs[4*i+2] = (char)(((a[offset+7*i+3] & 0xf0) >> 4) | ((a[offset+7*i+4] & 0xff) << 4) | ((a[offset+7*i+5] & 0x03) << 12));
			    coeffs[4*i+3] = (char)(((a[offset+7*i+5] & 0xfc) >> 2) | ((a[offset+7*i+6] & 0xff) << 6));
			}
		}

		public void tobytes(byte[] r, int offset)
		{
			int i;
			int t0,t1,t2,t3,m;
			int c;
			for (i = 0; i < PARAM_N/4; i++)
			{
			    t0 = barrett_reduce(coeffs[4*i+0]); //Make sure that coefficients have only 14 bits
			    t1 = barrett_reduce(coeffs[4*i+1]);
			    t2 = barrett_reduce(coeffs[4*i+2]);
			    t3 = barrett_reduce(coeffs[4*i+3]);
			
			    m = t0 - PARAM_Q;
			    c = m;
			    c >>= 15;
			    t0 = m ^ ((t0^m)&c); // <Make sure that coefficients are in [0,q]
			
			    m = t1 - PARAM_Q;
			    c = m;
			    c >>= 15;
			    t1 = m ^ ((t1^m)&c); // <Make sure that coefficients are in [0,q]
			
			    m = t2 - PARAM_Q;
			    c = m;
			    c >>= 15;
			    t2 = m ^ ((t2^m)&c); // <Make sure that coefficients are in [0,q]
			
			    m = t3 - PARAM_Q;
			    c = m;
			    c >>= 15;
			    t3 = m ^ ((t3^m)&c); // <Make sure that coefficients are in [0,q]
			
			    r[offset+7*i+0] = (byte)(t0 & 0xff);
			    r[offset+7*i+1] = (byte)((t0 >> 8) | (t1 << 6));
			    r[offset+7*i+2] = (byte)(t1 >> 2);
			    r[offset+7*i+3] = (byte)((t1 >> 10) | (t2 << 4));
			    r[offset+7*i+4] = (byte)(t2 >> 4);
			    r[offset+7*i+5] = (byte)((t2 >> 12) | (t3 << 2));
			    r[offset+7*i+6] = (byte)(t3 >> 6);
			}
		}

		public void getnoise(byte[] seed, byte nonce)
		{
		  byte[] buf = new byte [4*PARAM_N];
		  int /*t, d,*/ a, b;
		  int i/*,j*/;

		  try {
			  crypto_stream_chacha20(buf,0,4*PARAM_N,nonce,seed);
	
			  for(i=0;i<PARAM_N;i++)
			  {
				/*
				The original C reference code:
				
			    t = (buf[4*i] & 0xff) | (((buf[4*i+1]) & 0xff) << 8) | (((buf[4*i+2]) & 0xff) << 16) | (((buf[4*i+3]) & 0xff) << 24);
			    d = 0;
			    for(j=0;j<8;j++)
			      d += (t >>> j) & 0x01010101;
			    a = ((d >>> 8) & 0xff) + (d & 0xff);
			    b = (d >>> 24) + ((d >>> 16) & 0xff);
			    
			    What the above is doing is reading 32-bit words from buf and then
			    setting a and b to the number of 1 bits in the low and high 16 bits.
			    We instead use the following technique from "Bit Twiddling Hacks",
			    modified for 16-bit quantities:
			    
			    https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel
			    */
				a = (buf[4*i] & 0xff) | (((buf[4*i+1]) & 0xff) << 8);
				a = a - ((a >> 1) & 0x5555);
				a = (a & 0x3333) + ((a >> 2) & 0x3333);
				a = ((a >> 4) + a) & 0x0F0F;
				a = ((a >> 8) + a) & 0x00FF;

				b = (buf[4*i+2] & 0xff) | (((buf[4*i+3]) & 0xff) << 8);
				b = b - ((b >> 1) & 0x5555);
				b = (b & 0x3333) + ((b >> 2) & 0x3333);
				b = ((b >> 4) + b) & 0x0F0F;
				b = ((b >> 8) + b) & 0x00FF;

				coeffs[i] = (char)(a + PARAM_Q - b);
			  }
		  } finally {
			  Arrays.fill(buf, (byte)0);
		  }
		}
		
		public void pointwise(Poly a, Poly b)
		{
		  int i;
		  int t;
		  for(i=0;i<PARAM_N;i++)
		  {
		    t       = montgomery_reduce(3186*b.coeffs[i]); /* t is now in Montgomery domain */
		    coeffs[i] = (char)montgomery_reduce(a.coeffs[i] * t); /* coeffs[i] is back in normal domain */
		  }
		}

		public void add(Poly a, Poly b)
		{
		  int i;
		  for(i=0;i<PARAM_N;i++)
		    coeffs[i] = (char)barrett_reduce(a.coeffs[i] + b.coeffs[i]);
		}

		public void ntt()
		{
		  mul_coefficients(coeffs, psis_bitrev_montgomery);
		  ntt_global(coeffs, omegas_montgomery);
		}

		public void invntt()
		{
		  bitrev_vector(coeffs);
		  ntt_global(coeffs, omegas_inv_montgomery);
		  mul_coefficients(coeffs, psis_inv_montgomery);
		}
	}

	/**
	 * Derives the public "a" value from a 32-byte seed.
	 *  
	 * @param coeffs The 1024 16-bit coefficients of "a" on exit.
	 * @param seed The 32-byte seed to use to generate "a".
	 * 
	 * The base class implementation is not constant-time but usually
	 * this doesn't matter for the public "a" value.  However, as
	 * described in the New Hope paper, non constant-time generation
	 * of "a" can be a problem in anonymity networks like Tor.
	 * 
	 * This function can be overridden in subclasses to provide a
	 * different method for generating "a".  The NewHopeTor class
	 * provides such an example.
	 * 
	 * Reference: https://cryptojedi.org/papers/newhope-20160803.pdf
	 */
	protected void uniform(char[] coeffs, byte[] seed)
	{
	  int pos=0, ctr=0;
	  int val;
	  long[] state = new long [25];
	  int nblocks=14;
	  byte[] buf = new byte [SHAKE128_RATE*nblocks];

	  try {
		  shake128_absorb(state, seed, 0, SEEDBYTES);

		  shake128_squeezeblocks(buf, 0, nblocks, state);

		  while(ctr < PARAM_N)
		  {
		    val = ((buf[pos] & 0xff) | ((buf[pos+1] & 0xff) << 8));
		    if(val < 5*PARAM_Q)
		      coeffs[ctr++] = (char)val;
		    pos += 2;
		    if(pos > SHAKE128_RATE*nblocks-2)
		    {
		      nblocks=1;
		      shake128_squeezeblocks(buf,0,nblocks,state);
		      pos = 0;
		    }
		  }
	  } finally {
		  Arrays.fill(state, 0);
		  Arrays.fill(buf, (byte)0);
	  }
	}

	// -------------- reduce.c --------------

	private static final int qinv = 12287; // -inverse_mod(p,2^18)
	private static final int rlog = 18;

	private static int montgomery_reduce(int a)
	{
	  int u;

	  u = (a * qinv);
	  u &= ((1<<rlog)-1);
	  u *= PARAM_Q;
	  a = a + u;
	  return a >>> 18;
	}

	private static int barrett_reduce(int a)
	{
	  int u;
	  a &= 0xffff;
	  u = (a * 5) >> 16;
	  u *= PARAM_Q;
	  a -= u;
	  return a & 0xffff;
	}
	
	// -------------- error_correction.c --------------

	private static int abs(int v)
	{
	  int mask = v >> 31;
	  return (v ^ mask) - mask;
	}

	private static int f(int[] v0, int v0offset, int[] v1, int v1offset, int x)
	{
	  int xit, t, r, b;

	  // Next 6 lines compute t = x/PARAM_Q;
	  b = x*2730;
	  t = b >> 25;
	  b = x - t*12289;
	  b = 12288 - b;
	  b >>= 31;
	  t -= b;

	  r = t & 1;
	  xit = (t>>1);
	  v0[v0offset] = xit+r; // v0 = round(x/(2*PARAM_Q))

	  t -= 1;
	  r = t & 1;
	  v1[v1offset] = (t>>1)+r;

	  return abs(x-((v0[v0offset])*2*PARAM_Q));
	}

	private static int g(int x)
	{
	  int t,c,b;

	  // Next 6 lines compute t = x/(4*PARAM_Q);
	  b = x*2730;
	  t = b >> 27;
	  b = x - t*49156;
	  b = 49155 - b;
	  b >>= 31;
	  t -= b;

	  c = t & 1;
	  t = (t >> 1) + c; // t = round(x/(8*PARAM_Q))

	  t *= 8*PARAM_Q;

	  return abs(t - x);
	}

	private static int LDDecode(int xi0, int xi1, int xi2, int xi3)
	{
	  int t;

	  t  = g(xi0);
	  t += g(xi1);
	  t += g(xi2);
	  t += g(xi3);

	  t -= 8*PARAM_Q;
	  t >>= 31;
	  return t&1;
	}

	private static void helprec(Poly c, Poly v, byte[] seed, byte nonce)
	{
	  int[] v0 = new int [8];
	  int v_tmp0,v_tmp1,v_tmp2,v_tmp3;
	  int k;
	  int rbit;
	  byte[] rand = new byte [32];
	  int i;

	  try {
		  crypto_stream_chacha20(rand,0,32,((long)nonce) << 56,seed);
	
		  for(i=0; i<256; i++)
		  {
		    rbit = (rand[i>>3] >> (i&7)) & 1;
	
		    k  = f(v0,0, v0,4, 8*v.coeffs[  0+i] + 4*rbit);
		    k += f(v0,1, v0,5, 8*v.coeffs[256+i] + 4*rbit);
		    k += f(v0,2, v0,6, 8*v.coeffs[512+i] + 4*rbit);
		    k += f(v0,3, v0,7, 8*v.coeffs[768+i] + 4*rbit);
	
		    k = (2*PARAM_Q-1-k) >> 31;
	
		    v_tmp0 = ((~k) & v0[0]) ^ (k & v0[4]);
		    v_tmp1 = ((~k) & v0[1]) ^ (k & v0[5]);
		    v_tmp2 = ((~k) & v0[2]) ^ (k & v0[6]);
		    v_tmp3 = ((~k) & v0[3]) ^ (k & v0[7]);
	
		    c.coeffs[  0+i] = (char)((v_tmp0 -   v_tmp3) & 3);
		    c.coeffs[256+i] = (char)((v_tmp1 -   v_tmp3) & 3);
		    c.coeffs[512+i] = (char)((v_tmp2 -   v_tmp3) & 3);
		    c.coeffs[768+i] = (char)((   -k  + 2*v_tmp3) & 3);
		  }
	  } finally {
		  Arrays.fill(v0, 0);
		  Arrays.fill(rand, (byte)0);
	  }
	}

	private static void rec(byte[] key, Poly v, Poly c)
	{
	  int i;
	  int tmp0,tmp1,tmp2,tmp3;

	  for(i=0;i<32;i++)
	    key[i] = 0;

	  for(i=0; i<256; i++)
	  {
		char c768 = c.coeffs[768+i];
	    tmp0 = 16*PARAM_Q + 8*(int)v.coeffs[  0+i] - PARAM_Q * (2*c.coeffs[  0+i]+c768);
	    tmp1 = 16*PARAM_Q + 8*(int)v.coeffs[256+i] - PARAM_Q * (2*c.coeffs[256+i]+c768);
	    tmp2 = 16*PARAM_Q + 8*(int)v.coeffs[512+i] - PARAM_Q * (2*c.coeffs[512+i]+c768);
	    tmp3 = 16*PARAM_Q + 8*(int)v.coeffs[768+i] - PARAM_Q * (                  c768);

	    key[i>>3] |= LDDecode(tmp0, tmp1, tmp2, tmp3) << (i & 7);
	  }
	}
	
	// -------------- ntt.c --------------

	private static final int bitrev_table_combined[/*496*/] = {
		524289,262146,786435,131076,655365,393222,917511,65544,
		589833,327690,851979,196620,720909,458766,983055,32784,
		557073,294930,819219,163860,688149,426006,950295,98328,
		622617,360474,884763,229404,753693,491550,1015839,540705,
		278562,802851,147492,671781,409638,933927,81960,606249,
		344106,868395,213036,737325,475182,999471,573489,311346,
		835635,180276,704565,442422,966711,114744,639033,376890,
		901179,245820,770109,507966,1032255,532545,270402,794691,
		139332,663621,401478,925767,598089,335946,860235,204876,
		729165,467022,991311,565329,303186,827475,172116,696405,
		434262,958551,106584,630873,368730,893019,237660,761949,
		499806,1024095,548961,286818,811107,155748,680037,417894,
		942183,614505,352362,876651,221292,745581,483438,1007727,
		581745,319602,843891,188532,712821,450678,974967,647289,
		385146,909435,254076,778365,516222,1040511,528513,266370,
		790659,659589,397446,921735,594057,331914,856203,200844,
		725133,462990,987279,561297,299154,823443,168084,692373,
		430230,954519,626841,364698,888987,233628,757917,495774,
		1020063,544929,282786,807075,676005,413862,938151,610473,
		348330,872619,217260,741549,479406,1003695,577713,315570,
		839859,708789,446646,970935,643257,381114,905403,250044,
		774333,512190,1036479,536769,274626,798915,667845,405702,
		929991,602313,340170,864459,733389,471246,995535,569553,
		307410,831699,700629,438486,962775,635097,372954,897243,
		241884,766173,504030,1028319,553185,291042,815331,684261,
		422118,946407,618729,356586,880875,749805,487662,1011951,
		585969,323826,848115,717045,454902,979191,651513,389370,
		913659,782589,520446,1044735,526593,788739,657669,395526,
		919815,592137,329994,854283,723213,461070,985359,559377,
		297234,821523,690453,428310,952599,624921,362778,887067,
		755997,493854,1018143,543009,805155,674085,411942,936231,
		608553,346410,870699,739629,477486,1001775,575793,837939,
		706869,444726,969015,641337,379194,903483,772413,510270,
		1034559,534849,796995,665925,403782,928071,600393,862539,
		731469,469326,993615,567633,829779,698709,436566,960855,
		633177,371034,895323,764253,502110,1026399,551265,813411,
		682341,420198,944487,616809,878955,747885,485742,1010031,
		584049,846195,715125,452982,977271,649593,911739,780669,
		518526,1042815,530817,792963,661893,924039,596361,858507,
		727437,465294,989583,563601,825747,694677,432534,956823,
		629145,891291,760221,498078,1022367,547233,809379,678309,
		940455,612777,874923,743853,481710,1005999,580017,842163,
		711093,973239,645561,907707,776637,514494,1038783,539073,
		801219,670149,932295,604617,866763,735693,997839,571857,
		834003,702933,965079,637401,899547,768477,506334,1030623,
		555489,817635,686565,948711,621033,883179,752109,1014255,
		588273,850419,719349,981495,653817,915963,784893,1047039,
		787971,656901,919047,591369,853515,722445,984591,558609,
		820755,689685,951831,624153,886299,755229,1017375,804387,
		673317,935463,607785,869931,738861,1001007,837171,706101,
		968247,640569,902715,771645,1033791,796227,665157,927303,
		861771,730701,992847,829011,697941,960087,632409,894555,
		763485,1025631,812643,681573,943719,878187,747117,1009263,
		845427,714357,976503,910971,779901,1042047,792195,923271,
		857739,726669,988815,824979,693909,956055,890523,759453,
		1021599,808611,939687,874155,743085,1005231,841395,972471,
		906939,775869,1038015,800451,931527,865995,997071,833235,
		964311,898779,767709,1029855,816867,947943,882411,1013487,
		849651,980727,915195,1046271,921351,855819,986895,823059,
		954135,888603,1019679,937767,872235,1003311,970551,905019,
		1036095,929607,995151,962391,896859,1027935,946023,1011567,
		978807,1044351,991119,958359,1023903,1007535,1040319,1032159
	};

	// Modified version of bitrev_vector() from the C reference code
	// that reduces the number of array bounds checks on the bitrev_table
	// from 1024 to 496.  The values in the combined table are encoded
	// as (i + (r * PARAM_N)) where i and r are the indices to swap.
	// The pseudo-code to generate this combined table is:
    //     p = 0;
    //     for (i = 0; i < PARAM_N; i++) {
    //         r = bitrev_table[i];
    //         if (i < r)
    //             bitrev_table_combined[p++] = i + (r * PARAM_N);
    //     }
	private static void bitrev_vector(char[] poly)
	{
	    int i,r,p;
	    char tmp;

	    for(p = 0; p < 496; ++p)
	    {
	    	int indices = bitrev_table_combined[p];
	    	i = indices & 0x03FF;
	    	r = indices >> 10;
        	tmp = poly[i];
        	poly[i] = poly[r];
        	poly[r] = tmp;
	    }
	}

	private static void mul_coefficients(char[] poly, char[] factors)
	{
	    int i;

	    for(i = 0; i < PARAM_N; i++)
	      poly[i] = (char)montgomery_reduce((poly[i] * factors[i]));
	}
	
	/* GS_bo_to_no; omegas need to be in Montgomery domain */
	private static void ntt_global(char[] a, char[] omega)
	{
	  int i, start, j, jTwiddle, distance;
	  char temp, W;


	  for(i=0;i<10;i+=2)
	  {
	    // Even level
	    distance = (1<<i);
	    for(start = 0; start < distance;start++)
	    {
	      jTwiddle = 0;
	      for(j=start;j<PARAM_N-1;j+=2*distance)
	      {
	        W = omega[jTwiddle++];
	        temp = a[j];
	        a[j] = (char)(temp + a[j + distance]); // Omit reduction (be lazy)
	        a[j + distance] = (char)montgomery_reduce((W * ((int)temp + 3*PARAM_Q - a[j + distance])));
	      }
	    }

	    // Odd level
	    distance <<= 1;
	    for(start = 0; start < distance;start++)
	    {
	      jTwiddle = 0;
	      for(j=start;j<PARAM_N-1;j+=2*distance)
	      {
	        W = omega[jTwiddle++];
	        temp = a[j];
	        a[j] = (char)barrett_reduce((temp + a[j + distance]));
	        a[j + distance] = (char)montgomery_reduce((W * ((int)temp + 3*PARAM_Q - a[j + distance])));
	      }
	    }
	  }
	}

	// -------------- fips202.c --------------

	/* Based on the public domain implementation in
	 * crypto_hash/keccakc512/simple/ from http://bench.cr.yp.to/supercop.html
	 * by Ronny Van Keer 
	 * and the public domain "TweetFips202" implementation
	 * from https://twitter.com/tweetfips202
	 * by Gilles Van Assche, Daniel J. Bernstein, and Peter Schwabe */

	private static long ROL(long a, int offset)
	{
		return (a << offset) ^ (a >>> (64 - offset));
	}

	private static long load64(byte[] x, int offset)
	{
	  long r = 0;

	  for (int i = 0; i < 8; ++i) {
	    r |= ((long)(x[offset+i] & 0xff)) << (8 * i);
	  }
	  return r;
	}

	private static void store64(byte[] x, int offset, long u)
	{
	  int i;

	  for(i=0; i<8; ++i) {
	    x[offset+i] = (byte)u;
	    u >>= 8;
	  }
	}

	private static final long[] KeccakF_RoundConstants =
		{
		    0x0000000000000001L,
		    0x0000000000008082L,
		    0x800000000000808aL,
		    0x8000000080008000L,
		    0x000000000000808bL,
		    0x0000000080000001L,
		    0x8000000080008081L,
		    0x8000000000008009L,
		    0x000000000000008aL,
		    0x0000000000000088L,
		    0x0000000080008009L,
		    0x000000008000000aL,
		    0x000000008000808bL,
		    0x800000000000008bL,
		    0x8000000000008089L,
		    0x8000000000008003L,
		    0x8000000000008002L,
		    0x8000000000000080L,
		    0x000000000000800aL,
		    0x800000008000000aL,
		    0x8000000080008081L,
		    0x8000000000008080L,
		    0x0000000080000001L,
		    0x8000000080008008L
		};

	
	private static void KeccakF1600_StatePermute(long[] state)
	{
	    int round;

        long Aba, Abe, Abi, Abo, Abu;
        long Aga, Age, Agi, Ago, Agu;
        long Aka, Ake, Aki, Ako, Aku;
        long Ama, Ame, Ami, Amo, Amu;
        long Asa, Ase, Asi, Aso, Asu;
        long BCa, BCe, BCi, BCo, BCu;
        long Da, De, Di, Do, Du;
        long Eba, Ebe, Ebi, Ebo, Ebu;
        long Ega, Ege, Egi, Ego, Egu;
        long Eka, Eke, Eki, Eko, Eku;
        long Ema, Eme, Emi, Emo, Emu;
        long Esa, Ese, Esi, Eso, Esu;

        //copyFromState(A, state)
        Aba = state[ 0];
        Abe = state[ 1];
        Abi = state[ 2];
        Abo = state[ 3];
        Abu = state[ 4];
        Aga = state[ 5];
        Age = state[ 6];
        Agi = state[ 7];
        Ago = state[ 8];
        Agu = state[ 9];
        Aka = state[10];
        Ake = state[11];
        Aki = state[12];
        Ako = state[13];
        Aku = state[14];
        Ama = state[15];
        Ame = state[16];
        Ami = state[17];
        Amo = state[18];
        Amu = state[19];
        Asa = state[20];
        Ase = state[21];
        Asi = state[22];
        Aso = state[23];
        Asu = state[24];

        for( round = 0; round < 24; round += 2 )
        {
            //    prepareTheta
            BCa = Aba^Aga^Aka^Ama^Asa;
            BCe = Abe^Age^Ake^Ame^Ase;
            BCi = Abi^Agi^Aki^Ami^Asi;
            BCo = Abo^Ago^Ako^Amo^Aso;
            BCu = Abu^Agu^Aku^Amu^Asu;

            //thetaRhoPiChiIotaPrepareTheta(round  , A, E)
            Da = BCu^ROL(BCe, 1);
            De = BCa^ROL(BCi, 1);
            Di = BCe^ROL(BCo, 1);
            Do = BCi^ROL(BCu, 1);
            Du = BCo^ROL(BCa, 1);

            Aba ^= Da;
            BCa = Aba;
            Age ^= De;
            BCe = ROL(Age, 44);
            Aki ^= Di;
            BCi = ROL(Aki, 43);
            Amo ^= Do;
            BCo = ROL(Amo, 21);
            Asu ^= Du;
            BCu = ROL(Asu, 14);
            Eba =   BCa ^((~BCe)&  BCi );
            Eba ^= KeccakF_RoundConstants[round];
            Ebe =   BCe ^((~BCi)&  BCo );
            Ebi =   BCi ^((~BCo)&  BCu );
            Ebo =   BCo ^((~BCu)&  BCa );
            Ebu =   BCu ^((~BCa)&  BCe );

            Abo ^= Do;
            BCa = ROL(Abo, 28);
            Agu ^= Du;
            BCe = ROL(Agu, 20);
            Aka ^= Da;
            BCi = ROL(Aka,  3);
            Ame ^= De;
            BCo = ROL(Ame, 45);
            Asi ^= Di;
            BCu = ROL(Asi, 61);
            Ega =   BCa ^((~BCe)&  BCi );
            Ege =   BCe ^((~BCi)&  BCo );
            Egi =   BCi ^((~BCo)&  BCu );
            Ego =   BCo ^((~BCu)&  BCa );
            Egu =   BCu ^((~BCa)&  BCe );

            Abe ^= De;
            BCa = ROL(Abe,  1);
            Agi ^= Di;
            BCe = ROL(Agi,  6);
            Ako ^= Do;
            BCi = ROL(Ako, 25);
            Amu ^= Du;
            BCo = ROL(Amu,  8);
            Asa ^= Da;
            BCu = ROL(Asa, 18);
            Eka =   BCa ^((~BCe)&  BCi );
            Eke =   BCe ^((~BCi)&  BCo );
            Eki =   BCi ^((~BCo)&  BCu );
            Eko =   BCo ^((~BCu)&  BCa );
            Eku =   BCu ^((~BCa)&  BCe );

            Abu ^= Du;
            BCa = ROL(Abu, 27);
            Aga ^= Da;
            BCe = ROL(Aga, 36);
            Ake ^= De;
            BCi = ROL(Ake, 10);
            Ami ^= Di;
            BCo = ROL(Ami, 15);
            Aso ^= Do;
            BCu = ROL(Aso, 56);
            Ema =   BCa ^((~BCe)&  BCi );
            Eme =   BCe ^((~BCi)&  BCo );
            Emi =   BCi ^((~BCo)&  BCu );
            Emo =   BCo ^((~BCu)&  BCa );
            Emu =   BCu ^((~BCa)&  BCe );

            Abi ^= Di;
            BCa = ROL(Abi, 62);
            Ago ^= Do;
            BCe = ROL(Ago, 55);
            Aku ^= Du;
            BCi = ROL(Aku, 39);
            Ama ^= Da;
            BCo = ROL(Ama, 41);
            Ase ^= De;
            BCu = ROL(Ase,  2);
            Esa =   BCa ^((~BCe)&  BCi );
            Ese =   BCe ^((~BCi)&  BCo );
            Esi =   BCi ^((~BCo)&  BCu );
            Eso =   BCo ^((~BCu)&  BCa );
            Esu =   BCu ^((~BCa)&  BCe );

            //    prepareTheta
            BCa = Eba^Ega^Eka^Ema^Esa;
            BCe = Ebe^Ege^Eke^Eme^Ese;
            BCi = Ebi^Egi^Eki^Emi^Esi;
            BCo = Ebo^Ego^Eko^Emo^Eso;
            BCu = Ebu^Egu^Eku^Emu^Esu;

            //thetaRhoPiChiIotaPrepareTheta(round+1, E, A)
            Da = BCu^ROL(BCe, 1);
            De = BCa^ROL(BCi, 1);
            Di = BCe^ROL(BCo, 1);
            Do = BCi^ROL(BCu, 1);
            Du = BCo^ROL(BCa, 1);

            Eba ^= Da;
            BCa = Eba;
            Ege ^= De;
            BCe = ROL(Ege, 44);
            Eki ^= Di;
            BCi = ROL(Eki, 43);
            Emo ^= Do;
            BCo = ROL(Emo, 21);
            Esu ^= Du;
            BCu = ROL(Esu, 14);
            Aba =   BCa ^((~BCe)&  BCi );
            Aba ^= KeccakF_RoundConstants[round+1];
            Abe =   BCe ^((~BCi)&  BCo );
            Abi =   BCi ^((~BCo)&  BCu );
            Abo =   BCo ^((~BCu)&  BCa );
            Abu =   BCu ^((~BCa)&  BCe );

            Ebo ^= Do;
            BCa = ROL(Ebo, 28);
            Egu ^= Du;
            BCe = ROL(Egu, 20);
            Eka ^= Da;
            BCi = ROL(Eka, 3);
            Eme ^= De;
            BCo = ROL(Eme, 45);
            Esi ^= Di;
            BCu = ROL(Esi, 61);
            Aga =   BCa ^((~BCe)&  BCi );
            Age =   BCe ^((~BCi)&  BCo );
            Agi =   BCi ^((~BCo)&  BCu );
            Ago =   BCo ^((~BCu)&  BCa );
            Agu =   BCu ^((~BCa)&  BCe );

            Ebe ^= De;
            BCa = ROL(Ebe, 1);
            Egi ^= Di;
            BCe = ROL(Egi, 6);
            Eko ^= Do;
            BCi = ROL(Eko, 25);
            Emu ^= Du;
            BCo = ROL(Emu, 8);
            Esa ^= Da;
            BCu = ROL(Esa, 18);
            Aka =   BCa ^((~BCe)&  BCi );
            Ake =   BCe ^((~BCi)&  BCo );
            Aki =   BCi ^((~BCo)&  BCu );
            Ako =   BCo ^((~BCu)&  BCa );
            Aku =   BCu ^((~BCa)&  BCe );

            Ebu ^= Du;
            BCa = ROL(Ebu, 27);
            Ega ^= Da;
            BCe = ROL(Ega, 36);
            Eke ^= De;
            BCi = ROL(Eke, 10);
            Emi ^= Di;
            BCo = ROL(Emi, 15);
            Eso ^= Do;
            BCu = ROL(Eso, 56);
            Ama =   BCa ^((~BCe)&  BCi );
            Ame =   BCe ^((~BCi)&  BCo );
            Ami =   BCi ^((~BCo)&  BCu );
            Amo =   BCo ^((~BCu)&  BCa );
            Amu =   BCu ^((~BCa)&  BCe );

            Ebi ^= Di;
            BCa = ROL(Ebi, 62);
            Ego ^= Do;
            BCe = ROL(Ego, 55);
            Eku ^= Du;
            BCi = ROL(Eku, 39);
            Ema ^= Da;
            BCo = ROL(Ema, 41);
            Ese ^= De;
            BCu = ROL(Ese, 2);
            Asa =   BCa ^((~BCe)&  BCi );
            Ase =   BCe ^((~BCi)&  BCo );
            Asi =   BCi ^((~BCo)&  BCu );
            Aso =   BCo ^((~BCu)&  BCa );
            Asu =   BCu ^((~BCa)&  BCe );
        }

        //copyToState(state, A)
        state[ 0] = Aba;
        state[ 1] = Abe;
        state[ 2] = Abi;
        state[ 3] = Abo;
        state[ 4] = Abu;
        state[ 5] = Aga;
        state[ 6] = Age;
        state[ 7] = Agi;
        state[ 8] = Ago;
        state[ 9] = Agu;
        state[10] = Aka;
        state[11] = Ake;
        state[12] = Aki;
        state[13] = Ako;
        state[14] = Aku;
        state[15] = Ama;
        state[16] = Ame;
        state[17] = Ami;
        state[18] = Amo;
        state[19] = Amu;
        state[20] = Asa;
        state[21] = Ase;
        state[22] = Asi;
        state[23] = Aso;
        state[24] = Asu;
	}

	private static void keccak_absorb(long[] s, int r, byte[] m, int offset, int mlen, byte p)
	{
	  int i;
	  byte[] t = new byte [200];

	  try {
		  for (i = 0; i < 25; ++i)
		    s[i] = 0;

		  while (mlen >= r)
		  {
		    for (i = 0; i < r / 8; ++i)
		      s[i] ^= load64(m, offset + 8 * i);

		    KeccakF1600_StatePermute(s);
		    mlen -= r;
		    offset += r;
		  }

		  for (i = 0; i < r; ++i)
		    t[i] = 0;
		  for (i = 0; i < mlen; ++i)
		    t[i] = m[offset + i];
		  t[i] = p;
		  t[r - 1] |= 128;
		  for (i = 0; i < r / 8; ++i)
		    s[i] ^= load64(t, 8 * i);
	  } finally {
		  Arrays.fill(t, (byte)0);
	  }
	}

	private static void keccak_squeezeblocks(byte[] h, int offset, int nblocks, long [] s, int r)
	{
	  int i;
