
///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//                           Double Integer (uint128)                        //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#define HWSHIFT			// uncomment this flag to ensure hardware compatibility of bit shift
						// (number of shifts is masked to total number of bits)

/*
typedef char				s8;
typedef short				s16;
typedef long				s32;
typedef long long			s64;

typedef unsigned char		u8;
typedef unsigned short		u16;
typedef unsigned long		u32;
typedef unsigned long long	u64;
*/

#ifndef DBLUINT

#define DBLUINT uint128		// double integer class unsigned
#define BASEUINT u64		// base integer unsigned
#define HALFUINT u32		// half integer unsigned

#define DBLSINT sint128		// double integer class signed
#define BASESINT s64		// base integer signed
#define HALFSINT s32		// half integer signed

#define DBLFIX fix128		// double fixed point class

#define USEASM				// can use assembler

#endif // DBLUINT


#define BASEBITS ((int)sizeof(BASEUINT)*8)	// number of bits of base integer
#define BASELAST ((BASEUINT)1 << (BASEBITS-1)) // highest bit of base integer
#define HALFBITS ((int)sizeof(HALFUINT)*8)	// number of bits of half integer
#define HALFMASK (((BASEUINT)1 << HALFBITS)-1) // mask of half integer
#define BASENOT  (~(BASEUINT)0)			// all bits '1' (or ((BASEUINT)-1) )

#define DBLFIX_MAX ((double)(BASELAST-1) + (double)BASENOT/BASELAST/2) // max. decimal number as double
#define DBLFIX_MIN (-(double)BASELAST) // min. decimal number as double
#define DBLFIX_DEC ((int)(BASEBITS*0.301)) // default number of decimal places

#ifndef WIN64
#undef USEASM				// can use assembler
#endif // WIN64

// Approx. timmings of short functions:
//  AMD 4 GHz: uint128 (ASM) 1.5 ns, uint128ref (C++) 2.0 ns
//  Intel 3 GHz: uint128 (ASM): 2.5 ns, uint128ref (C++) 3.0 ns
//  Intel 2 GHz: uint128 (ASM): 3.0 ns, uint128ref (C++) 4.0 ns
#ifdef USEASM
#define SOFTASM				// uncomment this flag to use short assembler codes
							// tip: in some cases short C++ inline codes can be
							//      more effective than assembler functions
#endif // USEASM

// constants

// factorial table
DBLUINT* DBLUINT::Fact_;
DBLSINT* DBLSINT::Fact_;
DBLFIX* DBLFIX::Fact_;
DBLFIX* DBLFIX::FactRec_;

// initialize factorial tablbe
bool DBLUINT::FactOK = DBLUINT::FactInit();
bool DBLSINT::FactOK = DBLSINT::FactInit();
bool DBLFIX::FactOK = DBLFIX::FactInit();

// max. factorial
BASEUINT DBLUINT::FactMax = 0;
BASEUINT DBLSINT::FactMax = 0;
BASEUINT DBLFIX::FactMax = 0;

// Euler's constant
DBLFIX DBLFIX::Eul;
CText DBLFIX::EulTxt;

// Euler's constant ^2
DBLFIX DBLFIX::Eul2;
CText DBLFIX::Eul2Txt;

// Euler's constant ^4
DBLFIX DBLFIX::Eul4;
CText DBLFIX::Eul4Txt;

// Euler's constant ^8
DBLFIX DBLFIX::Eul8;
CText DBLFIX::Eul8Txt;

// 1/Euler's constant
DBLFIX DBLFIX::Eul1;
CText DBLFIX::Eul1Txt;

// 1/Euler's constant ^2
DBLFIX DBLFIX::Eul12;
CText DBLFIX::Eul12Txt;

// 1/Euler's constant ^4
DBLFIX DBLFIX::Eul14;
CText DBLFIX::Eul14Txt;

// 1/Euler's constant ^8
DBLFIX DBLFIX::Eul18;
CText DBLFIX::Eul18Txt;

// Pi constant
DBLFIX DBLFIX::Pi;
CText DBLFIX::PiTxt;

// Pi*2 constant
DBLFIX DBLFIX::Pi2;
CText DBLFIX::Pi2Txt;

// Pi/2 constant
DBLFIX DBLFIX::Pi12;
CText DBLFIX::Pi12Txt;

// Pi/4 constant
DBLFIX DBLFIX::Pi14;
CText DBLFIX::Pi14Txt;

// Pi/180 constant
DBLFIX DBLFIX::Pi180;
CText DBLFIX::Pi180Txt;

// 180/Pi constant
DBLFIX DBLFIX::Pi180R;
CText DBLFIX::Pi180RTxt;

// ln(10) constant
DBLFIX DBLFIX::Ln10;
CText DBLFIX::Ln10Txt;

// 1/ln(10) constant
DBLFIX DBLFIX::Ln110;
CText DBLFIX::Ln110Txt;

// ln(2) constant
DBLFIX DBLFIX::Ln2;
CText DBLFIX::Ln2Txt;

// 16*ln(2) constant
DBLFIX DBLFIX::Ln2_16;
CText DBLFIX::Ln2_16Txt;

// sqrt(2) constant
DBLFIX DBLFIX::Sqrt2;
CText DBLFIX::Sqrt2Txt;

// 1/sqrt(2) constant
DBLFIX DBLFIX::Sqrt12;
CText DBLFIX::Sqrt12Txt;

// initialize static constants
bool DBLFIX::InitOK = DBLFIX::ConstInit();

// assembler
#ifdef USEASM
extern "C" void AddU128A_x64(DBLUINT*, const DBLUINT*);
extern "C" void AddU128B_x64(DBLUINT*, BASEUINT);
extern "C" void AddU128C_x64(DBLUINT*, const DBLUINT*, const DBLUINT*);
extern "C" void AddU128D_x64(DBLUINT*, const DBLUINT*, BASEUINT);

extern "C" void SubU128A_x64(DBLUINT*, const DBLUINT*);
extern "C" void SubU128B_x64(DBLUINT*, BASEUINT);
extern "C" void SubU128C_x64(DBLUINT*, const DBLUINT*, const DBLUINT*);
extern "C" void SubU128D_x64(DBLUINT*, const DBLUINT*, BASEUINT);
extern "C" void SubU128E_x64(DBLUINT*, BASEUINT, const DBLUINT*);
extern "C" void SubU128F_x64(DBLUINT*, const DBLUINT*);
extern "C" void SubU128G_x64(DBLUINT*, BASEUINT);

extern "C" void IncU128A_x64(DBLUINT*);
extern "C" void IncU128B_x64(DBLUINT*, const DBLUINT*);
extern "C" void DecU128A_x64(DBLUINT*);
extern "C" void DecU128B_x64(DBLUINT*, const DBLUINT*);

extern "C" void MulU128A_x64(DBLUINT*, const DBLUINT*, const DBLUINT*);
extern "C" void MulU128B_x64(DBLUINT*, const DBLUINT*, BASEUINT);
extern "C" void MulU128C_x64(DBLUINT*, BASEUINT, BASEUINT);
extern "C" void MulU128D_x64(DBLUINT*, const DBLUINT*, const DBLUINT*, DBLUINT*);
extern "C" int MulU128E_x64(DBLUINT*, const DBLUINT*);

extern "C" void SqrU128A_x64(DBLUINT*, const DBLUINT*);
extern "C" void SqrU128B_x64(DBLUINT*, BASEUINT);

extern "C" void DivU128A_x64(DBLUINT*, const DBLUINT*, const DBLUINT*, DBLUINT*);
extern "C" BASEUINT DivU128B_x64(DBLUINT*, const DBLUINT*, BASEUINT);
extern "C" HALFUINT DivU128C_x64(DBLUINT*, const DBLUINT*, HALFUINT);
extern "C" void DivU128D_x64(DBLUINT*, const DBLUINT*, const DBLUINT*, const DBLUINT*, DBLUINT*);
extern "C" int DivU128E_x64(DBLUINT*, const DBLUINT*);

extern "C" void NegU128A_x64(DBLUINT*);
extern "C" void NegU128B_x64(DBLUINT*, const DBLUINT*);
extern "C" void NegU128C_x64(DBLUINT*, BASEUINT);

extern "C" s64 BitsU128_x64(const DBLUINT*);
#endif // USEASM

///////////////////////////////////////////////////////////////////////////////
// multiply 2 base entries, r1:r0 = a * b1:b0

void DBLUINT::MulBB(BASEUINT* r0, BASEUINT* r1, BASEUINT a, HALFUINT b0, HALFUINT b1)
{
	//       a1a0
	//     x b1b0
	// ----------
	//       a0b0 ...   t0
	//     a0b1   ...  t1
	//     a1b0   ...  t2
	//   a1b1     ... t3

	// multiply elements (max. FFFF * FFFF = FFFE0001)
	HALFUINT a0 = (HALFUINT)(a & HALFMASK);
	BASEUINT t0 = (BASEUINT)a0*b0;
	BASEUINT t1 = (BASEUINT)a0*b1;
	HALFUINT a1 = (HALFUINT)(a >> HALFBITS);
	BASEUINT t2 = (BASEUINT)a1*b0;
	BASEUINT t3 = (BASEUINT)a1*b1;

	// add elements
	t1 += t0 >> HALFBITS; // max. FFFE0001 + FFFE = FFFEFFFF, no carry yet
	t1 += t2; // max. FFFEFFFF + FFFE0001 = 1FFFD0000, it can carry
	if (t1 < t2) t3 += (BASEUINT)1 << HALFBITS; // add carry, FFFE0001 + 10000 = FFFF0001, no carry

	// result, max. FFFFFFFF * FFFFFFFF = FFFFFFFE:00000001
	*r0 = (t1 << HALFBITS) + (t0 & HALFMASK); // result low, FFFF0000 + FFFF = FFFFFFFF, no carry
	*r1 = t3 + (t1 >> HALFBITS); // result high, FFFF0001 + FFFD = FFFFFFFE, no carry 
}

///////////////////////////////////////////////////////////////////////////////
// multiply base entry by half entry, r1:r0 = a * 0:b0

void DBLUINT::MulBH(BASEUINT* r0, BASEUINT* r1, BASEUINT a, HALFUINT b0)
{
	//       a1a0
	//     x   b0
	// ----------
	//       a0b0 ...  t0
	//     a1b0   ... t1

	// multiply elements (max. FFFF * FFFF = FFFE0001)
	HALFUINT a0 = (HALFUINT)(a & HALFMASK);
	BASEUINT t0 = (BASEUINT)a0*b0;
	HALFUINT a1 = (HALFUINT)(a >> HALFBITS);
	BASEUINT t1 = (BASEUINT)a1*b0;

	// add elements
	t1 += t0 >> HALFBITS; // max. FFFE0001 + FFFE = FFFEFFFF, no carry yet

	// result, max. FFFFFFFF * 0000FFFF = 0000FFFE:FFFF0001
	*r0 = (t1 << HALFBITS) + (t0 & HALFMASK); // result low, FFFF0000 + FFFF = FFFFFFFF, no carry
	*r1 = t1 >> HALFBITS; // result high, max. 0000FFFE, no carry
}

///////////////////////////////////////////////////////////////////////////////
// initialize factorial table

bool DBLUINT::FactInit()
{
	// already initialized
	if (DBLUINT::Fact_ != NULL) return true;

	// factorial value 0
	DBLUINT::Fact_ = (DBLUINT*)malloc(sizeof(DBLUINT));
	if (DBLUINT::Fact_ == NULL) return false;
	DBLUINT::Fact_->Set0();

	// loop
	DBLUINT n, n2, n3;
	n.Set1();
	BASEUINT i;
	for (i = 1; ; i++)
	{
		// get new member
		n2 = n;
		n.Mul(i);
		n3.Div(n, i);
		if (n2.NEqu(n3)) break; // overflow

		// resize table
		DBLUINT::Fact_ = (DBLUINT*)realloc(DBLUINT::Fact_, (size_t)((i+1)*sizeof(DBLUINT)));
		if (DBLUINT::Fact_ == NULL)
		{
			DBLUINT::FactMax = 0;
			return false;
		}

		// store new member
		DBLUINT::FactMax = i;
		DBLUINT::Fact_[i].Set(n);
	}

	return true;
}

///////////////////////////////////////////////////////////////////////////////
// factorial this = num!

void DBLUINT::Fact(BASEUINT num)
{
	if (num <= this->FactMax)
		this->Set(this->Fact_[num]);
	else
		this->SetMax();
}

void DBLUINT::Fact()
{
	if (this->N1() == 0)
		this->Fact(this->N0());
	else
		this->SetMax();
}

///////////////////////////////////////////////////////////////////////////////
// exchange two numbers

void DBLUINT::Exchange(DBLUINT& num)
{
	BASEUINT r = num.N0();
	num.N0() = this->N0();
	this->N0() = r;

	r = num.N1();
	num.N1() = this->N1();
	this->N1() = r;
}

///////////////////////////////////////////////////////////////////////////////
// this += num

void DBLUINT::Add(const DBLUINT& num)
{
#ifdef SOFTASM
	AddU128A_x64(this, &num);
#else // SOFTASM
	BASEUINT r = this->N0() + num.N0();
	this->N1() = this->N1() + num.N1() + ((r < this->N0()) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this += 0:num

void DBLUINT::Add(BASEUINT num)
{
#ifdef SOFTASM
	AddU128B_x64(this, num);
#else // SOFTASM
	BASEUINT r = this->N0() + num;
	this->N1() = this->N1() + ((r < num) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 + num2

void DBLUINT::Add(const DBLUINT& num1, const DBLUINT& num2)
{
#ifdef SOFTASM
	AddU128C_x64(this, &num1, &num2);
#else // SOFTASM
	BASEUINT r = num1.N0() + num2.N0();
	this->N1() = num1.N1() + num2.N1() + ((r < num1.N0()) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 + 0:num2

void DBLUINT::Add(const DBLUINT& num1, BASEUINT num2)
{
#ifdef SOFTASM
	AddU128D_x64(this, &num1, num2);
#else // SOFTASM
	BASEUINT r = num1.N0() + num2;
	this->N1() = num1.N1() + ((r < num2) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this -= num

void DBLUINT::Sub(const DBLUINT& num)
{
#ifdef SOFTASM
	SubU128A_x64(this, &num);
#else // SOFTASM
	BASEUINT r = this->N0() - num.N0();
	this->N1() = this->N1() - num.N1() - ((this->N0() < num.N0()) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this -= 0:num

void DBLUINT::Sub(BASEUINT num)
{
#ifdef SOFTASM
	SubU128B_x64(this, num);
#else // SOFTASM
	BASEUINT r = this->N0() - num;
	this->N1() = this->N1() - ((this->N0() < num) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 - num2

void DBLUINT::Sub(const DBLUINT& num1, const DBLUINT& num2)
{
#ifdef SOFTASM
	SubU128C_x64(this, &num1, &num2);
#else // SOFTASM
	BASEUINT r = num1.N0() - num2.N0();
	this->N1() = num1.N1() - num2.N1() - ((num1.N0() < num2.N0()) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 - 0:num2

void DBLUINT::Sub(const DBLUINT& num1, BASEUINT num2)
{
#ifdef SOFTASM
	SubU128D_x64(this, &num1, num2);
#else // SOFTASM
	BASEUINT r = num1.N0() - num2;
	this->N1() = num1.N1() - ((num1.N0() < num2) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num1 - num2

void DBLUINT::Sub(BASEUINT num1, const DBLUINT& num2)
{
#ifdef SOFTASM
	SubU128E_x64(this, num1, &num2);
#else // SOFTASM
	BASEUINT r = num1 - num2.N0();
	this->N1() = 0 - num2.N1() - ((num1 < num2.N0()) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num - this

void DBLUINT::InvSub(const DBLUINT& num)
{
#ifdef SOFTASM
	SubU128F_x64(this, &num);
#else // SOFTASM
	BASEUINT r = num.N0() - this->N0();
	this->N1() = num.N1() - this->N1() - ((num.N0() < this->N0()) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num - this

void DBLUINT::InvSub(BASEUINT num)
{
#ifdef SOFTASM
	SubU128G_x64(this, num);
#else // SOFTASM
	BASEUINT r = num - this->N0();
	this->N1() = 0 - this->N1() - ((num < this->N0()) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this++

void DBLUINT::Inc()
{
#ifdef SOFTASM
	IncU128A_x64(this);
#else // SOFTASM
	BASEUINT r = this->N0() + 1;
	this->N0() = r;
	this->N1() = this->N1() + ((r == 0) ? 1 : 0);
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num + 1

void DBLUINT::Inc(const DBLUINT& num)
{
#ifdef SOFTASM
	IncU128B_x64(this, &num);
#else // SOFTASM
	BASEUINT r = num.N0() + 1;
	this->N0() = r;
	this->N1() = num.N1() + ((r == 0) ? 1 : 0);
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this--

void DBLUINT::Dec()
{
#ifdef SOFTASM
	DecU128A_x64(this);
#else // SOFTASM
	BASEUINT r = this->N0();
	this->N1() = this->N1() - ((r == 0) ? 1 : 0);
	this->N0() = r - 1;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num - 1

void DBLUINT::Dec(const DBLUINT& num)
{
#ifdef SOFTASM
	DecU128B_x64(this, &num);
#else // SOFTASM
	BASEUINT r = num.N0();
	this->N1() = num.N1() - ((r == 0) ? 1 : 0);
	this->N0() = r - 1;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * num2
// AMD 4 GHz: uint128 (ASM) 3 ns, uint128ref (C++) 10 ns
// Intel 3 GHz: uint128 (ASM): 4 ns, uint128ref (C++): 10 ns
// Intel 2 GHz: uint128 (ASM): 5 ns, uint128ref (C++): 11 ns

void DBLUINT::Mul(const DBLUINT& num1, const DBLUINT& num2)
{ 
#ifdef USEASM
	MulU128A_x64(this, &num1, &num2);
#else // USEASM

	// small number
	if (num2.N1() == 0)
	{
		this->Mul(num1, num2.N0());
		return;
	}

	// r = a * b

	//          a3a2a1a0
	//        x b3b2b1b0
	// ---------------------------
	//         |    a0b0    ...     r1r0
	//         |  a1b0      ...   r2r1
	//         |a2b0        ... r3r2
	//       a3|b0          ... r3
	//         |  a0b1      ...   r2r1
	//         |a1b1        ... r3r2
	//       a2|b1          ... r3
	//     a3b1| 
	//         |a0b2        ... r3r2
	//       a1|b2          ... r3
	//     a2b2|
	//   a3b2  |
	//       a0|b3          ... r3
	//     a1b3|
	//   a2b3  |
	// a3b3    |

	// result max. FFFF x FFFF = FFFE0001

	// load operands (value max FFFF)
	const HALFUINT a0 = ((const HALFUINT*)num1.m_N)[0];
	const HALFUINT a1 = ((const HALFUINT*)num1.m_N)[1];
	const HALFUINT a2 = ((const HALFUINT*)num1.m_N)[2];
	const HALFUINT a3 = ((const HALFUINT*)num1.m_N)[3];

	const HALFUINT b0 = ((const HALFUINT*)num2.m_N)[0];
	const HALFUINT b1 = ((const HALFUINT*)num2.m_N)[1];
	const HALFUINT b2 = ((const HALFUINT*)num2.m_N)[2];
	const HALFUINT b3 = ((const HALFUINT*)num2.m_N)[3];

	// rH (no carry required)
	BASEUINT rH = (((BASEUINT)a0*b3 + (BASEUINT)a1*b2 + (BASEUINT)a2*b1 + (BASEUINT)a3*b0) << HALFBITS)
		+ (BASEUINT)a0*b2 + (BASEUINT)a1*b1 + (BASEUINT)a2*b0;

	// rL
	BASEUINT rL = (BASEUINT)a0*b0;

	// middle
	BASEUINT r = (BASEUINT)a1*b0 + (rL >> HALFBITS); // max FFFEFFFF, no carry yet
	BASEUINT r2 = (BASEUINT)a0*b1;
	r += r2;
	rH += BASEUINT((r < r2) ? 1 : 0) << HALFBITS; // add carry

	// save result
	this->N0() = (rL & HALFMASK) | (r << HALFBITS);
	this->N1() = rH + (r >> HALFBITS);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * 0:num2
// AMD 4 GHz: uint128 (ASM) 2.5 ns, uint128ref (C++) 7 ns
// Intel 3 GHz: uint128 (ASM): 3 ns, uint128ref (C++): 8 ns
// Intel 2 GHz: uint128 (ASM): 4 ns, uint128ref (C++): 8 ns

void DBLUINT::Mul(const DBLUINT& num1, BASEUINT num2)
{
#ifdef USEASM
	MulU128B_x64(this, &num1, num2);
#else // USEASM
	// r = a * b

	//          a3a2a1a0
	//        x     b1b0
	// ---------------------------
	//         |    a0b0    ...     r1r0
	//         |  a1b0      ...   r2r1
	//         |a2b0        ... r3r2
	//       a3|b0          ... r3
	//         |  a0b1      ...   r2r1
	//         |a1b1        ... r3r2
	//       a2|b1          ... r3
	//     a3b1| 

	// result max. FFFF x FFFF = FFFE0001

	// small number
	if (num1.N1() == 0)
	{
		this->Mul(num1.N0(), num2);
		return;
	}

	// load operands (value max FFFF)
	const HALFUINT a0 = ((const HALFUINT*)num1.m_N)[0];
	const HALFUINT a1 = ((const HALFUINT*)num1.m_N)[1];
	const HALFUINT a2 = ((const HALFUINT*)num1.m_N)[2];
	const HALFUINT a3 = ((const HALFUINT*)num1.m_N)[3];

	const HALFUINT b0 = ((const HALFUINT*)&num2)[0];
	const HALFUINT b1 = ((const HALFUINT*)&num2)[1];

	// rH (no carry required)
	BASEUINT rH = (((BASEUINT)a2*b1 + (BASEUINT)a3*b0) << HALFBITS)
		+ (BASEUINT)a1*b1 + (BASEUINT)a2*b0;

	// rL
	BASEUINT rL = (BASEUINT)a0*b0;

	// middle
	BASEUINT r = (BASEUINT)a1*b0 + (rL >> HALFBITS); // max FFFEFFFF, no carry yet
	BASEUINT r2 = (BASEUINT)a0*b1;
	r += r2;
	rH += ((r < r2) ? (BASEUINT)1 : 0) << HALFBITS; // add carry

	// save result
	this->N0() = (rL & HALFMASK) | (r << HALFBITS);
	this->N1() = rH + (r >> HALFBITS);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num1 * 0:num2
// AMD 4 GHz: uint128 (ASM) 1.5 ns, uint128ref (C++) 5 ns
// Intel 3 GHz: uint128 (ASM): 2.5 ns, uint128ref (C++): 6 ns
// Intel 2 GHz: uint128 (ASM): 3.5 ns, uint128ref (C++): 7 ns

void DBLUINT::Mul(BASEUINT num1, BASEUINT num2)
{
#ifdef USEASM
	MulU128C_x64(this, num1, num2);
#else // USEASM
	// r = a * b

	//              a1a0
	//        x     b1b0
	// ---------------------------
	//         |    a0b0    ...     r1r0
	//         |  a1b0      ...   r2r1
	//         |  a0b1      ...   r2r1
	//         |a1b1        ... r3r2

	// result max. FFFF x FFFF = FFFE0001

	// load operands (value max FFFF)
	const HALFUINT a0 = ((const HALFUINT*)&num1)[0];
	const HALFUINT a1 = ((const HALFUINT*)&num1)[1];

	const HALFUINT b0 = ((const HALFUINT*)&num2)[0];
	const HALFUINT b1 = ((const HALFUINT*)&num2)[1];

	// rH (no carry required)
	BASEUINT rH = (BASEUINT)a1*b1;

	// rL
	BASEUINT rL = (BASEUINT)a0*b0;

	// middle
	BASEUINT r = (BASEUINT)a1*b0 + (rL >> HALFBITS); // max FFFEFFFF, no carry
	BASEUINT r2 = (BASEUINT)a0*b1;
	r += r2;
	rH += BASEUINT((r < r2) ? 1 : 0) << HALFBITS; // add carry

	// save result
	this->N0() = (rL & HALFMASK) | (r << HALFBITS);
	this->N1() = rH + (r >> HALFBITS);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * num2, with full range of result (high = destination of high part of result)

void DBLUINT::MulHigh(const DBLUINT& num1, const DBLUINT& num2, DBLUINT* high)
{
#ifdef USEASM
	MulU128D_x64(this, &num1, &num2, high);
#else // USEASM

	// prepare operands (r = a * b)
	const BASEUINT* a = &num1.N0();
	const BASEUINT* b = &num2.N0();

	// prepare result accumulator
	BASEUINT r[4];
	BASEUINT t, carry, r0, r1;
	HALFUINT b0, b1;

	// multiply b0 (store to result)
	t = b[0];
	b0 = (HALFUINT)t;
	b1 = (HALFUINT)(t >> HALFBITS);

	if (b1 == 0)
	{
		this->MulBH(&r[0], &carry, a[0], b0);
		this->MulBH(&r0, &r1, a[1], b0);
	}
	else
	{
		this->MulBB(&r[0], &carry, a[0], b0, b1);
		this->MulBB(&r0, &r1, a[1], b0, b1);

	}
	r0 += carry;
	carry = (r0 < carry) ? 1 : 0;
	carry += r1;
	r[1] = r0;
	r[2] = carry;

	// multiply b1 (add to result)
	t = b[1];
	b0 = (HALFUINT)t;
	b1 = (HALFUINT)(t >> HALFBITS);

	if (b1 == 0)
	{
		this->MulBH(&r0, &carry, a[0], b0);
		t = r[1];
		r0 += t;
		carry += (r0 < t) ? 1 : 0;
		r[1] = r0;

		this->MulBH(&r0, &r1, a[1], b0);
	}
	else
	{
		this->MulBB(&r0, &carry, a[0], b0, b1);
		t = r[1];
		r0 += t;
		carry += (r0 < t) ? 1 : 0;
		r[1] = r0;

		this->MulBB(&r0, &r1, a[1], b0, b1);
	}
	r0 += carry;
	carry = (r0 < carry) ? 1 : 0;
	carry += r1;
	t = r[2];
	r0 += t;
	carry += (r0 < t) ? 1 : 0;

	// save result
	this->N0() = r[0];
	this->N1() = r[1];
	high->N0() = r0;
	high->N1() = carry;
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num * 10, returns high carry (0..9)

int DBLUINT::Mul10(const DBLUINT& num)
{
#ifdef USEASM
	return MulU128E_x64(this, &num);
#else // USEASM

	const HALFUINT* a = (const HALFUINT*)num.m_N;

	BASEUINT r0 = a[0] * (BASEUINT)10;
	BASEUINT r1 = a[1] * (BASEUINT)10;
	r1 += (r0 >> HALFBITS);
	this->N0() = (r0 & HALFMASK) | (r1 << HALFBITS);

	BASEUINT r2 = a[2] * (BASEUINT)10;
	r2 += (r1 >> HALFBITS);
	BASEUINT r3 = a[3] * (BASEUINT)10;
	r3 += (r2 >> HALFBITS);
	this->N1() = (r2 & HALFMASK) | (r3 << HALFBITS);

	return (int)(r3 >> HALFBITS);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num * num
// AMD 4 GHz: uint128 (ASM) 2.5 ns, uint128ref (C++) 6 ns
// Intel 3 GHz: uint128 (ASM): 3 ns, uint128ref (C++): 6 ns
// Intel 2 GHz: uint128 (ASM): 4 ns, uint128ref (C++): 7 ns

void DBLUINT::Sqr(const DBLUINT& num)
{
#ifdef USEASM
	SqrU128A_x64(this, &num);
#else // USEASM
	// r = a * a

	//          a3a2a1a0
	//        x a3a2a1a0
	// ---------------------------
	//         |    a0a0    ...     r1r0
	//         |a1a1        ... r3r2
	//     a2a2|
	// a3a3    |

	// 2*
	//         |  a1a0      ...   r2r1
	//         |a2a0        ... r3r2
	//       a3|a0          ... r3
	//       a2|a1          ... r3
	//     a3a1| 
	//   a3a2  |

	// result max. FFFF x FFFF = FFFE0001

	// small number
	if (num.N1() == 0)
	{
		this->Sqr(num.N0());
		return;
	}

	// load operands (value max FFFF)
	const HALFUINT a0 = ((const HALFUINT*)num.m_N)[0];
	const HALFUINT a1 = ((const HALFUINT*)num.m_N)[1];
	const HALFUINT a2 = ((const HALFUINT*)num.m_N)[2];
	const HALFUINT a3 = ((const HALFUINT*)num.m_N)[3];

	// rH (no carry required)
	BASEUINT rH = (((((BASEUINT)a3*a0 + (BASEUINT)a2*a1) << HALFBITS)
		+ (BASEUINT)a2*a0) << 1) + (BASEUINT)a1*a1;

	// rL
	BASEUINT rL = (BASEUINT)a0*a0;

	// middle
	BASEUINT r2 = (BASEUINT)a1*a0;
	BASEUINT r = r2 + (rL >> HALFBITS); // max FFFEFFFF, no carry
	r += r2;
	rH += BASEUINT((r < r2) ? 1 : 0) << HALFBITS; // add carry

	// save result
	this->N0() = (rL & HALFMASK) | (r << HALFBITS);
	this->N1() = rH + (r >> HALFBITS);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num * 0:num
// AMD 4 GHz: uint128 (ASM) 1.5 ns, uint128ref (C++) 4 ns
// Intel 3 GHz: uint128 (ASM): 2.5 ns, uint128ref (C++): 4 ns
// Intel 2 GHz: uint128 (ASM): 3.5 ns, uint128ref (C++): 6 ns

void DBLUINT::Sqr(BASEUINT num)
{
#ifdef USEASM
	SqrU128B_x64(this, num);
#else // USEASM
	// r = a * a

	//          a1a0
	//        x a1a0
	// ---------------------------
	//         |    a0a0    ...     r1r0
	//         |a1a1        ... r3r2

	// 2*
	//         |  a1a0      ...   r2r1

	// result max. FFFF x FFFF = FFFE0001

	// load operands (value max FFFF)
	const HALFUINT a0 = ((const HALFUINT*)&num)[0];
	const HALFUINT a1 = ((const HALFUINT*)&num)[1];

	// rH
	BASEUINT rH = (BASEUINT)a1*a1;

	// rL
	BASEUINT rL = (BASEUINT)a0*a0;

	// middle
	BASEUINT r2 = (BASEUINT)a1*a0;
	BASEUINT r = r2 + (rL >> HALFBITS); // max FFFEFFFF, no carry
	r += r2;
	rH += BASEUINT((r < r2) ? 1 : 0) << HALFBITS; // add carry

	// save result
	this->N0() = (rL & HALFMASK) | (r << HALFBITS);
	this->N1() = rH + (r >> HALFBITS);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 / num2, with optional destination remainder (sources and destinations can be mixed)
// AMD 4 GHz: uint128 (ASM) 201 ns, uint128ref (C++) 336 ns
// Intel 3 GHz: uint128 (ASM): 280 ns, uint128ref (C++): 450 ns
// Intel 2 GHz: uint128 (ASM): 480 ns, uint128ref (C++): 500 ns

// this = num1 / (0:)num2
// AMD 4 GHz: uint128 (ASM) 30 ns, uint128ref (C++) 160 ns
// Intel 3 GHz: uint128 (ASM): 40 ns, uint128ref (C++): 215 ns
// Intel 2 GHz: uint128 (ASM): 64 ns, uint128ref (C++): 150 ns

void DBLUINT::Div(const DBLUINT& num1, const DBLUINT& num2, DBLUINT* rem /* = NULL */)
{
	// divide by small number
	if (num2.N1() == 0)
	{
		if (rem == NULL)
			this->Div(num1, *(const BASEUINT*)&num2);
		else
			rem->Set(this->Div(num1, *(const BASEUINT*)&num2));
		return;
	}

#ifdef USEASM
	DivU128A_x64(this, &num1, &num2, rem);
#else // USEASM

	// prepare dividend
	BASEUINT d[4];
	d[0] = num1.N0();
	d[1] = num1.N1();
	d[2] = 0;
	d[3] = 0;

	// prepare divisor
	BASEUINT n[2];
	n[0] = num2.N0();
	n[1] = num2.N1();

	// prepare accumulator
	BASEUINT a[2];
	a[0] = 0;
	a[1] = 0;

	// prepare mask
	BASEUINT m[2];
	m[0] = 0;
	m[1] = BASELAST;

	// divide
	int i;
	BASEUINT c;
	for (i = 2*BASEBITS; i > 0; i--)
	{
		// shift dividend
		d[3] = (d[3] << 1) | (d[2] >> (BASEBITS-1));
		d[2] = (d[2] << 1) | (d[1] >> (BASEBITS-1));
		d[1] = (d[1] << 1) | (d[0] >> (BASEBITS-1));
		d[0] <<= 1;

		// compare dividend with divisor
		if ((d[3] > n[1]) || ((d[3] == n[1]) && (d[2] >= n[0])))
		{
			// sub divisor from dividend
			c = (d[2] < n[0]) ? 1 : 0;
			d[2] -= n[0];
			d[3] = d[3] - n[1] - c;

			// add mask to accumulator
			a[0] |= m[0];
			a[1] |= m[1];
		}

		// shift mask
		m[0] = (m[0] >> 1) | (m[1] << (BASEBITS-1));
		m[1] >>= 1;
	}

	// save result
	this->N0() = a[0];
	this->N1() = a[1];

	// save remainder
	if (rem != NULL)
	{
		rem->N0() = d[2];
		rem->N1() = d[3];
	}
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 / 0:num2, returns remainder
// AMD 4 GHz: uint128 (ASM) 21.1 ns, uint128ref (C++) 151.8 ns
// Intel 3 GHz: uint128 (ASM): 36 ns, uint128ref (C++): 210 ns
// Intel 2 GHz: uint128 (ASM): 62 ns, uint128ref (C++): 150 ns

BASEUINT DBLUINT::Div(const DBLUINT& num1, BASEUINT num2)
{
	// divide by small number
	if (num2 <= HALFMASK) return this->Div(num1, *(const HALFUINT*)&num2);

#ifdef USEASM
	return DivU128B_x64(this, &num1, num2);
#else // USEASM

	// divide N1
	BASEUINT d = num1.N1(); // dividend
	BASEUINT a = 0; // accumulator
	BASEUINT m = BASELAST; // mask
	BASEUINT c = 0; // carry
	int n;
	if (d < num2)
	{
		c = d;
	}
	else
	{
		for (n = BASEBITS; n > 0; n--)
		{
			// add highest bit of dividend to carry
			c = (c << 1) | (d >> (BASEBITS-1));

			// compare, result bit is 1 or 0
			if (num2 <= c)
			{
				a |= m;
				c -= num2;
			}

			// shift dividend and mask
			d <<= 1;
			m >>= 1;
		}
	}
	this->N1() = a;
	
	// divide N0
	d = num1.N0(); // dividend
	a = 0; // accumulator
	m = BASELAST; // mask
	BASEUINT c2;

	if ((c == 0) && (d < num2))
	{
		c = d;
	}
	else
	{
		// divisor is full range of bits
		if (num2 >= BASELAST)
		{
			for (n = BASEBITS; n > 0; n--)
			{
				// add highest bit of dividend to carry
				c2 = c;
				c = (c << 1) | (d >> (BASEBITS-1));

				// compare, result bit is 1 or 0
				if ((num2 <= c) || ((c2 & BASELAST) != 0))
				{
					a |= m;
					c -= num2;
				}

				// shift dividend and mask
				d <<= 1;
				m >>= 1;
			}
		}

		// divisor is less than full range of bits
		else
		{
			for (n = BASEBITS; n > 0; n--)
			{
				// add highest bit of dividend to carry
				c = (c << 1) | (d >> (BASEBITS-1));

				// compare, result bit is 1 or 0
				if (num2 <= c)
				{
					a |= m;
					c -= num2;
				}

				// shift dividend and mask
				d <<= 1;
				m >>= 1;
			}
		}
	}
	this->N0() = a;

	// return carry
	return c;
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 / 0:0:0:num2, returns remainder
// AMD 4 GHz: uint128 (ASM) 21.1 ns, uint128ref (C++) 151.8 ns
// Intel 3 GHz: uint128 (ASM): 39 ns, uint128ref (C++): 54 ns
// Intel 2 GHz: uint128 (ASM): 62 ns, uint128ref (C++): 150 ns

HALFUINT DBLUINT::Div(const DBLUINT& num1, HALFUINT num2)
{
#ifdef USEASM
	return DivU128C_x64(this, &num1, num2);
#else // USEASM

	// zero divide
	if (num2 == 0)
	{
		this->SetMax();
		return 0;
	}

	// r3
	BASEUINT n = num1.N1();
	BASEUINT k1 = n >> HALFBITS;
	HALFUINT r1 = (HALFUINT)(k1/num2);
	k1 -= (BASEUINT)r1*num2;

	// r2
	BASEUINT k0 = (n & HALFMASK) | (k1 << HALFBITS);
	HALFUINT r0 = (HALFUINT)(k0/num2);
	k0 -= (BASEUINT)r0*num2;

	// save result N1
	this->N1() = ((BASEUINT)r1 << HALFBITS) | r0;

	// r1
	n = num1.N0();
	k1 = (n >> HALFBITS) | (k0 << HALFBITS);
	r1 = (HALFUINT)(k1/num2);
	k1 -= (BASEUINT)r1*num2;

	// r0
	k0 = (n & HALFMASK) | (k1 << HALFBITS);
	r0 = (HALFUINT)(k0/num2);
	k0 -= (BASEUINT)r0*num2;

	// save result N0
	this->N0() = ((BASEUINT)r1 << HALFBITS) | r0;

	// remainder
	return (HALFUINT)k0;
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num / 10, returns remainder 0..9

int DBLUINT::Div10(const DBLUINT& num)
{
#ifdef USEASM
	return DivU128E_x64(this, &num);
#else // USEASM

	// r3
	BASEUINT n = num.N1();
	BASEUINT k1 = n >> HALFBITS;
	HALFUINT r1 = (HALFUINT)(k1/10);
	k1 -= (BASEUINT)r1*10;

	// r2
	BASEUINT k0 = (n & HALFMASK) | (k1 << HALFBITS);
	HALFUINT r0 = (HALFUINT)(k0/10);
	k0 -= (BASEUINT)r0*10;

	// save result N1
	this->N1() = ((BASEUINT)r1 << HALFBITS) | r0;

	// r1
	n = num.N0();
	k1 = (n >> HALFBITS) | (k0 << HALFBITS);
	r1 = (HALFUINT)(k1/10);
	k1 -= (BASEUINT)r1*10;

	// r0
	k0 = (n & HALFMASK) | (k1 << HALFBITS);
	r0 = (HALFUINT)(k0/10);
	k0 -= (BASEUINT)r0*10;

	// save result N0
	this->N0() = ((BASEUINT)r1 << HALFBITS) | r0;

	// remainder
	return (int)k0;
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1H:num1L / num2, with full range of dividend

void DBLUINT::DivHigh(const DBLUINT& num1L, const DBLUINT& num1H, const DBLUINT& num2, DBLUINT* rem /* = NULL */)
{
	if (num1H >= num2)
	{
		this->SetMax();
		if (rem != NULL) rem->Set0();
		return;
	}

#ifdef USEASM
	DivU128D_x64(this, &num1L, &num1H, &num2, rem);
#else // USEASM

	// prepare dividend
	BASEUINT d[4];
	d[0] = num1L.N0();
	d[1] = num1L.N1();
	d[2] = num1H.N0();
	d[3] = num1H.N1();

	// prepare divisor
	BASEUINT n[2];
	n[0] = num2.N0();
	n[1] = num2.N1();

	// prepare accumulator
	BASEUINT a[2];
	a[0] = 0;
	a[1] = 0;

	// prepare mask
	BASEUINT m[2];
	m[0] = 0;
	m[1] = BASELAST;

	// divide
	int i;
	BASEUINT c;
	for (i = 2*BASEBITS; i > 0; i--)
	{
		// shift dividend
		c = d[3];
		d[3] = (d[3] << 1) | (d[2] >> (BASEBITS-1));
		d[2] = (d[2] << 1) | (d[1] >> (BASEBITS-1));
		d[1] = (d[1] << 1) | (d[0] >> (BASEBITS-1));
		d[0] <<= 1;

		// compare dividend with divisor
		if (((c & BASELAST) != 0) || (d[3] > n[1]) || ((d[3] == n[1]) && (d[2] >= n[0])))
		{
			// sub divisor from dividend
			c = (d[2] < n[0]) ? 1 : 0;
			d[2] -= n[0];
			d[3] = d[3] - n[1] - c;

			// add mask to accumulator
			a[0] |= m[0];
			a[1] |= m[1];
		}

		// shift mask
		m[0] = (m[0] >> 1) | (m[1] << (BASEBITS-1));
		m[1] >>= 1;
	}

	// save result
	this->N0() = a[0];
	this->N1() = a[1];

	// save remainder
	if (rem != NULL)
	{
		rem->N0() = d[2];
		rem->N1() = d[3];
	}

#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// precalculate multiplier for fast division, returns flags FASTDIV_*
// Divisor and "this" cannot be the same variables!

int DBLUINT::PreFastDiv(const DBLUINT& div)
{
	// special cases
	if (div.Equ0()) return FASTDIV_0;
	if (div.Equ1()) return FASTDIV_1;
	if (div.Equ(2)) return FASTDIV_2;
	if (div.IsPwr2()) return FASTDIV_SHIFT | ((div.Bits() - 1) << FASTDIV_BIT);
	if (div.GetBit(BASEBITS*2-1)==1)
	{
		this->Set(div);
		return FASTDIV_HIGH;
	}

	// number of valid bits of divisor
	DBLUINT n;
	n.Dec(div);
	int shift = n.Bits() - 1;

	// prepare dividend
	DBLUINT d[2];
	d[0].Set0();
	d[1].Set0();
	d[1].SetBit(shift);

	// divide
	this->DivHigh(d[0], d[1], div, &d[0]);

	// shift result left (forget highest bit, it has always value "1")
	this->LShift1();
	shift++;

	// shift dividend
	int c = d[0].GetBit(BASEBITS*2-1);
	d[0].LShift1();

	// add least bit
	if ((c == 1) || d[0].GreaterEqu(div))
	{
		// sub divisor from dividend
		d[0].Sub(div);

		// add mask to accumulator
		this->SetBit(0);
	}

	// check remainder, round result up if remainder is not zero
	if (d[0].NEqu0()) this->Inc();

	// reduction in case of least bit is 0
	int mode = FASTDIV_FULL;
	if ((shift > 0) && (this->GetBit(0) == 0))
	{
		// change to limited precision, we may use simple multiplication
		mode = FASTDIV_MUL;
		shift--;
		this->RShift1();
		this->SetBit(2*BASEBITS-1);

		// reduction of another "0" bits
		while ((shift > 0) && (this->GetBit(0) == 0))
		{
			shift--;
			this->RShift1();
		}
	}

	return mode | (shift << FASTDIV_BIT);
}

///////////////////////////////////////////////////////////////////////////////
// fast divide using precalculated multiplier with flags
// Dividend "num" and "this" cannot be the same variables!

void DBLUINT::FastDiv(const DBLUINT& num, const DBLUINT& mul, int flags)
{
	// special cases (without negative flag)
	switch (flags & FASTDIV_MASK2)
	{
	// divide by 0
	case FASTDIV_0:
		this->SetMax();
		return;

	// divide by 1
	case FASTDIV_1:
		this->Set(num);
		return;

	// divide by 2
	case FASTDIV_2:
		this->RShift1(num);
		return;

	// divide by power of 2
	case FASTDIV_SHIFT:
		this->RShift(num, flags >> FASTDIV_BIT);
		return;

	// divide by number with highest bit set (result is only 0 or 1)
	case FASTDIV_HIGH:
		this->Set(num.GreaterEqu(mul) ? 1 : 0);
		return;

	// divide using simple multiplication, limited precission
	case FASTDIV_MUL:
		{
			DBLUINT n;
			n.MulHigh(num, mul, this);
			this->RShift(flags >> FASTDIV_BIT);
		}
		return;

	// full precission
	//case FASTDIV_FULL:
	default:
		{
			// multiply this = num * mul
			DBLUINT n;
			n.MulHigh(num, mul, this);

			// add dividend (= multiply with highest hidden bit of mul)
			this->Add(num);
			bool c = this->Less(num); // carry

			// shift result right
			int shift = flags >> FASTDIV_BIT;
			this->RShift(shift);

			// add carry bit
			if (c) this->SetBit(2*BASEBITS - shift);
		}
		return;
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = -this

void DBLUINT::Neg()
{
#ifdef SOFTASM
	NegU128A_x64(this);
#else // SOFTASM
	BASEUINT r = 0 - this->N0();
	this->N1() = 0 - this->N1() - ((0 != r) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = -num

void DBLUINT::Neg(const DBLUINT& num)
{
#ifdef SOFTASM
	NegU128B_x64(this, &num);
#else // SOFTASM
	BASEUINT r = 0 - num.N0();
	this->N1() = 0 - num.N1() - ((0 != r) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// this = -0:num

void DBLUINT::Neg(BASEUINT num)
{
#ifdef SOFTASM
	NegU128C_x64(this, num);
#else // SOFTASM
	BASEUINT r = 0 - num;
	this->N1() = 0 - ((0 != r) ? 1 : 0);
	this->N0() = r;
#endif // SOFTASM
}

///////////////////////////////////////////////////////////////////////////////
// set bit

void DBLUINT::SetBit(int bit)
{
	if ((unsigned int)bit < 2*BASEBITS)
	{
		if (bit < BASEBITS)
			this->N0() |= (BASEUINT)1 << bit;
		else
			this->N1() |= (BASEUINT)1 << (bit - BASEBITS);
	}
}

///////////////////////////////////////////////////////////////////////////////
// reset bit

void DBLUINT::ResBit(int bit)
{
	if ((unsigned int)bit < 2*BASEBITS)
	{
		if (bit < BASEBITS)
			this->N0() &= ~((BASEUINT)1 << bit);
		else
			this->N1() &= ~((BASEUINT)1 << (bit - BASEBITS));
	}
}

///////////////////////////////////////////////////////////////////////////////
// flip bit

void DBLUINT::FlipBit(int bit)
{
	if ((unsigned int)bit < 2*BASEBITS)
	{
		if (bit < BASEBITS)
			this->N0() ^= ((BASEUINT)1 << bit);
		else
			this->N1() ^= ((BASEUINT)1 << (bit - BASEBITS));
	}
}

///////////////////////////////////////////////////////////////////////////////
// get bit (returns 0 or 1)

int DBLUINT::GetBit(int bit) const
{
	if ((unsigned int)bit < 2*BASEBITS)
	{
		if (bit < BASEBITS)
			return (int)((this->N0() >> bit) & 1);
		else
			return (int)((this->N1() >> (bit - BASEBITS)) & 1);
	}
	return 0;
}

///////////////////////////////////////////////////////////////////////////////
// get number of valid bits (=logarithm+1, returns position of highest bit + 1, 0..2*BASEBITS)

int DBLUINT::Bits() const
{
#ifdef USEASM
	return (int)BitsU128_x64(this);
#else // USEASM
	if (this->N1() == 0)
		return ::Bits64(this->N0());
	else
		return ::Bits64(this->N1()) + BASEBITS;
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// check if number is power of 2 (use Bits()-1 to determine 2nd logarithm)

bool DBLUINT::IsPwr2() const
{
	if (this->Equ0()) return false;
	DBLUINT n;
	n.Dec(*this);
	n.And(*this);
	return n.Equ0();
}

///////////////////////////////////////////////////////////////////////////////
// this <<= 1

void DBLUINT::LShift1()
{
	this->N1() = (this->N1() << 1) | (this->N0() >> (BASEBITS - 1));
	this->N0() <<= 1;
}

///////////////////////////////////////////////////////////////////////////////
// this <<= shift

void DBLUINT::LShift(int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to right or no shift
	if (shift <= 0)
	{
		if (shift < 0) this->RShift(-shift);
		return;
	}

	// shift to higher word or overload
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N1() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '<<' < BASEBITS
			this->N1() = this->N0() << (shift - BASEBITS);

		this->N0() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
	this->N1() = (this->N1() << shift) | (this->N0() >> (BASEBITS - shift));
	this->N0() <<= shift;
}

///////////////////////////////////////////////////////////////////////////////
// this = num << 1

void DBLUINT::LShift1(const DBLUINT& num)
{
	this->N1() = (num.N1() << 1) | (num.N0() >> (BASEBITS - 1));
	this->N0() = num.N0() << 1;
}

///////////////////////////////////////////////////////////////////////////////
// this = num << shift

void DBLUINT::LShift(const DBLUINT& num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to right or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->RShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	// shift to higher word
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N1() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '<<' < BASEBITS
			this->N1() = num.N0() << (shift - BASEBITS);

		this->N0() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
	this->N1() = (num.N1() << shift) | (num.N0() >> (BASEBITS - shift));
	this->N0() = num.N0() << shift;
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num << 1

void DBLUINT::LShift1(BASEUINT num)
{
	this->N1() = num >> (BASEBITS - 1);
	this->N0() = num << 1;
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num << shift

void DBLUINT::LShift(BASEUINT num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to right or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->RShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	// shift to higher word
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N1() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '<<' < BASEBITS
			this->N1() = num << (shift - BASEBITS);

		this->N0() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
	this->N1() = num >> (BASEBITS - shift);
	this->N0() = num << shift;
}

///////////////////////////////////////////////////////////////////////////////
// this >>= 1

void DBLUINT::RShift1()
{
	this->N0() = (this->N0() >> 1) | (this->N1() << (BASEBITS - 1));
	this->N1() >>= 1;
}

///////////////////////////////////////////////////////////////////////////////
// this >>= shift

void DBLUINT::RShift(int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to left or no shift
	if (shift <= 0)
	{
		if (shift < 0) this->LShift(-shift);
		return;
	}

	// shift to lower word or overload
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N0() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '>>' < BASEBITS
			this->N0() = this->N1() >> (shift - BASEBITS);

		this->N1() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
	this->N0() = (this->N0() >> shift) | (this->N1() << (BASEBITS - shift));
	this->N1() >>= shift;
}

///////////////////////////////////////////////////////////////////////////////
// this = num >> 1

void DBLUINT::RShift1(const DBLUINT& num)
{
	this->N0() = (num.N0() >> 1) | (num.N1() << (BASEBITS - 1));
	this->N1() = num.N1() >> 1;
}

///////////////////////////////////////////////////////////////////////////////
// this = num >> shift

void DBLUINT::RShift(const DBLUINT& num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to left or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->LShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	// shift to lower word or overload
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N0() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '>>' < BASEBITS
			this->N0() = num.N1() >> (shift - BASEBITS);

		this->N1() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
	this->N0() = (num.N0() >> shift) | (num.N1() << (BASEBITS - shift));
	this->N1() = num.N1() >> shift;
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num >> 1

void DBLUINT::RShift1(BASEUINT num)
{
	this->N0() = num >> 1;
	this->N1() = 0;
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num >> shift

void DBLUINT::RShift(BASEUINT num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to left or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->LShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	// overload
	if (shift >= BASEBITS)
	{
		this->N0() = 0;
		this->N1() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
	this->N0() = num >> shift;
	this->N1() = 0;
}

///////////////////////////////////////////////////////////////////////////////
// convert to ASCIIZ text into buffer of size approx. 1 + 2.41*size in bytes
// (uint128: 39+1, uin64: 20+1, uint32: 10+1), returns number of digits
// sign = convert number as signed, with extra '-' character (reserve one more place in buffer)

int DBLUINT::ToTextBuf(char* buf, int bufsize, bool sign /* = false */) const
{
	// invalid buffer
	if ((buf == NULL) || (bufsize <= 0)) return 0;

	if (bufsize == 1)
	{
		buf[0] = 0;
		return 0;
	}

	// number is zero
	if (this->Equ0())
	{
		buf[0] = '0';
		buf[1] = 0;
		return 1;
	}

	// copy number
	DBLUINT s = *this;

	// signed
	int n2 = 0;
	if (sign && s.IsNeg())
	{
		s.Neg();
		buf[0] = '-';
		buf[1] = 0;
		buf++;
		bufsize--;
		n2++;
		if (bufsize == 1) return 1;
	}

	// decode number (to end of buffer)
	bufsize--;
	int n = 0;
	char* d = &buf[bufsize-1];
	while (s.NEqu0() && (n < bufsize))
	{
		*d = s.Div10() + '0';
		d--;
		n++;
	}

	// shift text to start of buffer
	memmove(buf, &buf[bufsize - n], n);
	buf[n] = 0;

	return n + n2;
}

CText DBLUINT::ToText(bool sign /* = false */) const
{
	CText txt;
	int len = (int)(BASEBITS/8*2*2.5+3);
	if (!txt.SetLength(len)) return txt;
	len = this->ToTextBuf(txt.DataData(), len, sign);
	txt.SetLength(len);
	return txt;
}

///////////////////////////////////////////////////////////////////////////////
// convert number from text (len=text length, -1=unlimited or terminated by 0, returns number of processed characters)

int DBLUINT::FromTextBuf(const char* text, int len /* = -1 */)
{
	// clear accumulator
	this->Set0();
	char ch;
	int len0 = len;

	// skip spaces and unary sign
	bool neg = false;
	while ((len > 0) && ((ch = *text) != 0))
	{
		if ((ch == ' ') || (ch == 9) || (ch == '+'))
		{
			text++;
			len--;
		}
		else
		{
			if (ch != '-') break;
			text++;
			len--;
			neg = !neg;
		}
	}

	// load digits
	while ((len > 0) && ((ch = *text++) != 0))
	{
		if ((ch < '0') || (ch > '9')) break;
		len--;
		this->Mul10();
		this->Add((BASEUINT)(ch - '0'));
	}

	// negate
	if (neg) this->Neg();

	return len0 - len;
}

// convert number from text (returns new position)
int DBLUINT::FromText(const CText& text, int pos /* = 0 */)
{
	int len = text.Length();
	if (pos < 0) pos = 0;
	if (pos >= len)
	{
		this->Set0();
		return len;
	}
	return this->FromTextBuf(text.DataData() + pos, len - pos) + pos;
}

///////////////////////////////////////////////////////////////////////////////
// integer square root with optional remainder, this = sqrt(num)
// - using abacus algorithm
// http://freaknet.org/martin/tape/gos/misc/personal/msc/sqrt/sqrt.c
// http://www.codecodex.com/wiki/Calculate_an_integer_square_root

void DBLUINT::Sqrt(const DBLUINT& num, DBLUINT* rem /* = NULL */)
{
	// special cases
	if (num.LessEqu((BASEUINT)3))
	{
		// zero
		if (num.IsZero())
		{
			this->Set0();
			if (rem != NULL) rem->Set0();
			return;
		}

		// 1..3
		this->Set1();
		if (rem != NULL) rem->Dec(num);
		return;
	}

	// prepare registers
	DBLUINT r, p, t; // remainder, place (mask), temporary
	this->Set0(); // root = 0
	r.Set(num); // remainder = n
	p.Set0(); // place = 0 (clear bit mask)
	p.SetBit((num.Bits() - 1) & ~1); // place = bit mask (must be <= number, power of 4)

	// search square root
	do {
		t.Add(*this, p); // t = root + place
		if (r.GreaterEqu(t)) // if (remainder >= root + place)
		{
			r.Sub(t); // remainder = remainder - root - place
			t.LShift1(p); // t = (place << 1)
			this->Add(t); // root = root + (place << 1)
		}

		this->RShift1(); // root = root >> 1
		p.RShift(2); // place = place >> 2
	} while(p.IsNotZero()); // while (place > 0)

	// save remainder
	if (rem != NULL) rem->Set(r);
}

///////////////////////////////////////////////////////////////////////////////
// initialize factorial table

bool DBLSINT::FactInit()
{
	// already initialized
	if (DBLSINT::Fact_ != NULL) return true;

	// factorial value 0
	DBLSINT::Fact_ = (DBLSINT*)malloc(sizeof(DBLSINT));
	if (DBLSINT::Fact_ == NULL) return false;
	DBLSINT::Fact_->Set0();

	// loop
	DBLSINT n, n2, n3;
	n.Set1();
	BASEUINT i;
	for (i = 1; ; i++)
	{
		// get new member
		n2 = n;
		n.Mul((BASESINT)i);
		if (n.IsNeg()) break; // overflow
		n3.Div(n, (BASESINT)i);
		if (n2.NEqu(n3)) break; // overflow

		// resize table
		DBLSINT::Fact_ = (DBLSINT*)realloc(DBLSINT::Fact_, (size_t)((i+1)*sizeof(DBLSINT)));
		if (DBLSINT::Fact_ == NULL)
		{
			DBLSINT::FactMax = 0;
			return false;
		}

		// store new member
		DBLSINT::FactMax = i;
		DBLSINT::Fact_[i].Set(n);
	}

	return true;
}

///////////////////////////////////////////////////////////////////////////////
// factorial this = num!

void DBLSINT::Fact(BASEUINT num)
{
	if (num <= this->FactMax)
		this->Set(this->Fact_[num]);
	else
		this->SetMax();
}

void DBLSINT::Fact()
{
	if (this->N1() == 0)
		this->Fact(this->N0());
	else
		if (this->N1() >= 0)
			this->SetMax();
		else
			this->Set0();
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * num2

void DBLSINT::Mul(const DBLSINT& num1, const DBLSINT& num2)
{
	if (num1.IsNeg())
	{
		DBLUINT n1;
		n1.Neg(*(const DBLUINT*)&num1);

		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			DBLUINT::Mul(n1, n2);
		}
		else
		{
			DBLUINT::Mul(n1, *(const DBLUINT*)&num2);
			this->Neg();
		}
	}
	else
	{
		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			DBLUINT::Mul(*(const DBLUINT*)&num1, n2);
			this->Neg();
		}
		else
		{
			DBLUINT::Mul(*(const DBLUINT*)&num1, *(const DBLUINT*)&num2);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * s:num2

void DBLSINT::Mul(const DBLSINT& num1, BASESINT num2)
{
	if (num1.IsNeg())
	{
		DBLUINT n1;
		n1.Neg(*(const DBLUINT*)&num1);

		if (num2 < 0)
		{
			DBLUINT::Mul(n1, (BASEUINT)-num2);
		}
		else
		{
			DBLUINT::Mul(n1, (BASEUINT)num2);
			this->Neg();
		}
	}
	else
	{
		if (num2 < 0)
		{
			DBLUINT::Mul(*(const DBLUINT*)&num1, (BASEUINT)-num2);
			this->Neg();
		}
		else
		{
			DBLUINT::Mul(*(const DBLUINT*)&num1, (BASEUINT)num2);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = s:num1 * s:num2

void DBLSINT::Mul(BASESINT num1, BASESINT num2)
{
	if (num1 < 0)
	{
		if (num2 < 0)
		{
			DBLUINT::Mul((BASEUINT)-num1, (BASEUINT)-num2);
		}
		else
		{
			DBLUINT::Mul((BASEUINT)-num1, (BASEUINT)num2);
			this->Neg();
		}
	}
	else
	{
		if (num2 < 0)
		{
			DBLUINT::Mul((BASEUINT)num1, (BASEUINT)-num2);
			this->Neg();
		}
		else
		{
			DBLUINT::Mul((BASEUINT)num1, (BASEUINT)num2);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * num2, with full range of result (high = destination of high part of result)
// "this" result is unsigned! (sign is in high part of result)

void DBLSINT::MulHigh(const DBLSINT& num1, const DBLSINT& num2, DBLSINT* high)
{
	if (num1.IsNeg())
	{
		DBLUINT n1;
		n1.Neg(*(const DBLUINT*)&num1);

		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			DBLUINT::MulHigh(n1, n2, high);
		}
		else
		{
			DBLUINT::MulHigh(n1, *(const DBLUINT*)&num2, high);
			this->Neg();
			high->Neg();
			if (this->IsNotZero()) high->Dec();
		}
	}
	else
	{
		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			DBLUINT::MulHigh(*(const DBLUINT*)&num1, n2, high);
			this->Neg();
			high->Neg();
			if (this->IsNotZero()) high->Dec();
		}
		else
		{
			DBLUINT::MulHigh(*(const DBLUINT*)&num1, *(const DBLUINT*)&num2, high);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = num * 10, returns high carry (+-0..9)

int DBLSINT::Mul10(const DBLSINT& num)
{
	if (num.IsNeg())
	{
		DBLUINT n;
		n.Neg(*(const DBLUINT*)&num);
		int res = DBLUINT::Mul10(n);
		this->Neg();
		return this->IsZero() ? -res : (-res-1);
	}
	else
		return DBLUINT::Mul10(*(const DBLUINT*)&num);
}

///////////////////////////////////////////////////////////////////////////////
// this = num * num

void DBLSINT::Sqr(const DBLSINT& num)
{
	if (num.IsNeg())
	{
		DBLUINT n;
		n.Neg(*(const DBLUINT*)&num);
		DBLUINT::Sqr(n);
	}
	else
		DBLUINT::Sqr(*(const DBLUINT*)&num);
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 / num2, with optional destination remainder (sources and destinations can be mixed)

void DBLSINT::Div(const DBLSINT& num1, const DBLSINT& num2, DBLSINT* rem /* = NULL */)
{
	if (num1.IsNeg())
	{
		DBLUINT n1;
		n1.Neg(*(const DBLUINT*)&num1);

		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			DBLUINT::Div(n1, n2, (DBLUINT*)rem);
		}
		else
		{
			DBLUINT::Div(n1, *(const DBLUINT*)&num2, (DBLUINT*)rem);
			this->Neg();
		}

		if (rem != NULL) rem->Neg();
	}
	else
	{
		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			DBLUINT::Div(*(const DBLUINT*)&num1, n2, (DBLUINT*)rem);
			this->Neg();
		}
		else
		{
			DBLUINT::Div(*(const DBLUINT*)&num1, *(const DBLUINT*)&num2, (DBLUINT*)rem);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 / s:num2, returns remainder

BASESINT DBLSINT::Div(const DBLSINT& num1, BASESINT num2)
{
	if (num1.IsNeg())
	{
		DBLUINT n1;
		n1.Neg(*(const DBLUINT*)&num1);

		if (num2 < 0)
		{
			return -(BASESINT)DBLUINT::Div(n1, (BASEUINT)-num2);
		}
		else
		{
			BASESINT res = (BASESINT)DBLUINT::Div(n1, (BASEUINT)num2);
			this->Neg();
			return -res;
		}
	}
	else
	{
		if (num2 < 0)
		{
			BASESINT res = (BASESINT)DBLUINT::Div(*(const DBLUINT*)&num1, (BASEUINT)-num2);
			this->Neg();
			return res;
		}
		else
		{
			return (BASESINT)DBLUINT::Div(*(const DBLUINT*)&num1, (BASEUINT)num2);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 / s:s:s:num2, returns remainder

HALFSINT DBLSINT::Div(const DBLSINT& num1, HALFSINT num2)
{
	if (num1.IsNeg())
	{
		DBLUINT n1;
		n1.Neg(*(const DBLUINT*)&num1);

		if (num2 < 0)
		{
			return -(HALFSINT)DBLUINT::Div(n1, (HALFUINT)-num2);
		}
		else
		{
			HALFSINT res = (HALFSINT)DBLUINT::Div(n1, (HALFUINT)num2);
			this->Neg();
			return -res;
		}
	}
	else
	{
		if (num2 < 0)
		{
			HALFSINT res = (HALFSINT)DBLUINT::Div(*(const DBLUINT*)&num1, (HALFUINT)-num2);
			this->Neg();
			return res;
		}
		else
		{
			return (HALFSINT)DBLUINT::Div(*(const DBLUINT*)&num1, (HALFUINT)num2);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = num1H:num1L / num2, with full range of dividend
// remainder will have the same sign as dividend

void DBLSINT::DivHigh(const DBLUINT& num1L, const DBLSINT& num1H, const DBLSINT& num2, DBLSINT* rem /* = NULL */)
{
	if (num1H.IsNeg())
	{
		DBLUINT n1L, n1H;
		n1L.Neg(num1L);
		n1H.Neg(*(const DBLUINT*)&num1H);
		if (n1L.IsNotZero()) n1H.Dec();

		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			if (n1H >= n2)
			{
				this->SetMax();
				if (rem != NULL) rem->Set0();
				return;
			}
			DBLUINT::DivHigh(n1L, n1H, n2, (DBLUINT*)rem);
		}
		else
		{
			if (n1H >= *(const DBLUINT*)&num2)
			{
				this->SetMin();
				if (rem != NULL) rem->Set0();
				return;
			}
			DBLUINT::DivHigh(n1L, n1H, *(const DBLUINT*)&num2, (DBLUINT*)rem);
			this->Neg();
		}

		if (rem != NULL) rem->Neg();
	}
	else
	{
		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			if (*(const DBLUINT*)&num1H >= n2)
			{
				this->SetMin();
				if (rem != NULL) rem->Set0();
				return;
			}
			DBLUINT::DivHigh(num1L, *(const DBLUINT*)&num1H, n2, (DBLUINT*)rem);
			this->Neg();
		}
		else
		{
			if (*(const DBLUINT*)&num1H >= *(const DBLUINT*)&num2)
			{
				this->SetMax();
				if (rem != NULL) rem->Set0();
				return;
			}
			DBLUINT::DivHigh(num1L, *(const DBLUINT*)&num1H, *(const DBLUINT*)&num2, (DBLUINT*)rem);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = num / 10, returns remainder +-0..9

int DBLSINT::Div10(const DBLSINT& num)
{
	if (num.IsNeg())
	{
		DBLUINT n;
		n.Neg(*(const DBLUINT*)&num);
		int res = DBLUINT::Div10(n);
		this->Neg();
		return -res;
	}
	else
		return DBLUINT::Div10(*(const DBLUINT*)&num);
}

///////////////////////////////////////////////////////////////////////////////
// precalculate multiplier for fast division, returns flags FASTDIV_*
// Divisor and "this" cannot be the same variables!

int DBLSINT::PreFastDiv(const DBLSINT& div)
{
	if (div.IsNeg())
	{
		DBLUINT n;
		n.Neg(*(const DBLUINT*)&div);
		return DBLUINT::PreFastDiv(n) | FASTDIV_NEG;
	}
	else
		return DBLUINT::PreFastDiv(div);
}

///////////////////////////////////////////////////////////////////////////////
// fast divide using precalculated multiplier with flags
// Dividend "num" and "this" cannot be the same variables!

void DBLSINT::FastDiv(const DBLSINT& num, const DBLSINT& mul, int flags)
{
	if (num.IsNeg())
	{
		DBLUINT n;
		n.Neg(*(const DBLUINT*)&num);
		DBLUINT::FastDiv(n, *(const DBLUINT*)&mul, flags);
		if ((flags & FASTDIV_NEG) == 0) this->Neg();
	}
	else
	{
		DBLUINT::FastDiv(*(const DBLUINT*)&num, *(const DBLUINT*)&mul, flags);
		if ((flags & FASTDIV_NEG) != 0) this->Neg();
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = -s:num

void DBLSINT::Neg(BASESINT num)
{
	BASEUINT r = 0 - num;
	BASESINT k = (num < 0) ? -1 : 0;
	this->N1() = 0 - k - ((0 != r) ? 1 : 0);
	this->N0() = r;
}

///////////////////////////////////////////////////////////////////////////////
// get number of valid bits of absolute value (=logarithm+1, returns position of highest bit + 1, 0..2*BASEBITS)

int DBLSINT::Bits() const
{
	if (this->IsNeg())
	{
		DBLUINT num;
		num.Neg(*(const DBLUINT*)this);
		return num.Bits();
	}
	else
		return DBLUINT::Bits();
}

///////////////////////////////////////////////////////////////////////////////
// check if absolute value of number is power of 2 (use Bits()-1 to determine 2nd logarithm)

bool DBLSINT::IsPwr2() const
{
	if (this->IsNeg())
	{
		DBLUINT num;
		num.Neg(*(const DBLUINT*)this);
		return num.IsPwr2();
	}
	else
		return DBLUINT::IsPwr2();
}

///////////////////////////////////////////////////////////////////////////////
// this <<= shift
// Note: Do not use DBLUINT::LShift, this one is calling DBLSINT functions.

void DBLSINT::LShift(int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to right or no shift
	if (shift <= 0)
	{
		if (shift < 0) this->RShift(-shift);
		return;
	}

	// shift to higher word or overload
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N1() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '<<' < BASEBITS
			this->N1() = this->N0() << (shift - BASEBITS);

		this->N0() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
	this->N1() = (this->N1() << shift) | (this->N0() >> (BASEBITS - shift));
	this->N0() <<= shift;
}

///////////////////////////////////////////////////////////////////////////////
// this = num << shift
// Note: Do not use DBLUINT::LShift, this one is calling DBLSINT functions.

void DBLSINT::LShift(const DBLSINT& num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to right or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->RShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	// shift to higher word
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N1() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '<<' < BASEBITS
			this->N1() = num.N0() << (shift - BASEBITS);

		this->N0() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
	this->N1() = (num.N1() << shift) | (num.N0() >> (BASEBITS - shift));
	this->N0() = num.N0() << shift;
}

///////////////////////////////////////////////////////////////////////////////
// this = s:num << 1

void DBLSINT::LShift1(BASESINT num)
{
	BASESINT k = (num < 0) ? -1 : 0;
	this->N1() = ((BASEUINT)num >> (BASEBITS - 1)) | (k << 1);
	this->N0() = num << 1;
}

///////////////////////////////////////////////////////////////////////////////
// this = s:num << shift

void DBLSINT::LShift(BASESINT num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to right or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->RShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	// shift to higher word
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N1() = 0;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '<<' < BASEBITS
			this->N1() = num << (shift - BASEBITS);

		this->N0() = 0;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
	BASESINT k = (num < 0) ? -1 : 0;
	this->N1() = ((BASEUINT)num >> (BASEBITS - shift)) | (k << shift);
	this->N0() = num << shift;
}

///////////////////////////////////////////////////////////////////////////////
// this >>= 1

void DBLSINT::RShift1()
{
	this->N0() = (this->N0() >> 1) | (this->N1() << (BASEBITS - 1));
	BASESINT k = this->IsNeg() ? ((BASEUINT)1 << (BASEBITS-1)) : 0;
	this->N1() = ((BASEUINT)this->N1() >> 1) | k;
}

///////////////////////////////////////////////////////////////////////////////
// this >>= shift

void DBLSINT::RShift(int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to left or no shift
	if (shift <= 0)
	{
		if (shift < 0) this->LShift(-shift);
		return;
	}

	BASESINT k = this->IsNeg() ? -1 : 0;

	// shift to lower word or overload
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N0() = k;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '>>' < BASEBITS
			if (shift == BASEBITS)
				this->N0() = (BASEUINT)this->N1();
			else
				this->N0() = ((BASEUINT)this->N1() >> (shift - BASEBITS)) | (k << (2*BASEBITS - shift));

		this->N1() = k;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
	this->N0() = (this->N0() >> shift) | (this->N1() << (BASEBITS - shift));
	this->N1() = ((BASEUINT)this->N1() >> shift) | (k << (BASEBITS - shift));
}

///////////////////////////////////////////////////////////////////////////////
// this = num >> 1

void DBLSINT::RShift1(const DBLSINT& num)
{
	this->N0() = (num.N0() >> 1) | (num.N1() << (BASEBITS - 1));
	BASESINT k = num.IsNeg() ? ((BASEUINT)1 << (BASEBITS-1)) : 0;
	this->N1() = ((BASEUINT)num.N1() >> 1) | k;
}

///////////////////////////////////////////////////////////////////////////////
// this = num >> shift

void DBLSINT::RShift(const DBLSINT& num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to left or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->LShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	BASESINT k = num.IsNeg() ? -1 : 0;

	// shift to lower word or overload
	if (shift >= BASEBITS)
	{
		// overload
		if (shift >= 2*BASEBITS)
			this->N0() = k;
		else
			// here: BASEBITS <= shift < 2*BASEBITS
			//       0 <= '>>' < BASEBITS
			if (shift == BASEBITS)
				this->N0() = (BASEUINT)num.N1();
			else
				this->N0() = ((BASEUINT)num.N1() >> (shift - BASEBITS)) | (k << (2*BASEBITS - shift));

		this->N1() = k;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
	this->N0() = (num.N0() >> shift) | (num.N1() << (BASEBITS - shift));
	this->N1() = ((BASEUINT)num.N1() >> shift) | (k << (BASEBITS - shift));
}

///////////////////////////////////////////////////////////////////////////////
// this = s:num >> 1

void DBLSINT::RShift1(BASESINT num)
{
	BASESINT k = (num < 0) ? -1 : 0;
	this->N0() = ((BASEUINT)num >> 1) | (k << (BASEBITS - 1));
	this->N1() = k;
}

///////////////////////////////////////////////////////////////////////////////
// this = s:num >> shift

void DBLSINT::RShift(BASESINT num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (2*BASEBITS - 1);
#endif // HWSHIFT

	// shift to left or no shift
	if (shift <= 0)
	{
		if (shift < 0)
			this->LShift(num, -shift);
		else
			this->Set(num);
		return;
	}

	BASESINT k = (num < 0) ? -1 : 0;

	// overload
	if (shift >= BASEBITS)
	{
		this->N0() = k;
		this->N1() = k;
		return;
	}

	// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
	this->N0() = ((BASEUINT)num >> shift) | (k << (BASEBITS - shift));
	this->N1() = k;
}

///////////////////////////////////////////////////////////////////////////////
// initialize factorial table

bool DBLFIX::FactInit()
{
	// already initialized
	if (DBLFIX::Fact_ != NULL) return true;

	// factorial value 0
	DBLFIX::Fact_ = (DBLFIX*)malloc(sizeof(DBLFIX));
	if (DBLFIX::Fact_ == NULL) return false;
	DBLFIX::Fact_->Set0();

	DBLFIX::FactRec_ = (DBLFIX*)malloc(sizeof(DBLFIX));
	if (DBLFIX::FactRec_ == NULL) return false;
	DBLFIX::FactRec_->SetMax();

	// loop
	DBLFIX n, n2, n3;
	n.Set1();
	BASEUINT i;
	for (i = 1; ; i++)
	{
		// get new member
		n2 = n;
		n.Mul((BASESINT)i);
		if (n.IsNeg()) break; // overflow
		n3.Div(n, (BASESINT)i);
		if (n2.NEqu(n3)) break; // overflow

		// resize table
		DBLFIX::Fact_ = (DBLFIX*)realloc(DBLFIX::Fact_, (size_t)((i+1)*sizeof(DBLFIX)));
		DBLFIX::FactRec_ = (DBLFIX*)realloc(DBLFIX::FactRec_, (size_t)((i+1)*sizeof(DBLFIX)));
		if ((DBLFIX::Fact_ == NULL) || (DBLFIX::FactRec_ == NULL))
		{
			DBLFIX::FactMax = 0;
			return false;
		}

		// store new member
		DBLFIX::FactMax = i;
		DBLFIX::Fact_[i].Set(n);
		DBLFIX::FactRec_[i].Recip(n);
	}

	return true;
}

///////////////////////////////////////////////////////////////////////////////
// factorial this = num!

void DBLFIX::Fact(BASEUINT num)
{
	if (num <= this->FactMax)
		this->Set(this->Fact_[num]);
	else
		this->SetMax();
}

///////////////////////////////////////////////////////////////////////////////
// reciprocal factorial this = 1/num!

void DBLFIX::FactRec(BASEUINT num)
{
	if (num <= this->FactMax)
		this->Set(this->FactRec_[num]);
	else
		this->Set0();
}

///////////////////////////////////////////////////////////////////////////////
// export to double number

double DBLFIX::Double()
{
	return (double)this->N1() + (double)this->N0()/BASELAST/2;
}

///////////////////////////////////////////////////////////////////////////////
// export to float number

float DBLFIX::Float()
{
	return (float)this->N1() + (float)this->N0()/BASELAST/2;
}

///////////////////////////////////////////////////////////////////////////////
// this = double

void DBLFIX::Set(double num)
{
	bool neg = false;

	// negative number
	if (num < 0)
	{
		// negative overflow
		if (num <= DBLFIX_MIN)
		{
			this->SetMin();
			return;
		}

		// absolute value
		neg = true;
		num = -num;
	}
	else
	{
		// positive overflow
		if (num >= DBLFIX_MAX)
		{
			this->SetMax();
			return;
		}
	}

	// N1
	BASEUINT k = (BASEUINT)num;
	if (k >= BASELAST) k = BASELAST-1;
	num -= k;
	if (num < 0)
	{
		num += 1;
		k--;
	}
	this->N1() = (BASESINT)k;

	// N0
	// Warning: std software library limits conversion from double to s64; in Win64 mode it works OK
	num *= BASELAST;
	num *= 2;
	k = (BASEUINT)num;
	num -= k;
	BASEUINT k2 = (BASEUINT)num;
	this->N0() = k + k2;

	// negate
	if (neg) this->Neg();
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * num2

void DBLFIX::Mul(const DBLFIX& num1, const DBLFIX& num2)
{
	DBLSINT n1, n2;
	n1.MulHigh(*(const DBLSINT*)&num1, *(const DBLSINT*)&num2, &n2);
	this->N0() = n1.N1();
	this->N1() = n2.N0();
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 / num2

void DBLFIX::Div(const DBLFIX& num1, const DBLFIX& num2)
{
	DBLSINT n1, n2;
	n2.N1() = num1.IsNeg() ? -1 : 0;
	n2.N0() = num1.N1();
	n1.N1() = num1.N0();
	n1.N0() = 0;
	DBLSINT::DivHigh(n1, n2, num2);
}

///////////////////////////////////////////////////////////////////////////////
// convert to ASCIIZ text into buffer of size approx. 2 + 2.41*size in bytes
// (uint128: 39+2, uin64: 20+2, uint32: 10+2), returns number of characters
// dec = max. number of decimal places (max. DBLFIX_DEC+2)
//   default precision: fix128: 18, fix64: 8, fix32: 3
//   max. precision: fix128: 20, fix64: 10, fix32: 5 (last digit can be incorrect, previous one can be incorrectly rounded)

int DBLFIX::ToTextBuf(char* buf, int bufsize, int dec /* = DBLFIX_DEC */) const
{
	// invalid buffer (size 0 or 1)
	if ((buf == NULL) || (bufsize <= 0)) return 0;

	if (bufsize == 1)
	{
		buf[0] = 0;
		return 0;
	}

	// number is zero
	if (this->Equ0())
	{
		buf[0] = '0';
		buf[1] = 0;
		return 1;
	}

	// limit number of decimal places
	if (dec < 0) dec = 0;
	if (dec > DBLFIX_DEC+1) dec = DBLFIX_DEC+1;

	// copy number
	DBLUINT s = *(DBLUINT*)this;

	// sign
	int n2 = 0;
	if (s.IsNeg())
	{
		s.Neg();
		buf[0] = '-';
		buf[1] = 0;
		buf++;
		bufsize--;
		n2++;
		if (bufsize == 1) return 1;
	}

	// decode integer part of the number (to end of buffer)
	BASEUINT s2 = s.N1();
	bufsize--; // reserve space for terminating zero
	int n = 0; // counter of digits
	char* d = &buf[bufsize];
	BASEUINT k;
	while (n < bufsize)
	{
		d--;
		k = s2/10;
		*d = (char)(s2 - k*10 + '0');
		s2 = k;
		n++;
		if (s2 == 0) break;
	}

	// shift text to start of buffer
	memmove(buf, d, n);

	// max. length of fractional part
	dec++; // including decimal point
	if (dec > bufsize - n) dec = bufsize - n;

	// store decimal point
	int n3 = 0;
	d = &buf[n];
	if (dec > 0)
	{
		*d++ = '.';
		n3++;
	}

	// decode decimal part of the number
	while (n3 < dec)
	{
		s.N1() = 0;
		s.Mul10();
		*d++ = (char)(s.N1() + '0');
		n3++;
	}

	// round number up
	if ((BASESINT)s.N0() < 0)
	{
		int i;
		for (i = n + n3; i > 0; i--)
		{
			d--;
			if (*d == '.') continue;
			*d = *d + 1;
			if (*d <= '9') break;
			*d = '0';
		}

		// overflow, add new '1' digit
		if (i == 0)
		{
			n++;
			if (n3 >= bufsize - n) n3--;
			memmove(buf+1, buf, n3 + n);
			buf[0] = '1';
		}
	}

	// reduce trailing zeroes and '.'
	while ((n3 > 0) && ((buf[n+n3-1] == '.') || (buf[n+n3-1] == '0'))) n3--;

	// mark end of text
	buf[n+n3] = 0;

	// return length of text
	return n2 + n + n3;
}

CText DBLFIX::ToText(int dec /* = DBLFIX_DEC */) const
{
	CText txt;
	int len = (int)(BASEBITS/8*2*2.5+5);
	if (!txt.SetLength(len)) return txt;
	len = this->ToTextBuf(txt.DataData(), len, dec);
	txt.SetLength(len);
	return txt;
}

///////////////////////////////////////////////////////////////////////////////
// convert number from text (len=text length, -1=unlimited or terminated by 0, returns number of processed characters)

int DBLFIX::FromTextBuf(const char* text, int len /* = -1 */)
{
	// clear accumulator
	this->Set0();
	char ch;
	int len0 = len;

	// skip spaces and unary sign
	bool neg = false;
	while ((len > 0) && ((ch = *text) != 0))
	{
		if ((ch == ' ') || (ch == 9) || (ch == '+'))
		{
			text++;
			len--;
		}
		else
		{
			if (ch != '-') break;
			text++;
			len--;
			neg = !neg;
		}
	}

	// load digits of integer part
	BASEUINT s = 0;
	ch = 0;
	while ((len > 0) && ((ch = *text++) != 0))
	{
		if ((ch < '0') || (ch > '9')) break;
		len--;
		s = s*10 + (ch - '0');
	}
	
	// load decimal part
	if ((len > 0) && ((ch == '.') || (ch == ',')))
	{
		len--;

		// get decimal part length
		int n = 0;
		while ((len > 0) && ((ch = *text) != 0))
		{
			if ((ch < '0') || (ch > '9')) break;
			len--;
			text++;
			n++;
		}

		// decode decimal part
		for (; n > 0; n--)
		{
			text--;
			ch = *text;
			this->N1() = ch - '0';
			this->Div10();
		}
	}

	// store integer part
	this->N1() = s;

	// negate
	if (neg) this->Neg();

	return len0 - len;
}

// convert number from text (returns new position)
int DBLFIX::FromText(const CText& text, int pos /* = 0 */)
{
	int len = text.Length();
	if (pos < 0) pos = 0;
	if (pos >= len)
	{
		this->Set0();
		return len;
	}
	return this->FromTextBuf(text.DataData() + pos, len - pos) + pos;
}

///////////////////////////////////////////////////////////////////////////////
// initialize static constants

bool DBLFIX::ConstInit()
{
	BASESINT i;
	DBLFIX x, y, z;

	// Euler's constant, e = 1 + 1/1! + 1/2! + 1/3! + 1/4! + ...
	// e = 2.7182 81828 4590452353 6028747135266249775724709369995957496696762772407663035354759...
	// https://www.math.utah.edu/~pa/math/e.html
#define SHFT (BASEBITS/2) // higher precision
	x = (BASESINT)((BASESINT)2 << SHFT); // accumulator
	y = (BASESINT)((BASESINT)1 << SHFT); // member of sequence
	i = 1; // counter
	for (;;) // loops fix32: 11, fix64: 17, fix128: 28
	{
		i++;
		y.Div(i);
		if (y.IsZero()) break;
		x.Add(y);
	}

	DBLFIX::Eul.RShift(x,SHFT - SHFT/2); // << SHFT/4
	DBLFIX::Eul2.Sqr(DBLFIX::Eul); // << SHFT/2 * 2
	DBLFIX::Eul4.RShift(DBLFIX::Eul2, SHFT/2);
	DBLFIX::Eul4.Sqr();// << SHFT/2 * 2
	DBLFIX::Eul8.RShift(DBLFIX::Eul4, SHFT/2);
	DBLFIX::Eul8.Sqr(); // << SHFT/2 * 2

	DBLFIX::Eul.RShift(SHFT/2);
	DBLFIX::Eul2.RShift(SHFT/2 * 2);
	DBLFIX::Eul4.RShift(SHFT/2 * 2);
	DBLFIX::Eul8.RShift(SHFT/2 * 2);

	DBLFIX::Eul1.Recip(DBLFIX::Eul);
	DBLFIX::Eul12.Recip(DBLFIX::Eul2);
	DBLFIX::Eul14.Recip(DBLFIX::Eul4);
	DBLFIX::Eul18.Recip(DBLFIX::Eul8);

	DBLFIX::EulTxt = DBLFIX::Eul.ToText(); // decimal places fix32: 4, fix64: 9, fix128: 19
	// fix32: 2.7183
	// fix64: 2.718281828
	// fix128: 2.7182818284590452353
	DBLFIX::Eul2Txt = DBLFIX::Eul2.ToText();
	DBLFIX::Eul4Txt = DBLFIX::Eul4.ToText();
	DBLFIX::Eul8Txt = DBLFIX::Eul8.ToText();
	DBLFIX::Eul1Txt = DBLFIX::Eul1.ToText();
	DBLFIX::Eul12Txt = DBLFIX::Eul12.ToText();
	DBLFIX::Eul14Txt = DBLFIX::Eul14.ToText();
	DBLFIX::Eul18Txt = DBLFIX::Eul18.ToText();

#undef SHFT

	// Pi constant, pi = 3 + 1/(2^3) + 3*3/(2^5*4*5) + 3*3*5/(2^7*4*6*7) + 3*3*5*7/(2^9*4*6*8*9)...
	// Isaac Newton method, using formula pi = 6*arcsin(1/2)
	// pi = 3.1415 92653 5897932384 6264338327950288419716939937510582097494459230781640628620899...
	// https://www.math.utah.edu/~pa/math/pi.html
#define SHFT (BASEBITS/2) // higher precision
	x = (BASESINT)0; // accumulator
	y = (BASESINT)((BASESINT)3 << SHFT); // member of sequence
	i = 1; // counter
	for (;;) // loops fix32: 21, fix64: 43, fix128: 89
	{
		// add member to accumulator
		x.Add(y);

		// multiply member by i*i (square of counter)
		// ... number of loops is low enough to not overflow BASESINT
		y.Mul((BASESINT)(i*i));

		// increase counter
		i += 2;

		// divide member by (i+1)*(i+2)*4
		// ... number of loops is low enough to not overflow BASESINT
		y.Div((BASESINT)(i*(i-1)*4));

		// end of aproximation
		if (y.IsZero()) break;
	}

	// Pi/180
	DBLFIX::Pi180.Div(x, (BASESINT)180);
	DBLFIX::Pi180.RShift(SHFT);
	DBLFIX::Pi180Txt = DBLFIX::Pi180.ToText();

	// 180/Pi
	DBLFIX::Pi180R.Div((BASESINT)((BASESINT)180 << SHFT), x);
	DBLFIX::Pi180RTxt = DBLFIX::Pi180R.ToText();

	// Pi*2
	DBLFIX::Pi2.RShift(x, SHFT-1);
	DBLFIX::Pi2Txt = DBLFIX::Pi2.ToText();

	// Pi
	DBLFIX::Pi.RShift(x, SHFT);
	DBLFIX::PiTxt = DBLFIX::Pi.ToText(); // decimal places fix32: 4, fix64: 9, fix128: 19
	// fix32: 3.1416
	// fix64: 3.141592653
	// fix128: 3.1415926535897932385

	// Pi/2
	DBLFIX::Pi12.RShift(x, SHFT+1);
	DBLFIX::Pi12Txt = DBLFIX::Pi12.ToText();

	// Pi/4
	DBLFIX::Pi14.RShift(x, SHFT+2);
	DBLFIX::Pi14Txt = DBLFIX::Pi14.ToText();

#undef SHFT

	// ln(10) - our calculation is not precise enough, so we must use precalculated constant
	// https://oeis.org/A002392
	CText t = "2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599";
	DBLFIX::Ln10.FromText(t);
	DBLFIX::Ln10Txt = DBLFIX::Ln10.ToText();

	DBLFIX::Ln110.Recip(DBLFIX::Ln10);
	DBLFIX::Ln110Txt = DBLFIX::Ln110.ToText();

	// 16*ln(2) - our calculation is not precise enough, so we must use precalculated constant
	t = "11.090354888959124950675713943330825089208002";
	DBLFIX::Ln2_16.FromText(t);
	DBLFIX::Ln2_16Txt = DBLFIX::Ln2_16.ToText();

	// ln(2)
	DBLFIX::Ln2.RShift(DBLFIX::Ln2_16, 4);
	DBLFIX::Ln2Txt = DBLFIX::Ln2.ToText();

	DBLFIX::Sqrt2.Set((BASESINT)2);
	DBLFIX::Sqrt2.Sqrt();
	DBLFIX::Sqrt2Txt = DBLFIX::Sqrt2.ToText();

	DBLFIX::Sqrt12.Recip((BASESINT)2);
	DBLFIX::Sqrt12.Sqrt();
	DBLFIX::Sqrt12Txt = DBLFIX::Sqrt12.ToText();

	return true;
}

///////////////////////////////////////////////////////////////////////////////
// square root, this = sqrt(num)
// - using Newton's method

void DBLFIX::Sqrt(const DBLFIX& num)
{
	// zero or negative
	if (num.IsZero() || num.IsNeg())
	{
		this->Set0();
		return;
	}

	// prepare registers
	DBLFIX n, x, d;
	n.Set(num); // copy number (needed because 'this' can be equal to 'num')

	// first guess of result
	this->Set0();
	this->SetBit((n.Bits() + BASEBITS)/2);

	// find root (usually around 5 iterations)
	for (;;)
	{
		// one step of iteration
		x.Div(n, *this); // x = n/this
		x.Add(*this); // x = this + n/this
		x.RShift1(); // x = (this + n/this)/2

		// difference
		d.Sub(x, *this);

		// end of iteration
		this->Set(x);
		if (d.IsZero()) break;
	}
}

///////////////////////////////////////////////////////////////////////////////
// normalize angle in radians into range 0..Pi*2 (including 0, excluding Pi*2)

void DBLFIX::NormRad(const DBLFIX& num)
{
	// get number
	this->Set(num);

	// small overflow (-10*Pi .. 10*Pi)
	if ((this->N1() < 3*11) && (this->N1() > -3*11))
	{
		while (this->IsNeg()) this->Add(DBLFIX::Pi2);
		while (this->GreaterEqu(DBLFIX::Pi2)) this->Sub(DBLFIX::Pi2);
	}
	else
	{
		this->Mod(DBLFIX::Pi2);
		if (this->IsNeg()) this->Add(DBLFIX::Pi2);
	}
}

///////////////////////////////////////////////////////////////////////////////
// normalize angle in degrees into range 0..360 (including 0, excluding 360)

void DBLFIX::NormDeg(const DBLFIX& num)
{
	// get number
	this->Set(num);

	// small overflow (-10*360 .. 10*360)
	if ((this->N1() < 360*10) && (this->N1() > -360*10))
	{
		while (this->IsNeg()) this->Add((BASESINT)360);
		while (this->GreaterEqu((BASESINT)360)) this->Sub((BASESINT)360);
	}
	else
	{
		this->Mod((BASESINT)360);
		if (this->IsNeg()) this->Add((BASESINT)360);
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = sin(num), sine in radiands, result -1..+1
//   SIN x = x/1! - x^3/3! + x^5/5! - x^7/7! + x^9/9! - ...

void DBLFIX::Sin(const DBLFIX& num)
{
	// normalize angle
	DBLFIX n;
	n.NormRad(num);

	// normalize to range 0..Pi
	bool neg = false;
	if (n.Greater(DBLFIX::Pi))
	{
		n.InvSub(DBLFIX::Pi2);
		neg = true;
	}

	// square
	DBLFIX sqr;
	sqr.Sqr(n);

	// prepare accumulator
	this->Set(n);

	// loop (1 to 8 loops)
	BASESINT i;
	for (i = 2; ;)
	{
		// multiply increment by square
		n.Mul(sqr);

		// divide increment by factorial
		n.Div((BASESINT)(i*(i+1)));

		// increment = 0
		if (n.IsZero()) break;

		// add increment to accumulator
		this->Sub(n);

		// increase counter
		i += 2;

		// multiply increment by square
		n.Mul(sqr);

		// divide increment by factorial
		n.Div((BASESINT)(i*(i+1)));

		// increment = 0
		if (n.IsZero()) break;

		// add increment to accumulator
		this->Add(n);

		// increase counter
		i += 2;
	}

	// negate result
	if (neg) this->Neg();
}

///////////////////////////////////////////////////////////////////////////////
// this = arcsin(num), input -1..+1, result in radians
//   ARCSIN = x + x^3/2/3 + 3*x^5/2/4/5 + 3*5*x^7/2/4/6/7 +

void DBLFIX::ASin0(const DBLFIX& num, bool arccos)
{
	// get number
	DBLFIX n;
	n.Set(num);

	// do arcsin/arccos correction
	bool arccos2 = arccos;

	// absolute value
	bool neg = false;
	if (n.IsNeg())
	{
		n.Neg();
		neg = true;
	}
	
	// overflow or 1
	if (n.GreaterEqu((BASESINT)1))
	{
		// equ 1, result is pi/2
		if (n.Equ1())
		{
			if (arccos2)
			{
				if (neg)
					this->Set(this->Pi);
				else
					this->Set0();
			}
			else
			{
				this->Set(this->Pi12);
				if (neg) this->Neg();
			}
		}

		// overflow
		else
		{
			if (neg)
				this->SetMin();
			else
				this->SetMax();
		}
		return;
	}

	// square
	DBLFIX sqr;
	sqr.Sqr(n);

	// too many loops, convert to similar function y = sqrt(1 - x*x)
	bool sim = false;
	if (n.Greater(DBLFIX::Sqrt12)) // aprox. 0.7 (can be 0.5, too)
	{
		sqr.Sub((BASESINT)1, sqr);
		n.Sqrt(sqr);
		sim = true;
		arccos2 = !arccos2;
	}

	// prepare accumulator
	this->Set(n);

	// loop (up to 100 loops)
	BASESINT i;
	DBLFIX tmp;
	for (i = 2; ; i += 2)
	{
		// multiply increment by square
		n.Mul(sqr);

		// multiply increment by index-1
		n.Mul((BASESINT)(i-1));

		// divide increment by index
		n.Div(i);

		// divide increment by index+1
		tmp.Div(n, (BASESINT)(i+1));

		// increment = 0
		if (tmp.IsZero()) break;

		// add increment to accumulator
		this->Add(tmp);
	}

	// negate correction
	if (neg) this->Neg();

	// arcsin/arccos correcton
	if (arccos2) this->Sub(DBLFIX::Pi12, *this);

	// similar correction
	if (sim && neg)
	{
		if (arccos)
			this->Add(DBLFIX::Pi);
		else
			this->Sub(DBLFIX::Pi);
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = arctan(num), result in radians
//   ARCTG = x - x^3/3 + x^5/5 - x^7/7 + ....

void DBLFIX::ATan0(const DBLFIX& num, bool cotg)
{
	// get number
	DBLFIX n;
	n.Set(num);

	// absolute value
	bool neg = false;
	if (n.IsNeg())
	{
		n.Neg();
		neg = true;
	}

	// input is 1, result is PI/4
	if (n.Equ1())
	{
		this->Set(this->Pi14);
		if (neg) this->Neg();
		return;
	}

	// input > 1
	if (n.Greater((BASESINT)1))
	{
		cotg = !cotg;
		n.Recip();
	}

	// prepare accumulator
	this->Set(n);

	// -square
	DBLFIX sqr;
	sqr.Sqr(n);
	sqr.Neg();

	// loop
	BASESINT i;
	DBLFIX tmp;
	for (i = 3; ; i += 2)
	{
		// multiply increment by square
		n.Mul(sqr);

		// divide increment by index
		tmp.Div(n, i);

		// increment = 0
		if (tmp.IsZero()) break;

		// add increment to accumulator
		this->Add(tmp);
	}

	// cotg correction (Pi/2 - this)
	if (cotg) this->Sub(this->Pi12, *this);

	// negate correction
	if (neg) this->Neg();
}

///////////////////////////////////////////////////////////////////////////////
// this = log(num), natural logarithm
//    D = (x-1)/(x+1), x will be normalized into range 1..3
// LN x = 2*D + 2*D/3*D^2 + 2*D/5*D^4 + 2*D/7*D^6 + 2*D/9*D^8 + ...

void DBLFIX::Log(const DBLFIX& num)
{
	// zero or negative
	DBLFIX k;
	k.Set(num);
	if (num.IsZero() || num.IsNeg())
	{
		this->SetMin();
		return;
	}

	// 1
	if (num.Equ1())
	{
		this->Set0();
		return;
	}

	// normalize number to range 1..2, to speed up approximation
	int i = k.Bits() - BASEBITS;
	if (i < 0) k.LShift(-i);
	if (i > 0) k.RShift(i);

	// prepare coefficient
	DBLFIX d, m;
	d.Dec(k);
	k.Inc();
	d.Div(k);

	// prepare increment
	k.Set(d);
	k.Add(k);

	// prepare accumulator
	this->Set(k);

	// add normalization correction
	while (i <= -16) { this->Sub(this->Ln2_16); i += 16; }
	while (i <= -1) { this->Sub(this->Ln2); i++; }
	while (i >= 16) { this->Add(this->Ln2_16); i -= 16; }
	while (i >= 1) { this->Add(this->Ln2); i--; }

	// square of coefficient
	d.Sqr();

	// loop
	for (i = 3; ; i += 2)
	{
		// add coefficient
		k.Mul(d);

		// add counter
		m.Div(k, (BASESINT)i);

		// end of loop
		if (m.IsZero()) break;

		// add increment to accumulator
		this->Add(m);
	}
}

///////////////////////////////////////////////////////////////////////////////
// this = exp(num), natural exponent
//   e^x = 1 + x/1! + x^2/2! + x^3/3! + x^4/4! + x^5/5! + ...

void DBLFIX::Exp(const DBLFIX& num)
{
	// save input variable
	DBLFIX x;
	x.Set(num);

	// prepare increment
	DBLFIX k;
	k.Set(x);

	// prepare accumulator
	this->Inc(k);

	// loop
	BASESINT i;
	for (i = 2; ; i++)
	{
		// multiply increment
		k.Mul(x);
		k.Div(i);

		// check end of approximation
		if (k.IsZero() || ((DBLSINT*)&k)->EquM1()) break;

		// add increment to accumulator
		this->Add(k);
	}
}

#undef DBLUINT
#undef BASEUINT
#undef HALFUINT

#undef DBLSINT
#undef BASESINT
#undef HALFSINT

#undef DBLFIX

#undef USEASM
#undef SOFTASM

#undef BASEBITS
#undef BASELAST
#undef HALFBITS
#undef HALFMASK
#undef BASENOT

#undef DBLFIX_MAX
#undef DBLFIX_MIN
#undef DBLFIX_DEC
