/*

* Copyright (c) 2003-2009 Nokia Corporation and/or its subsidiary(-ies).

* All rights reserved.

* This component and the accompanying materials are made available

* under the terms of the License "Eclipse Public License v1.0"

* which accompanies this distribution, and is available

* at the URL "http://www.eclipse.org/legal/epl-v10.html".

*

* Initial Contributors:

* Nokia Corporation - initial contribution.

*

* Contributors:

*

* Description: 

*

*/

#include "words.h"

#include "algorithms.h"

word Add(word *C, const word *A, const word *B, unsigned int N)

{

	assert (N%2 == 0);

	word carry = 0;

	for (unsigned int i = 0; i < N; i+=2)

	{

		dword u = (dword) carry + A[i] + B[i];

		C[i] = LOW_WORD(u);

		u = (dword) HIGH_WORD(u) + A[i+1] + B[i+1];

		C[i+1] = LOW_WORD(u);

		carry = HIGH_WORD(u);

	}

	return carry;

}

word Subtract(word *C, const word *A, const word *B, unsigned int N)

{

	assert (N%2 == 0);

	word borrow=0;

	for (unsigned i = 0; i < N; i+=2)

	{

		dword u = (dword) A[i] - B[i] - borrow;

		C[i] = LOW_WORD(u);

		u = (dword) A[i+1] - B[i+1] - (word)(0-HIGH_WORD(u));

		C[i+1] = LOW_WORD(u);

		borrow = 0-HIGH_WORD(u);

	}

	return borrow;

}

int Compare(const word *A, const word *B, unsigned int N)

{

	while (N--)

		if (A[N] > B[N])

			return 1;

		else if (A[N] < B[N])

			return -1;

	return 0;

}

// It is the job of the calling code to ensure that this won't carry.

// If you aren't sure, use the next version that will tell you if you need to

// grow your integer.

// Having two of these creates ever so slightly more code but avoids having

// ifdefs all over the rest of the code checking the following type stuff which

// causes warnings in certain compilers about unused parameters in release

// builds.  We can't have that can we!

/*

Allows avoid this all over bigint.cpp and primes.cpp

ifdef _DEBUG

	TUint carry = Increment(Ptr(), Size());

	assert(!carry);

else

	Increment(Ptr(), Size())

endif

*/

void IncrementNoCarry(word *A, unsigned int N, word B)

{

	assert(N);

	word t = A[0];

	A[0] = t+B;

	if (A[0] >= t)

		return;

	for (unsigned i=1; i<N; i++)

		if (++A[i])

			return;

	assert(0);

}

word Increment(word *A, unsigned int N, word B)

{

	assert(N);

	word t = A[0];

	A[0] = t+B;

	if (A[0] >= t)

		return 0;

	for (unsigned i=1; i<N; i++)

		if (++A[i])

			return 0;

	return 1;

}

//See commments above about IncrementNoCarry

void DecrementNoCarry(word *A, unsigned int N, word B)

{

	assert(N);

	word t = A[0];

	A[0] = t-B;

	if (A[0] <= t)

		return;

	for (unsigned i=1; i<N; i++)

		if (A[i]--)

			return;

	assert(0);

}

word Decrement(word *A, unsigned int N, word B)

{

	assert(N);

	word t = A[0];

	A[0] = t-B;

	if (A[0] <= t)

		return 0;

	for (unsigned i=1; i<N; i++)

		if (A[i]--)

			return 0;

	return 1;

}

void TwosComplement(word *A, unsigned int N)

{

	Decrement(A, N);

	for (unsigned i=0; i<N; i++)

		A[i] = ~A[i];

}

static word LinearMultiply(word *C, const word *A, word B, unsigned int N)

{

	word carry=0;

	for(unsigned i=0; i<N; i++)

	{

		dword p = (dword)A[i] * B + carry;

		C[i] = LOW_WORD(p);

		carry = HIGH_WORD(p);

	}

	return carry;

}

static void AtomicMultiply(word *C, const word *A, const word *B)

{

/*

	word s;

	dword d;

	if (A1 >= A0)

		if (B0 >= B1)

		{

			s = 0;

			d = (dword)(A1-A0)*(B0-B1);

		}

		else

		{

			s = (A1-A0);

			d = (dword)s*(word)(B0-B1);

		}

	else

		if (B0 > B1)

		{

			s = (B0-B1);

			d = (word)(A1-A0)*(dword)s;

		}

		else

		{

			s = 0;

			d = (dword)(A0-A1)*(B1-B0);

		}

*/

	// this segment is the branchless equivalent of above

	word D[4] = {A[1]-A[0], A[0]-A[1], B[0]-B[1], B[1]-B[0]};

	unsigned int ai = A[1] < A[0];

	unsigned int bi = B[0] < B[1];

	unsigned int di = ai & bi;

	dword d = (dword)D[di]*D[di+2];

	D[1] = D[3] = 0;

	unsigned int si = ai + !bi;

	word s = D[si];

	dword A0B0 = (dword)A[0]*B[0];

	C[0] = LOW_WORD(A0B0);

	dword A1B1 = (dword)A[1]*B[1];

	dword t = (dword) HIGH_WORD(A0B0) + LOW_WORD(A0B0) + LOW_WORD(d) + LOW_WORD(A1B1);

	C[1] = LOW_WORD(t);

	t = A1B1 + HIGH_WORD(t) + HIGH_WORD(A0B0) + HIGH_WORD(d) + HIGH_WORD(A1B1) - s;

	C[2] = LOW_WORD(t);

	C[3] = HIGH_WORD(t);

}

static word AtomicMultiplyAdd(word *C, const word *A, const word *B)

{

	word D[4] = {A[1]-A[0], A[0]-A[1], B[0]-B[1], B[1]-B[0]};

	unsigned int ai = A[1] < A[0];

	unsigned int bi = B[0] < B[1];

	unsigned int di = ai & bi;

	dword d = (dword)D[di]*D[di+2];

	D[1] = D[3] = 0;

	unsigned int si = ai + !bi;

	word s = D[si];

	dword A0B0 = (dword)A[0]*B[0];

	dword t = A0B0 + C[0];

	C[0] = LOW_WORD(t);

	dword A1B1 = (dword)A[1]*B[1];

	t = (dword) HIGH_WORD(t) + LOW_WORD(A0B0) + LOW_WORD(d) + LOW_WORD(A1B1) + C[1];

	C[1] = LOW_WORD(t);

	t = (dword) HIGH_WORD(t) + LOW_WORD(A1B1) + HIGH_WORD(A0B0) + HIGH_WORD(d) + HIGH_WORD(A1B1) - s + C[2];

	C[2] = LOW_WORD(t);

	t = (dword) HIGH_WORD(t) + HIGH_WORD(A1B1) + C[3];

	C[3] = LOW_WORD(t);

	return HIGH_WORD(t);

}

static inline void AtomicMultiplyBottom(word *C, const word *A, const word *B)

{

	dword t = (dword)A[0]*B[0];

	C[0] = LOW_WORD(t);

	C[1] = HIGH_WORD(t) + A[0]*B[1] + A[1]*B[0];

}

#define MulAcc(x, y)								\

	p = (dword)A[x] * B[y] + c; 					\

	c = LOW_WORD(p);								\

	p = (dword)d + HIGH_WORD(p);					\

	d = LOW_WORD(p);								\

	e += HIGH_WORD(p);

#define SaveMulAcc(s, x, y) 						\

	R[s] = c;										\

	p = (dword)A[x] * B[y] + d; 					\

	c = LOW_WORD(p);								\

	p = (dword)e + HIGH_WORD(p);					\

	d = LOW_WORD(p);								\

	e = HIGH_WORD(p);

#define MulAcc1(x, y)								\

	p = (dword)A[x] * A[y] + c; 					\

	c = LOW_WORD(p);								\

	p = (dword)d + HIGH_WORD(p);					\

	d = LOW_WORD(p);								\

	e += HIGH_WORD(p);

#define SaveMulAcc1(s, x, y) 						\

	R[s] = c;										\

	p = (dword)A[x] * A[y] + d; 					\

	c = LOW_WORD(p);								\

	p = (dword)e + HIGH_WORD(p);					\

	d = LOW_WORD(p);								\

	e = HIGH_WORD(p);

#define SquAcc(x, y)								\

	p = (dword)A[x] * A[y];	\

	p = p + p + c; 					\

	c = LOW_WORD(p);								\

	p = (dword)d + HIGH_WORD(p);					\

	d = LOW_WORD(p);								\

	e += HIGH_WORD(p);

#define SaveSquAcc(s, x, y) 						\

	R[s] = c;										\

	p = (dword)A[x] * A[y];	\

	p = p + p + d; 					\

	c = LOW_WORD(p);								\

	p = (dword)e + HIGH_WORD(p);					\

	d = LOW_WORD(p);								\

	e = HIGH_WORD(p);

// VC60 workaround: MSVC 6.0 has an optimization problem that makes

// (dword)A*B where either A or B has been cast to a dword before

// very expensive. Revisit a CombaSquare4() function when this

// problem is fixed.               

// WARNING: KeithR.  05/08/03 This routine doesn't work with gcc on hardware

// either.  I've completely removed it.  It may be worth looking into sometime

// in the future.

/*#ifndef __WINS__

static void CombaSquare4(word *R, const word *A)

{

	dword p;

	word c, d, e;

	p = (dword)A[0] * A[0];

	R[0] = LOW_WORD(p);

	c = HIGH_WORD(p);

	d = e = 0;

	SquAcc(0, 1);

	SaveSquAcc(1, 2, 0);

	MulAcc1(1, 1);

	SaveSquAcc(2, 0, 3);

	SquAcc(1, 2);

	SaveSquAcc(3, 3, 1);

	MulAcc1(2, 2);

	SaveSquAcc(4, 2, 3);

	R[5] = c;

	p = (dword)A[3] * A[3] + d;

	R[6] = LOW_WORD(p);

	R[7] = e + HIGH_WORD(p);

}

#endif */

static void CombaMultiply4(word *R, const word *A, const word *B)

{

	dword p;

	word c, d, e;

	p = (dword)A[0] * B[0];

	R[0] = LOW_WORD(p);

	c = HIGH_WORD(p);

	d = e = 0;

	MulAcc(0, 1);

	MulAcc(1, 0);

	SaveMulAcc(1, 2, 0);

	MulAcc(1, 1);

	MulAcc(0, 2);

	SaveMulAcc(2, 0, 3);

	MulAcc(1, 2);

	MulAcc(2, 1);

	MulAcc(3, 0);

	SaveMulAcc(3, 3, 1);

	MulAcc(2, 2);

	MulAcc(1, 3);

	SaveMulAcc(4, 2, 3);

	MulAcc(3, 2);

	R[5] = c;

	p = (dword)A[3] * B[3] + d;

	R[6] = LOW_WORD(p);

	R[7] = e + HIGH_WORD(p);

}

static void CombaMultiply8(word *R, const word *A, const word *B)

{

	dword p;

	word c, d, e;

	p = (dword)A[0] * B[0];

	R[0] = LOW_WORD(p);

	c = HIGH_WORD(p);

	d = e = 0;

	MulAcc(0, 1);

	MulAcc(1, 0);

	SaveMulAcc(1, 2, 0);

	MulAcc(1, 1);

	MulAcc(0, 2);

	SaveMulAcc(2, 0, 3);

	MulAcc(1, 2);

	MulAcc(2, 1);

	MulAcc(3, 0);

	SaveMulAcc(3, 0, 4);

	MulAcc(1, 3);

	MulAcc(2, 2);

	MulAcc(3, 1);

	MulAcc(4, 0);

	SaveMulAcc(4, 0, 5);

	MulAcc(1, 4);

	MulAcc(2, 3);

	MulAcc(3, 2);

	MulAcc(4, 1);

	MulAcc(5, 0);

	SaveMulAcc(5, 0, 6);

	MulAcc(1, 5);

	MulAcc(2, 4);

	MulAcc(3, 3);

	MulAcc(4, 2);

	MulAcc(5, 1);

	MulAcc(6, 0);

	SaveMulAcc(6, 0, 7);

	MulAcc(1, 6);

	MulAcc(2, 5);

	MulAcc(3, 4);

	MulAcc(4, 3);

	MulAcc(5, 2);

	MulAcc(6, 1);

	MulAcc(7, 0);

	SaveMulAcc(7, 1, 7);

	MulAcc(2, 6);

	MulAcc(3, 5);

	MulAcc(4, 4);

	MulAcc(5, 3);

	MulAcc(6, 2);

	MulAcc(7, 1);

	SaveMulAcc(8, 2, 7);

	MulAcc(3, 6);

	MulAcc(4, 5);

	MulAcc(5, 4);

	MulAcc(6, 3);

	MulAcc(7, 2);

	SaveMulAcc(9, 3, 7);

	MulAcc(4, 6);

	MulAcc(5, 5);

	MulAcc(6, 4);

	MulAcc(7, 3);

	SaveMulAcc(10, 4, 7);

	MulAcc(5, 6);

	MulAcc(6, 5);

	MulAcc(7, 4);

	SaveMulAcc(11, 5, 7);

	MulAcc(6, 6);

	MulAcc(7, 5);

	SaveMulAcc(12, 6, 7);

	MulAcc(7, 6);

	R[13] = c;

	p = (dword)A[7] * B[7] + d;

	R[14] = LOW_WORD(p);

	R[15] = e + HIGH_WORD(p);

}

static void CombaMultiplyBottom4(word *R, const word *A, const word *B)

{

	dword p;

	word c, d, e;

	p = (dword)A[0] * B[0];

	R[0] = LOW_WORD(p);

	c = HIGH_WORD(p);

	d = e = 0;

	MulAcc(0, 1);

	MulAcc(1, 0);

	SaveMulAcc(1, 2, 0);

	MulAcc(1, 1);

	MulAcc(0, 2);

	R[2] = c;

	R[3] = d + A[0] * B[3] + A[1] * B[2] + A[2] * B[1] + A[3] * B[0];

}

static void CombaMultiplyBottom8(word *R, const word *A, const word *B)

{

	dword p;

	word c, d, e;

	p = (dword)A[0] * B[0];

	R[0] = LOW_WORD(p);

	c = HIGH_WORD(p);

	d = e = 0;

	MulAcc(0, 1);

	MulAcc(1, 0);

	SaveMulAcc(1, 2, 0);

	MulAcc(1, 1);

	MulAcc(0, 2);

	SaveMulAcc(2, 0, 3);

	MulAcc(1, 2);

	MulAcc(2, 1);

	MulAcc(3, 0);

	SaveMulAcc(3, 0, 4);

	MulAcc(1, 3);

	MulAcc(2, 2);

	MulAcc(3, 1);

	MulAcc(4, 0);

	SaveMulAcc(4, 0, 5);

	MulAcc(1, 4);

	MulAcc(2, 3);

	MulAcc(3, 2);

	MulAcc(4, 1);

	MulAcc(5, 0);

	SaveMulAcc(5, 0, 6);

	MulAcc(1, 5);

	MulAcc(2, 4);

	MulAcc(3, 3);

	MulAcc(4, 2);

	MulAcc(5, 1);

	MulAcc(6, 0);

	R[6] = c;

	R[7] = d + A[0] * B[7] + A[1] * B[6] + A[2] * B[5] + A[3] * B[4] +

				A[4] * B[3] + A[5] * B[2] + A[6] * B[1] + A[7] * B[0];

}

#undef MulAcc

#undef SaveMulAcc

static void AtomicInverseModPower2(word *C, word A0, word A1)

{

	assert(A0%2==1);

	dword A=MAKE_DWORD(A0, A1), R=A0%8;

	for (unsigned i=3; i<2*WORD_BITS; i*=2)

		R = R*(2-R*A);

	assert(R*A==1);

	C[0] = LOW_WORD(R);

	C[1] = HIGH_WORD(R);

}

// ********************************************************

#define A0		A

#define A1		(A+N2)

#define B0		B

#define B1		(B+N2)

#define T0		T

#define T1		(T+N2)

#define T2		(T+N)

#define T3		(T+N+N2)

#define R0		R

#define R1		(R+N2)

#define R2		(R+N)

#define R3		(R+N+N2)

// R[2*N] - result = A*B

// T[2*N] - temporary work space

// A[N] --- multiplier

// B[N] --- multiplicant

void RecursiveMultiply(word *R, word *T, const word *A, const word *B, unsigned int N)

{

	assert(N>=2 && N%2==0);

	if (N==2)

		AtomicMultiply(R, A, B);

	else if (N==4)

		CombaMultiply4(R, A, B);

	else if (N==8)

		CombaMultiply8(R, A, B);

	else

	{

		const unsigned int N2 = N/2;

		int carry;

		int aComp = Compare(A0, A1, N2);

		int bComp = Compare(B0, B1, N2);

		switch (2*aComp + aComp + bComp)

		{

		case -4:

			Subtract(R0, A1, A0, N2);

			Subtract(R1, B0, B1, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			Subtract(T1, T1, R0, N2);

			carry = -1;

			break;

		case -2:

			Subtract(R0, A1, A0, N2);

			Subtract(R1, B0, B1, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			carry = 0;

			break;

		case 2:

			Subtract(R0, A0, A1, N2);

			Subtract(R1, B1, B0, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			carry = 0;

			break;

		case 4:

			Subtract(R0, A1, A0, N2);

			Subtract(R1, B0, B1, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			Subtract(T1, T1, R1, N2);

			carry = -1;

			break;

		default:

			SetWords(T0, 0, N);

			carry = 0;

		}

		RecursiveMultiply(R0, T2, A0, B0, N2);

		RecursiveMultiply(R2, T2, A1, B1, N2);

		// now T[01] holds (A1-A0)*(B0-B1), R[01] holds A0*B0, R[23] holds A1*B1

		carry += Add(T0, T0, R0, N);

		carry += Add(T0, T0, R2, N);

		carry += Add(R1, R1, T0, N);

		assert (carry >= 0 && carry <= 2);

		Increment(R3, N2, carry);

	}

}

// R[2*N] - result = A*A

// T[2*N] - temporary work space

// A[N] --- number to be squared

void RecursiveSquare(word *R, word *T, const word *A, unsigned int N)

{

	assert(N && N%2==0);

	if (N==2)

		AtomicMultiply(R, A, A);

	else if (N==4)

	{

		// VC60 workaround: MSVC 6.0 has an optimization problem that makes

		// (dword)A*B where either A or B has been cast to a dword before

		// very expensive. Revisit a CombaSquare4() function when this

		// problem is fixed.               

// WARNING: KeithR. 05/08/03 This routine doesn't work with gcc on hardware

// either.  I've completely removed it.  It may be worth looking into sometime

// in the future.  Therefore, we use the CombaMultiply4 on all targets.

//#ifdef __WINS__

		CombaMultiply4(R, A, A);

/*#else

		CombaSquare4(R, A);

#endif*/

	}

	else

	{

		const unsigned int N2 = N/2;

		RecursiveSquare(R0, T2, A0, N2);

		RecursiveSquare(R2, T2, A1, N2);

		RecursiveMultiply(T0, T2, A0, A1, N2);

		word carry = Add(R1, R1, T0, N);

		carry += Add(R1, R1, T0, N);

		Increment(R3, N2, carry);

	}

}

// R[N] - bottom half of A*B

// T[N] - temporary work space

// A[N] - multiplier

// B[N] - multiplicant

void RecursiveMultiplyBottom(word *R, word *T, const word *A, const word *B, unsigned int N)

{

	assert(N>=2 && N%2==0);

	if (N==2)

		AtomicMultiplyBottom(R, A, B);

	else if (N==4)

		CombaMultiplyBottom4(R, A, B);

	else if (N==8)

		CombaMultiplyBottom8(R, A, B);

	else

	{

		const unsigned int N2 = N/2;

		RecursiveMultiply(R, T, A0, B0, N2);

		RecursiveMultiplyBottom(T0, T1, A1, B0, N2);

		Add(R1, R1, T0, N2);

		RecursiveMultiplyBottom(T0, T1, A0, B1, N2);

		Add(R1, R1, T0, N2);

	}

}

// R[N] --- upper half of A*B

// T[2*N] - temporary work space

// L[N] --- lower half of A*B

// A[N] --- multiplier

// B[N] --- multiplicant

void RecursiveMultiplyTop(word *R, word *T, const word *L, const word *A, const word *B, unsigned int N)

{

	assert(N>=2 && N%2==0);

	if (N==2)

	{

		AtomicMultiply(T, A, B);

		((dword *)R)[0] = ((dword *)T)[1];

	}

	else if (N==4)

	{

		CombaMultiply4(T, A, B);

		((dword *)R)[0] = ((dword *)T)[2];

		((dword *)R)[1] = ((dword *)T)[3];

	}

	else

	{

		const unsigned int N2 = N/2;

		int carry;

		int aComp = Compare(A0, A1, N2);

		int bComp = Compare(B0, B1, N2);

		switch (2*aComp + aComp + bComp)

		{

		case -4:

			Subtract(R0, A1, A0, N2);

			Subtract(R1, B0, B1, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			Subtract(T1, T1, R0, N2);

			carry = -1;

			break;

		case -2:

			Subtract(R0, A1, A0, N2);

			Subtract(R1, B0, B1, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			carry = 0;

			break;

		case 2:

			Subtract(R0, A0, A1, N2);

			Subtract(R1, B1, B0, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			carry = 0;

			break;

		case 4:

			Subtract(R0, A1, A0, N2);

			Subtract(R1, B0, B1, N2);

			RecursiveMultiply(T0, T2, R0, R1, N2);

			Subtract(T1, T1, R1, N2);

			carry = -1;

			break;

		default:

			SetWords(T0, 0, N);

			carry = 0;

		}

		RecursiveMultiply(T2, R0, A1, B1, N2);

		// now T[01] holds (A1-A0)*(B0-B1), T[23] holds A1*B1

		CopyWords(R0, L+N2, N2);

		word c2 = Subtract(R0, R0, L, N2);

		c2 += Subtract(R0, R0, T0, N2);

		word t = (Compare(R0, T2, N2) == -1);

		carry += t;

		carry += Increment(R0, N2, c2+t);

		carry += Add(R0, R0, T1, N2);

		carry += Add(R0, R0, T3, N2);

		CopyWords(R1, T3, N2);

		assert (carry >= 0 && carry <= 2);

		Increment(R1, N2, carry);

	}

}

// R[NA+NB] - result = A*B

// T[NA+NB] - temporary work space

// A[NA] ---- multiplier

// B[NB] ---- multiplicant

void AsymmetricMultiply(word *R, word *T, const word *A, unsigned int NA, const word *B, unsigned int NB)

{

	if (NA == NB)

	{

		if (A == B)

			RecursiveSquare(R, T, A, NA);

		else

			RecursiveMultiply(R, T, A, B, NA);

		return;

	}

	if (NA > NB)

	{

		TClassSwap(A, B);

		TClassSwap(NA, NB);

		//std::swap(A, B);

		//std::swap(NA, NB);

	}

	assert(NB % NA == 0);

	assert((NB/NA)%2 == 0); 	// NB is an even multiple of NA

	if (NA==2 && !A[1])

	{

		switch (A[0])

		{

		case 0:

			SetWords(R, 0, NB+2);

			return;

		case 1:

			CopyWords(R, B, NB);

			R[NB] = R[NB+1] = 0;

			return;

		default:

			R[NB] = LinearMultiply(R, B, A[0], NB);

			R[NB+1] = 0;

			return;

		}

	}

	RecursiveMultiply(R, T, A, B, NA);

	CopyWords(T+2*NA, R+NA, NA);

	unsigned i;

	for (i=2*NA; i<NB; i+=2*NA)

		RecursiveMultiply(T+NA+i, T, A, B+i, NA);

	for (i=NA; i<NB; i+=2*NA)

		RecursiveMultiply(R+i, T, A, B+i, NA);

	if (Add(R+NA, R+NA, T+2*NA, NB-NA))

		Increment(R+NB, NA);

}

// R[N] ----- result = A inverse mod 2**(WORD_BITS*N)

// T[3*N/2] - temporary work space

// A[N] ----- an odd number as input

void RecursiveInverseModPower2(word *R, word *T, const word *A, unsigned int N)

{

	if (N==2)

		AtomicInverseModPower2(R, A[0], A[1]);

	else

	{

		const unsigned int N2 = N/2;

		RecursiveInverseModPower2(R0, T0, A0, N2);

		T0[0] = 1;

		SetWords(T0+1, 0, N2-1);

		RecursiveMultiplyTop(R1, T1, T0, R0, A0, N2);

		RecursiveMultiplyBottom(T0, T1, R0, A1, N2);

		Add(T0, R1, T0, N2);

		TwosComplement(T0, N2);

		RecursiveMultiplyBottom(R1, T1, R0, T0, N2);

	}

}

#undef A0

#undef A1

#undef B0

#undef B1

#undef T0

#undef T1

#undef T2

#undef T3

#undef R0

#undef R1

#undef R2

#undef R3

// R[N] --- result = X/(2**(WORD_BITS*N)) mod M

// T[3*N] - temporary work space

// X[2*N] - number to be reduced

// M[N] --- modulus

// U[N] --- multiplicative inverse of M mod 2**(WORD_BITS*N)

void MontgomeryReduce(word *R, word *T, const word *X, const word *M, const word *U, unsigned int N)

{

	RecursiveMultiplyBottom(R, T, X, U, N);

	RecursiveMultiplyTop(T, T+N, X, R, M, N);

	if (Subtract(R, X+N, T, N))

	{

#ifdef _DEBUG

		word carry = Add(R, R, M, N);

		assert(carry);

#else

		Add(R, R, M, N);

#endif

	}

}

// do a 3 word by 2 word divide, returns quotient and leaves remainder in A

static word SubatomicDivide(word *A, word B0, word B1)

{

	// assert {A[2],A[1]} < {B1,B0}, so quotient can fit in a word

	assert(A[2] < B1 || (A[2]==B1 && A[1] < B0));

	dword p, u;

	word Q;

	// estimate the quotient: do a 2 word by 1 word divide

	if (B1+1 == 0)

		Q = A[2];

	else

		Q = word(MAKE_DWORD(A[1], A[2]) / (B1+1));

	// now subtract Q*B from A

	p = (dword) B0*Q;

	u = (dword) A[0] - LOW_WORD(p);

	A[0] = LOW_WORD(u);

	u = (dword) A[1] - HIGH_WORD(p) - (word)(0-HIGH_WORD(u)) - (dword)B1*Q;

	A[1] = LOW_WORD(u);

	A[2] += HIGH_WORD(u);

	// Q <= actual quotient, so fix it

	while (A[2] || A[1] > B1 || (A[1]==B1 && A[0]>=B0))

	{

		u = (dword) A[0] - B0;

		A[0] = LOW_WORD(u);

		u = (dword) A[1] - B1 - (word)(0-HIGH_WORD(u));

		A[1] = LOW_WORD(u);

		A[2] += HIGH_WORD(u);

		Q++;

		assert(Q);	// shouldn't overflow

	}

	return Q;

}

// do a 4 word by 2 word divide, returns 2 word quotient in Q0 and Q1

static inline void AtomicDivide(word *Q, const word *A, const word *B)

{

	if (!B[0] && !B[1]) // if divisor is 0, we assume divisor==2**(2*WORD_BITS)

	{

		Q[0] = A[2];

		Q[1] = A[3];

	}

	else

	{

		word T[4];

		T[0] = A[0]; T[1] = A[1]; T[2] = A[2]; T[3] = A[3];

		Q[1] = SubatomicDivide(T+1, B[0], B[1]);

		Q[0] = SubatomicDivide(T, B[0], B[1]);

#ifdef _DEBUG

		// multiply quotient and divisor and add remainder, make sure it equals dividend

		assert(!T[2] && !T[3] && (T[1] < B[1] || (T[1]==B[1] && T[0]<B[0])));

		word P[4];

		AtomicMultiply(P, Q, B);

		Add(P, P, T, 4);

		assert(Mem::Compare((TUint8*)P, 4*WORD_SIZE, (TUint8*)A, 4*WORD_SIZE)==0);

#endif

	}

}

// for use by Divide(), corrects the underestimated quotient {Q1,Q0}

static void CorrectQuotientEstimate(word *R, word *T, word *Q, const word *B, unsigned int N)

{

	assert(N && N%2==0);

	if (Q[1])

	{

		T[N] = T[N+1] = 0;

		unsigned i;

		for (i=0; i<N; i+=4)

			AtomicMultiply(T+i, Q, B+i);

		for (i=2; i<N; i+=4)

			if (AtomicMultiplyAdd(T+i, Q, B+i))

				T[i+5] += (++T[i+4]==0);

	}

	else

	{

		T[N] = LinearMultiply(T, B, Q[0], N);

		T[N+1] = 0;

	}

#ifdef _DEBUG

	word borrow = Subtract(R, R, T, N+2);

	assert(!borrow && !R[N+1]);

#else

	Subtract(R, R, T, N+2);

#endif

	while (R[N] || Compare(R, B, N) >= 0)

	{

		R[N] -= Subtract(R, R, B, N);

		Q[1] += (++Q[0]==0);

		assert(Q[0] || Q[1]); // no overflow

	}

}

// R[NB] -------- remainder = A%B

// Q[NA-NB+2] --- quotient	= A/B

// T[NA+2*NB+4] - temp work space

// A[NA] -------- dividend

// B[NB] -------- divisor

void Divide(word *R, word *Q, word *T, const word *A, unsigned int NA, const word *B, unsigned int NB)

{

	assert(NA && NB && NA%2==0 && NB%2==0);

	assert(B[NB-1] || B[NB-2]);

	assert(NB <= NA);