// integer.cpp - written and placed in the public domain by Wei Dai
// contains public domain code contributed by Alister Lee and Leonard Janke

#include "pch.h"

#ifndef CRYPTOPP_IMPORTS

#include "integer.h"
#include "modarith.h"
#include "nbtheory.h"
#include "asn.h"
#include "oids.h"
#include "words.h"
#include "algparam.h"
#include "pubkey.h"		// for P1363_KDF2
#include "sha.h"

#include <iostream>

#ifdef SSE2_INTRINSICS_AVAILABLE
	#ifdef __GNUC__
		#include <xmmintrin.h>
		#include <signal.h>
		#include <setjmp.h>
		#ifdef CRYPTOPP_MEMALIGN_AVAILABLE
			#include <malloc.h>
		#else
			#include <stdlib.h>
		#endif
	#else
		#include <emmintrin.h>
	#endif
#elif defined(_MSC_VER) && defined(_M_IX86)
	#pragma message("You do not seem to have the Visual C++ Processor Pack installed, so use of SSE2 intrinsics will be disabled.")
#elif defined(__GNUC__) && defined(__i386__)
	#warning "You do not have GCC 3.3 or later, or did not specify -msse2 compiler option, so use of SSE2 intrinsics will be disabled."
#endif

NAMESPACE_BEGIN(CryptoPP)

bool FunctionAssignIntToInteger(const std::type_info &valueType, void *pInteger, const void *pInt)
{
	if (valueType != typeid(Integer))
		return false;
	*reinterpret_cast<Integer *>(pInteger) = *reinterpret_cast<const int *>(pInt);
	return true;
}

static const char s_RunAtStartup = (AssignIntToInteger = FunctionAssignIntToInteger, 0);

#ifdef SSE2_INTRINSICS_AVAILABLE
template <class T>
CPP_TYPENAME AllocatorBase<T>::pointer AlignedAllocator<T>::allocate(size_type n, const void *)
{
	CheckSize(n);
	if (n == 0)
		return NULL;
	if (n >= 4)
	{
		void *p;
	#ifdef CRYPTOPP_MM_MALLOC_AVAILABLE
		while (!(p = _mm_malloc(sizeof(T)*n, 16)))
	#elif defined(CRYPTOPP_MEMALIGN_AVAILABLE)
		while (!(p = memalign(16, sizeof(T)*n)))
	#elif defined(CRYPTOPP_MALLOC_ALIGNMENT_IS_16)
		while (!(p = malloc(sizeof(T)*n)))
	#else
		while (!(p = (byte *)malloc(sizeof(T)*n + 8)))	// assume malloc alignment is at least 8
	#endif
			CallNewHandler();

	#ifdef CRYPTOPP_NO_ALIGNED_ALLOC
		assert(m_pBlock == NULL);
		m_pBlock = p;
		if (!IsAlignedOn(p, 16))
		{
			assert(IsAlignedOn(p, 8));
			p = (byte *)p + 8;
		}
	#endif

		assert(IsAlignedOn(p, 16));
		return (T*)p;
	}
	return new T[n];
}

template <class T>
void AlignedAllocator<T>::deallocate(void *p, size_type n)
{
	memset(p, 0, n*sizeof(T));
	if (n >= 4)
	{
		#ifdef CRYPTOPP_MM_MALLOC_AVAILABLE
			_mm_free(p);
		#elif defined(CRYPTOPP_NO_ALIGNED_ALLOC)
			assert(m_pBlock == p || (byte *)m_pBlock+8 == p);
			free(m_pBlock);
			m_pBlock = NULL;
		#else
			free(p);
		#endif
	}
	else
		delete [] (T *)p;
}
#endif

static 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;
}

static word Increment(word *A, unsigned int N, word B=1)
{
	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;
}

static word Decrement(word *A, unsigned int N, word B=1)
{
	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;
}

static void TwosComplement(word *A, unsigned int N)
{
	Decrement(A, N);
	for (unsigned i=0; i<N; i++)
		A[i] = ~A[i];
}

static word AtomicInverseModPower2(word A)
{
	assert(A%2==1);

	word R=A%8;

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

	assert(R*A==1);
	return R;
}

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

class DWord
{
public:
	DWord() {}

#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
	explicit DWord(word low)
	{
		m_whole = low;
	}
#else
	explicit DWord(word low)
	{
		m_halfs.low = low;
		m_halfs.high = 0;
	}
#endif

	DWord(word low, word high)
	{
		m_halfs.low = low;
		m_halfs.high = high;
	}

	static DWord Multiply(word a, word b)
	{
		DWord r;
		#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
			r.m_whole = (dword)a * b;
		#elif defined(__alpha__)
			r.m_halfs.low = a*b; __asm__("umulh %1,%2,%0" : "=r" (r.m_halfs.high) : "r" (a), "r" (b));
		#elif defined(__ia64__)
			r.m_halfs.low = a*b; __asm__("xmpy.hu %0=%1,%2" : "=f" (r.m_halfs.high) : "f" (a), "f" (b));
		#elif defined(_ARCH_PPC64)
			r.m_halfs.low = a*b; __asm__("mulhdu %0,%1,%2" : "=r" (r.m_halfs.high) : "r" (a), "r" (b) : "cc");
		#elif defined(__x86_64__)
			__asm__("mulq %3" : "=d" (r.m_halfs.high), "=a" (r.m_halfs.low) : "a" (a), "rm" (b) : "cc");
		#elif defined(__mips64)
			__asm__("dmultu %2,%3" : "=h" (r.m_halfs.high), "=l" (r.m_halfs.low) : "r" (a), "r" (b));
		#elif defined(_M_IX86)
			// for testing
			word64 t = (word64)a * b;
			r.m_halfs.high = ((word32 *)(&t))[1];
			r.m_halfs.low = (word32)t;
		#else
			#error can not implement DWord
		#endif
		return r;
	}

	static DWord MultiplyAndAdd(word a, word b, word c)
	{
		DWord r = Multiply(a, b);
		return r += c;
	}

	DWord & operator+=(word a)
	{
		#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
			m_whole = m_whole + a;
		#else
			m_halfs.low += a;
			m_halfs.high += (m_halfs.low < a);
		#endif
		return *this;
	}

	DWord operator+(word a)
	{
		DWord r;
		#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
			r.m_whole = m_whole + a;
		#else
			r.m_halfs.low = m_halfs.low + a;
			r.m_halfs.high = m_halfs.high + (r.m_halfs.low < a);
		#endif
		return r;
	}

	DWord operator-(DWord a)
	{
		DWord r;
		#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
			r.m_whole = m_whole - a.m_whole;
		#else
			r.m_halfs.low = m_halfs.low - a.m_halfs.low;
			r.m_halfs.high = m_halfs.high - a.m_halfs.high - (r.m_halfs.low > m_halfs.low);
		#endif
		return r;
	}

	DWord operator-(word a)
	{
		DWord r;
		#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
			r.m_whole = m_whole - a;
		#else
			r.m_halfs.low = m_halfs.low - a;
			r.m_halfs.high = m_halfs.high - (r.m_halfs.low > m_halfs.low);
		#endif
		return r;
	}

	// returns quotient, which must fit in a word
	word operator/(word divisor);

	word operator%(word a);

	bool operator!() const
	{
	#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
		return !m_whole;
	#else
		return !m_halfs.high && !m_halfs.low;
	#endif
	}

	word GetLowHalf() const {return m_halfs.low;}
	word GetHighHalf() const {return m_halfs.high;}
	word GetHighHalfAsBorrow() const {return 0-m_halfs.high;}

private:
	union
	{
	#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
		dword m_whole;
	#endif
		struct
		{
		#ifdef IS_LITTLE_ENDIAN
			word low;
			word high;
		#else
			word high;
			word low;
		#endif
		} m_halfs;
	};
};

class Word
{
public:
	Word() {}

	Word(word value)
	{
		m_whole = value;
	}

	Word(hword low, hword high)
	{
		m_whole = low | (word(high) << (WORD_BITS/2));
	}

	static Word Multiply(hword a, hword b)
	{
		Word r;
		r.m_whole = (word)a * b;
		return r;
	}

	Word operator-(Word a)
	{
		Word r;
		r.m_whole = m_whole - a.m_whole;
		return r;
	}

	Word operator-(hword a)
	{
		Word r;
		r.m_whole = m_whole - a;
		return r;
	}

	// returns quotient, which must fit in a word
	hword operator/(hword divisor)
	{
		return hword(m_whole / divisor);
	}

	bool operator!() const
	{
		return !m_whole;
	}

	word GetWhole() const {return m_whole;}
	hword GetLowHalf() const {return hword(m_whole);}
	hword GetHighHalf() const {return hword(m_whole>>(WORD_BITS/2));}
	hword GetHighHalfAsBorrow() const {return 0-hword(m_whole>>(WORD_BITS/2));}
	
private:
	word m_whole;
};

// do a 3 word by 2 word divide, returns quotient and leaves remainder in A
template <class S, class D>
S DivideThreeWordsByTwo(S *A, S B0, S B1, D *dummy=NULL)
{
	// assert {A[2],A[1]} < {B1,B0}, so quotient can fit in a S
	assert(A[2] < B1 || (A[2]==B1 && A[1] < B0));

	// estimate the quotient: do a 2 S by 1 S divide
	S Q;
	if (S(B1+1) == 0)
		Q = A[2];
	else
		Q = D(A[1], A[2]) / S(B1+1);

	// now subtract Q*B from A
	D p = D::Multiply(B0, Q);
	D u = (D) A[0] - p.GetLowHalf();
	A[0] = u.GetLowHalf();
	u = (D) A[1] - p.GetHighHalf() - u.GetHighHalfAsBorrow() - D::Multiply(B1, Q);
	A[1] = u.GetLowHalf();
	A[2] += u.GetHighHalf();

	// Q <= actual quotient, so fix it
	while (A[2] || A[1] > B1 || (A[1]==B1 && A[0]>=B0))
	{
		u = (D) A[0] - B0;
		A[0] = u.GetLowHalf();
		u = (D) A[1] - B1 - u.GetHighHalfAsBorrow();
		A[1] = u.GetLowHalf();
		A[2] += u.GetHighHalf();
		Q++;
		assert(Q);	// shouldn't overflow
	}

	return Q;
}

// do a 4 word by 2 word divide, returns 2 word quotient in Q0 and Q1
template <class S, class D>
inline D DivideFourWordsByTwo(S *T, const D &Al, const D &Ah, const D &B)
{
	if (!B) // if divisor is 0, we assume divisor==2**(2*WORD_BITS)
		return D(Ah.GetLowHalf(), Ah.GetHighHalf());
	else
	{
		S Q[2];
		T[0] = Al.GetLowHalf();
		T[1] = Al.GetHighHalf(); 
		T[2] = Ah.GetLowHalf();
		T[3] = Ah.GetHighHalf();
		Q[1] = DivideThreeWordsByTwo<S, D>(T+1, B.GetLowHalf(), B.GetHighHalf());
		Q[0] = DivideThreeWordsByTwo<S, D>(T, B.GetLowHalf(), B.GetHighHalf());
		return D(Q[0], Q[1]);
	}
}

// returns quotient, which must fit in a word
inline word DWord::operator/(word a)
{
	#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
		return word(m_whole / a);
	#else
		hword r[4];
		return DivideFourWordsByTwo<hword, Word>(r, m_halfs.low, m_halfs.high, a).GetWhole();
	#endif
}

inline word DWord::operator%(word a)
{
	#ifdef CRYPTOPP_NATIVE_DWORD_AVAILABLE
		return word(m_whole % a);
	#else
		if (a < (word(1) << (WORD_BITS/2)))
		{
			hword h = hword(a);
			word r = m_halfs.high % h;
			r = ((m_halfs.low >> (WORD_BITS/2)) + (r << (WORD_BITS/2))) % h;
			return hword((hword(m_halfs.low) + (r << (WORD_BITS/2))) % h);
		}
		else
		{
			hword r[4];
			DivideFourWordsByTwo<hword, Word>(r, m_halfs.low, m_halfs.high, a);
			return Word(r[0], r[1]).GetWhole();
		}
	#endif
}

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

class Portable
{
public:
	static word Add(word *C, const word *A, const word *B, unsigned int N);
	static word Subtract(word *C, const word *A, const word *B, unsigned int N);

	static inline void Multiply2(word *C, const word *A, const word *B);
	static inline word Multiply2Add(word *C, const word *A, const word *B);
	static void Multiply4(word *C, const word *A, const word *B);
	static void Multiply8(word *C, const word *A, const word *B);
	static inline unsigned int MultiplyRecursionLimit() {return 8;}

	static inline void Multiply2Bottom(word *C, const word *A, const word *B);
	static void Multiply4Bottom(word *C, const word *A, const word *B);
	static void Multiply8Bottom(word *C, const word *A, const word *B);
	static inline unsigned int MultiplyBottomRecursionLimit() {return 8;}

	static void Square2(word *R, const word *A);
	static void Square4(word *R, const word *A);
	static void Square8(word *R, const word *A) {assert(false);}
	static inline unsigned int SquareRecursionLimit() {return 4;}
};

word Portable::Add(word *C, const word *A, const word *B, unsigned int N)
{
	assert (N%2 == 0);

	DWord u(0, 0);
	for (unsigned int i = 0; i < N; i+=2)
	{
		u = DWord(A[i]) + B[i] + u.GetHighHalf();
		C[i] = u.GetLowHalf();
		u = DWord(A[i+1]) + B[i+1] + u.GetHighHalf();
		C[i+1] = u.GetLowHalf();
	}
	return u.GetHighHalf();
}

word Portable::Subtract(word *C, const word *A, const word *B, unsigned int N)
{
	assert (N%2 == 0);

	DWord u(0, 0);
	for (unsigned int i = 0; i < N; i+=2)
	{
		u = (DWord) A[i] - B[i] - u.GetHighHalfAsBorrow();
		C[i] = u.GetLowHalf();
		u = (DWord) A[i+1] - B[i+1] - u.GetHighHalfAsBorrow();
		C[i+1] = u.GetLowHalf();
	}
	return 0-u.GetHighHalf();
}

void Portable::Multiply2(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::Multiply(D[di], D[di+2]);
	D[1] = D[3] = 0;
	unsigned int si = ai + !bi;
	word s = D[si];

	DWord A0B0 = DWord::Multiply(A[0], B[0]);
	C[0] = A0B0.GetLowHalf();

	DWord A1B1 = DWord::Multiply(A[1], B[1]);
	DWord t = (DWord) A0B0.GetHighHalf() + A0B0.GetLowHalf() + d.GetLowHalf() + A1B1.GetLowHalf();
	C[1] = t.GetLowHalf();

	t = A1B1 + t.GetHighHalf() + A0B0.GetHighHalf() + d.GetHighHalf() + A1B1.GetHighHalf() - s;
	C[2] = t.GetLowHalf();
	C[3] = t.GetHighHalf();
}

inline void Portable::Multiply2Bottom(word *C, const word *A, const word *B)
{
	DWord t = DWord::Multiply(A[0], B[0]);
	C[0] = t.GetLowHalf();
	C[1] = t.GetHighHalf() + A[0]*B[1] + A[1]*B[0];
}

word Portable::Multiply2Add(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;