}

			// x[2] = x[0] - x[1] * q (mod b0)
			SubMulMod( x1, q, x0, b0, x2 );

			// {q, b} = b / a
			udivide( b, a, q, b );

			// if b = 0, return x[2]
			if ( isZero( b ) )
			{
				set ( ap, x2 )

					;
				return ;
			}

			// x[0] = x[1] - x[2] * q (mod b0)
			SubMulMod( x2, q, x1, b0, x0 );
		}
	}

	// indirectly computes  n' s.t. 1 = r' r - n' n = GCD(r, n)
	BIGONETYPE void computeNRinverse( T &n0, T &np )
	{
		BIGDOUBLESIZE( T, r );
		BIGDOUBLESIZE( T, n );

		// r' = (1 + n' n) / r
		computeRinverse( n0, np );

		// n' = (r' r - 1) / n
		sethigh( r, np ); // special case of r = 2^k
		decrement( r );
		usetlow( n, n0 );
		udivide( r, n, n, r );
		takelow( np, n );
	}

	/*
	// indirectly computes n' s.t. 1 = r' r - n' n = GCD(r, n)
	BIGONETYPE void computeNRinverse2(T &n0, T &np)
	{
	BIGDOUBLESIZE(T, r);
	BIGDOUBLESIZE(T, n);

	// r' = (1 + n' n) / r
	computeRinverse2(n0, np);

	// n' = (r' r - 1) / n
	sethigh(r, np); // special case of r = 2^k
	decrement(r);
	usetlow(n, n0);
	udivide(r, n, n, r);
	takelow(np, n);
	}
	*/
	// Montgomery product u = a * b (mod n)
	BIGONETYPE void MonPro( T &ap, T &bp, T &n, T &np, T &u_out )
	{
		BIGDOUBLESIZE( T, t );
		BIGDOUBLESIZE( T, u );
		T m;

		// t = a' b'
		umultiply( ap, bp, t );

		// m = (low half of t)*np (mod r)
		takelow( m, t );
		umultiply( m, np );

		// u = (t + m*n), u_out = u / r = high half of u
		umultiply( m, n, u );

		add ( u, t )

			;
		takehigh( u_out, u );

		// if u >= n, return u - n, else u
		if ( ugreaterOrEqual( u_out, n ) )
			subtract( u_out, n );
	}

	// indirectly calculates x = M^e (mod n)
	BIGONETYPE void MonModExp( T &x, T &M, T &e, T &n, T &np, T &xp0 )
	{
		// x' = xp0

		set ( x, xp0 )

			;

		// find M' = M r (mod n)
		BIGDOUBLESIZE( T, dM );

		BIGDOUBLESIZE( T, dn );

		T Mp;

		sethigh( dM, M ); // dM = M r

		usetlow( dn, n );

		umodulo( dM, dn, dM ); // dM = dM (mod n)

		takelow( Mp, dM ); // M' = M r (mod n)

		/*	i may be wrong, but it seems to me that the squaring
		results in a constant until we hit the first set bit
		this could save a lot of time, but it needs to be proven
		*/

		s32 ii, bc;

		word e_i;

		// for i = k - 1 down to 0 do
		for ( ii = BIGWORDCOUNT( T ) - 1; ii >= 0; --ii )
		{
			e_i = e[ ii ];
			bc = WORDBITS;

			while ( bc-- )
			{
				// if e_i = 1, x = MonPro(M', x')

				if ( e_i & WORDHIGHBIT )
					goto start_squaring;

				e_i <<= 1;
			}
		}

		for ( ; ii >= 0; --ii )
		{
			e_i = e[ ii ];
			bc = WORDBITS;

			while ( bc-- )
			{
				// x' = MonPro(x', x')
				MonPro( x, x, n, np, x );

				// if e_i = 1, x = MonPro(M', x')

				if ( e_i & WORDHIGHBIT )
				{

start_squaring:
					MonPro( Mp, x, n, np, x );
				}

				e_i <<= 1;
			}
		}

		// x = MonPro(x', 1)
		T one;

		usetw( one, 1 );

		MonPro( x, one, n, np, x );
	}

	// indirectly calculates x = C ^ d (mod n) using the Chinese Remainder Thm
#ifdef _MSC_VER
	#pragma warning( disable : 4100 ) // warning C4100: <variable name> : unreferenced formal parameter
#endif
	BIGTWOTYPES void CRTModExp( Bigger &x, Bigger &C, Bigger &d, T &p, T &q, T &pInverse, T &pnp, T &pxp, T &qnp, T &qxp )
	{
		// d1 = d mod (p - 1)
		Bigger dd1;
		T d1;
		usetlow( dd1, p );
		decrement( dd1 );
		umodulo( d, dd1, dd1 );
		takelow( d1, dd1 );

		// M1 = C1^d1 (mod p)
		Bigger dp, dC1;
		T M1, C1;
		usetlow( dp, p );
		umodulo( C, dp, dC1 );
		takelow( C1, dC1 );
		simpleModExp( C1, d1, p, M1 );
		//MonModExp(M1, C1, d1, p, pnp, pxp);

		// d2 = d mod (q - 1)
		Bigger dd2;
		T d2;
		usetlow( dd2, q );
		decrement( dd2 );
		umodulo( d, dd2, dd2 );
		takelow( d2, dd2 );

		// M2 = C2^d2 (mod q)
		Bigger dq, dC2;
		T M2, C2;
		usetlow( dq, q );
		umodulo( C, dq, dC2 );
		takelow( C2, dC2 );
		simpleModExp( C2, d2, q, M2 );
		//MonModExp(M2, C2, d2, q, qnp, qxp);

		// x = M1 + p * ((M2 - M1)(p^-1 mod q) mod q)

		if ( ugreater( M2, M1 ) )
		{
			subtract( M2, M1 );
		}

		else
		{
			subtract( M2, M1 );

			add ( M2, q )

				;
		}

		// x = M1 + p * (( M2 )(p^-1 mod q) mod q)

		umultiply( M2, pInverse, x );

		// x = M1 + p * (( x ) mod q)

		umodulo( x, dq, x );

		// x = M1 + p * ( x )

		umultiply( x, dp );

		// x = M1 + ( x )

		Bigger dM1;

		usetlow( dM1, M1 );

		// x = ( dM1 ) + ( x )

		add ( x, dM1 )

			;
	}

	// generates a suitable public exponent s.t. 4 < e << phi, GCD(e, phi) = 1
	BIGONETYPE void computePublicExponent( T &phi, T &e )
	{
		T r, one, two;
		usetw( one, 1 );
		usetw( two, 2 );
		usetw( e, 65537 - 2 );

		if ( ugreater( e, phi ) )
			usetw( e, 5 - 2 );

		do
		{
			add ( e, two )

				;

			GCD( phi, e, r );
		}

		while ( !equal( r, one ) );
	}

	// directly computes private exponent
	BIGONETYPE void computePrivateExponent( T &e, T &phi, T &d )
	{
		// d = e^-1 (mod phi), 1 < e << phi
		computeModularInverse( e, phi, d );
	}

#ifdef RSASUPPORTGENPRIME

	static const u16 PRIME_TABLE[ 256 ] =
	{
		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, 257, 263, 269,
			271, 277, 281, 283, 293, 307, 311, 313,
			317, 331, 337, 347, 349, 353, 359, 367,
			373, 379, 383, 389, 397, 401, 409, 419,
			421, 431, 433, 439, 443, 449, 457, 461,
			463, 467, 479, 487, 491, 499, 503, 509,
			521, 523, 541, 547, 557, 563, 569, 571,
			577, 587, 593, 599, 601, 607, 613, 617,
			619, 631, 641, 643, 647, 653, 659, 661,
			673, 677, 683, 691, 701, 709, 719, 727,
			733, 739, 743, 751, 757, 761, 769, 773,
			787, 797, 809, 811, 821, 823, 827, 829,
			839, 853, 857, 859, 863, 877, 881, 883,
			887, 907, 911, 919, 929, 937, 941, 947,
			953, 967, 971, 977, 983, 991, 997, 1009,
			1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051,
			1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103,
			1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171,
			1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229,
			1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289,
			1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327,
			1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427,
			1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471,
			1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523,
			1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579,
			1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621
	};

	/*modified Rabin-Miller primality test (added small primes)

	When picking a value for insurance, note that the probability of failure
	of the test to detect a composite number is at most 4^(-insurance), so:
	insurance max. probability of failure
	3   1.56%
	4   0.39%
	5   0.098% <-- default
	6   0.024%
	...
	*/
	BIGONETYPE bool RabinMillerPrimalityTest( T &n, u32 insurance )
	{
		// check divisibility by small primes <= 1621 (speeds up computation)
		T temp;

		for ( u32 ii = 0; ii < 256; ++ii )
		{
			usetw( temp, PRIME_TABLE[ ii++ ] );

			umodulo( n, temp, temp );

			if ( isZero( temp ) )
				return false;
		}

		// n1 = n - 1
		T n1;

		set ( n1, n )

			;
		decrement( n1 );

		// write r 2^s = n - 1, r is odd
		T r;

		u32 s = 0;

		set ( r, n1 )

			;
		while ( !( r[ 0 ] & 1 ) )
		{
			ushiftRight1( r );
			++s;
		}

		// one = 1
		T one;

		usetw( one, 1 );

		// cache n -> dn
		BIGDOUBLESIZE( T, dy );

		BIGDOUBLESIZE( T, dn );

		usetlow( dn, n );

		while ( insurance-- )
		{
			// choose random integer a s.t. 1 < a < n - 1
			T a;
			int index;

			for ( index = 0; index < (int) sizeof( a ) / (int) sizeof( a[ 0 ] ); index++ )
				a[ index ] = randomMT();

			umodulo( a, n1, a );

			// compute y = a ^ r (mod n)
			T y;

			simpleModExp( a, r, n, y );

			if ( !equal( y, one ) && !equal( y, n1 ) )
			{
				u32 j = s;

				while ( ( j-- > 1 ) && !equal( y, n1 ) )
				{
					umultiply( y, y, dy );
					umodulo( dy, dn, dy );
					takelow( y, dy );

					if ( equal( y, one ) )
						return false;
				}

				if ( !equal( y, n1 ) )
					return false;
			}
		}

		return true;
	}

	// generates a strong pseudo-prime
	BIGONETYPE void generateStrongPseudoPrime( T &n )
	{
		do
		{
			int index;

			for ( index = 0; index < (int) sizeof( n ) / (int) sizeof( n[ 0 ] ); index++ )
				n[ index ] = randomMT();

			n[ BIGWORDCOUNT( T ) - 1 ] |= WORDHIGHBIT;

			//n[BIGWORDCOUNT(T) - 1] &= ~WORDHIGHBIT; n[BIGWORDCOUNT(T) - 1] |= WORDHIGHBIT >> 1;
			n[ 0 ] |= 1;
		}

		while ( !RabinMillerPrimalityTest( n, 5 ) );
	}

#endif // RSASUPPORTGENPRIME


	//////// RSACrypt class ////////

	BIGONETYPE class RAK_DLL_EXPORT RSACrypt
	{
		// public key
		T e, n;
		T np, xp;

		// private key
		bool factorsAvailable;
		T d, phi;
		BIGHALFSIZE( T, p );
		BIGHALFSIZE( T, pnp );
		BIGHALFSIZE( T, pxp );
		BIGHALFSIZE( T, q );
		BIGHALFSIZE( T, qnp );
		BIGHALFSIZE( T, qxp );
		BIGHALFSIZE( T, pInverse );

	public:
		RSACrypt()
		{
			reset();
		}

		~RSACrypt()
		{
			reset();
		}

	public:
		void reset()
		{
			zero( d );
			zero( p );
			zero( q );
			zero( pInverse );
			factorsAvailable = false;
		}

#ifdef RSASUPPORTGENPRIME

		void generateKeys()
		{
			BIGHALFSIZE( T, p0 );
			BIGHALFSIZE( T, q0 );

			generateStrongPseudoPrime( p0 );
			generateStrongPseudoPrime( q0 );

			setPrivateKey( p0, q0 );
		}

#endif // RSASUPPORTGENPRIME

		BIGSMALLTYPE void setPrivateKey( Smaller &c_p, Smaller &c_q )
		{
			factorsAvailable = true;

			// re-order factors s.t. q > p

			if ( ugreater( c_p, c_q ) )
			{
				set ( q, c_p )

					;
				set ( p, c_q )

					;
			}

			else
			{
				set ( p, c_p )

					;
				set ( q, c_q )

					;
			}

			// phi = (p - 1)(q - 1)
			BIGHALFSIZE( T, p1 );

			BIGHALFSIZE( T, q1 );

			set ( p1, p )

				;
			decrement( p1 );

			set ( q1, q )

				;
			decrement( q1 );

			umultiply( p1, q1, phi );

			// compute e
			computePublicExponent( phi, e );

			// compute d
			computePrivateExponent( e, phi, d );

			// compute p^-1 mod q
			computeModularInverse( p, q, pInverse );

			// compute n = pq
			umultiply( p, q, n );

			// find n'
			computeNRinverse( n, np );

			// x' = 1*r (mod n)
			rModn( n, xp );

			// find pn'
			computeNRinverse( p, pnp );

			//   computeNRinverse2(p, pnp);

			// px' = 1*r (mod p)
			rModn( p, pxp );

			// find qn'
			computeNRinverse( q, qnp );

			// qx' = 1*r (mod q)
			rModn( q, qxp );
		}

		void setPublicKey( u32 c_e, T &c_n )
		{
			reset(); // in case we knew a private key

			usetw( e, c_e );

			set ( n, c_n )

				;

			// find n'
			computeNRinverse( n, np );

			// x' = 1*r (mod n)
			rModn( n, xp );
		}

	public:
		void getPublicKey( u32 &c_e, T &c_n )
		{
			c_e = e[ 0 ];

			set ( c_n, n )

				;
		}

		BIGSMALLTYPE void getPrivateKey( Smaller &c_p, Smaller &c_q )
		{
			set ( c_p, p )

				;
			set ( c_q, q )

				;
		}

	public:
		void encrypt( T &M, T &x )
		{
			if ( factorsAvailable )
				CRTModExp( x, M, e, p, q, pInverse, pnp, pxp, qnp, qxp );
			else
				simpleModExp( M, e, n, x );
		}

		void decrypt( T &C, T &x )
		{
			if ( factorsAvailable )
				CRTModExp( x, C, d, p, q, pInverse, pnp, pxp, qnp, qxp );
		}
	};
}

#ifdef _MSC_VER
#pragma warning( pop )
#endif

#endif // #if !defined(_COMPATIBILITY_1)

#endif // RSACRYPT_H