
///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//                        Quadruple Integer (uint256)                        //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

//#define HWSHIFT		// uncomment this flag to ensure hardware compatibility of bit shift
						// (number of shifts is masked to 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 QUADUINT

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

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

#define QUADFIX fix256		// quadruple fixed point class

#define USEASM				// can use assembler

#endif // QUADUINT


#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 QUADFIX_MAX ((double)(BASELAST-1)*BASELAST*2 + (double)BASENOT + \
	(double)BASENOT/BASELAST/2 + (double)BASENOT/BASELAST/2/BASELAST/2) // max. decimal number as double
#define QUADFIX_MIN (-(double)BASELAST*BASELAST*2) // min. decimal number as double
#define QUADFIX_DEC ((int)(BASEBITS*2*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
QUADUINT* QUADUINT::Fact_;
QUADSINT* QUADSINT::Fact_;
QUADFIX* QUADFIX::Fact_;
QUADFIX* QUADFIX::FactRec_;

// initialize factorial tablbe
bool QUADUINT::FactOK = QUADUINT::FactInit();
bool QUADSINT::FactOK = QUADSINT::FactInit();
bool QUADFIX::FactOK = QUADFIX::FactInit();

// max. factorial
BASEUINT QUADUINT::FactMax = 0;
BASEUINT QUADSINT::FactMax = 0;
BASEUINT QUADFIX::FactMax = 0;

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

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

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

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

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

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

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

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

// Pi constant
QUADFIX QUADFIX::Pi;
CText QUADFIX::PiTxt;

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

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

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

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

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

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

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

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

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

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

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

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

// assembler
#ifdef USEASM
extern "C" void AddU256A_x64(QUADUINT*, const QUADUINT*);
extern "C" void AddU256B_x64(QUADUINT*, const DBLUINT*);
extern "C" void AddU256C_x64(QUADUINT*, BASEUINT);
extern "C" void AddU256D_x64(QUADUINT*, const QUADUINT*, const QUADUINT*);
extern "C" void AddU256E_x64(QUADUINT*, const QUADUINT*, const DBLUINT*);
extern "C" void AddU256F_x64(QUADUINT*, const QUADUINT*, BASEUINT);

extern "C" void SubU256A_x64(QUADUINT*, const QUADUINT*);
extern "C" void SubU256B_x64(QUADUINT*, const DBLUINT*);
extern "C" void SubU256C_x64(QUADUINT*, BASEUINT);
extern "C" void SubU256D_x64(QUADUINT*, const QUADUINT*, const QUADUINT*);
extern "C" void SubU256E_x64(QUADUINT*, const QUADUINT*, const DBLUINT*);
extern "C" void SubU256F_x64(QUADUINT*, const DBLUINT*, const QUADUINT*);
extern "C" void SubU256G_x64(QUADUINT*, const QUADUINT*, BASEUINT);
extern "C" void SubU256H_x64(QUADUINT*, BASEUINT, const QUADUINT*);
extern "C" void SubU256I_x64(QUADUINT*, const QUADUINT*);
extern "C" void SubU256J_x64(QUADUINT*, const DBLUINT*);
extern "C" void SubU256K_x64(QUADUINT*, BASEUINT);

extern "C" void IncU256A_x64(QUADUINT*);
extern "C" void IncU256B_x64(QUADUINT*, const QUADUINT*);
extern "C" void DecU256A_x64(QUADUINT*);
extern "C" void DecU256B_x64(QUADUINT*, const QUADUINT*);

extern "C" void MulU256A_x64(QUADUINT*, const QUADUINT*, const QUADUINT*);
extern "C" void MulU256B_x64(QUADUINT*, const QUADUINT*, const DBLUINT*);
extern "C" void MulU256C_x64(QUADUINT*, const QUADUINT*, BASEUINT);
extern "C" void MulU256D_x64(QUADUINT*, const DBLUINT*, const DBLUINT*);
extern "C" void MulU256E_x64(QUADUINT*, const DBLUINT*, BASEUINT);
extern "C" void MulU256F_x64(QUADUINT*, BASEUINT, BASEUINT);
extern "C" void MulU256G_x64(QUADUINT*, const QUADUINT*, const QUADUINT*, QUADUINT*);
extern "C" int MulU256H_x64(QUADUINT*, const QUADUINT*);

extern "C" void SqrU256A_x64(QUADUINT*, const QUADUINT*);
extern "C" void SqrU256B_x64(QUADUINT*, const DBLUINT*);
extern "C" void SqrU256C_x64(QUADUINT*, BASEUINT);

extern "C" void DivU256A_x64(QUADUINT*, const QUADUINT*, const QUADUINT*, QUADUINT*);
extern "C" void DivU256B_x64(QUADUINT*, const QUADUINT*, const DBLUINT*, QUADUINT*);
extern "C" BASEUINT DivU256C_x64(QUADUINT*, const QUADUINT*, BASEUINT);
extern "C" HALFUINT DivU256D_x64(QUADUINT*, const QUADUINT*, HALFUINT);
extern "C" void DivU256E_x64(QUADUINT*, const QUADUINT*, const QUADUINT*, const QUADUINT*, QUADUINT*);
extern "C" int DivU256F_x64(QUADUINT*, const QUADUINT*);

extern "C" void NegU256A_x64(QUADUINT*);
extern "C" void NegU256B_x64(QUADUINT*, const QUADUINT*);
extern "C" void NegU256C_x64(QUADUINT*, const DBLUINT*);
extern "C" void NegU256D_x64(QUADUINT*, BASEUINT);

extern "C" s64 BitsU256_x64(const QUADUINT*);
#endif // USEASM

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

void QUADUINT::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 QUADUINT::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
}

///////////////////////////////////////////////////////////////////////////////
// square base entry, r1:r0 = a * a

void QUADUINT::SqrB(BASEUINT* r0, BASEUINT* r1, BASEUINT a)
{
	//       a1a0
	//     x a1a0
	// ----------
	//       a0a0 ...   t0
	//     a0a1   ...  t1
	//     a1a0   ...  t2
	//   a1a1     ... t3

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

	// 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 
}

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

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

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

	// loop
	QUADUINT 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
		QUADUINT::Fact_ = (QUADUINT*)realloc(QUADUINT::Fact_, (size_t)((i+1)*sizeof(QUADUINT)));
		if (QUADUINT::Fact_ == NULL)
		{
			QUADUINT::FactMax = 0;
			return false;
		}

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

	return true;
}

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

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

void QUADUINT::Fact()
{
	if ((this->N1() == 0) && (this->N2() == 0) && (this->N3() == 0))
		this->Fact(this->N0());
	else
		this->SetMax();
}

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

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

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

	r = num.N2();
	num.N2() = this->N2();
	this->N2() = r;

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

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

void QUADUINT::Add(const QUADUINT& num)
{
#ifdef SOFTASM
	AddU256A_x64(this, &num);
#else // SOFTASM
	BASEUINT a = this->N0();
	BASEUINT b = num.N0();
	a += b;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a;

	a = this->N1() + c;
	c = (a < c) ? 1 : 0;
	b = num.N1();
	a += b;
	c += (a < b) ? 1 : 0;
	this->N1() = a;

	a = this->N2() + c;
	c = (a < c) ? 1 : 0;
	b = num.N2();
	a += b;
	c += (a < b) ? 1 : 0;
	this->N2() = a;

	this->N3() = this->N3() + num.N3() + c;
#endif // SOFTASM
}

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

void QUADUINT::Add(const DBLUINT& num)
{
#ifdef SOFTASM
	AddU256B_x64(this, &num);
#else // SOFTASM
	BASEUINT a = this->N0();
	BASEUINT b = num.N0();
	a += b;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a;

	a = this->N1() + c;
	c = (a < c) ? 1 : 0;
	b = num.N1();
	a += b;
	c += (a < b) ? 1 : 0;
	this->N1() = a;

	if (c)
	{
		a = this->N2() + c;
		c = (a < c) ? 1 : 0;
		this->N2() = a;

		if (c) this->N3() = this->N3() + c;
	}
#endif // SOFTASM
}

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

void QUADUINT::Add(BASEUINT num)
{
#ifdef SOFTASM
	AddU256C_x64(this, num);
#else // SOFTASM
	BASEUINT a = this->N0();
	BASEUINT b = num;
	a += b;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a;

	if (c)
	{
		a = this->N1() + c;
		c = (a < c) ? 1 : 0;
		this->N1() = a;

		if (c)
		{
			a = this->N2() + c;
			c = (a < c) ? 1 : 0;
			this->N2() = a;

			if (c) this->N3() = this->N3() + c;
		}
	}
#endif // SOFTASM
}

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

void QUADUINT::Add(const QUADUINT& num1, const QUADUINT& num2)
{
#ifdef SOFTASM
	AddU256D_x64(this, &num1, &num2);
#else // SOFTASM
	BASEUINT a = num1.N0();
	BASEUINT b = num2.N0();
	a += b;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a;

	a = num1.N1() + c;
	c = (a < c) ? 1 : 0;
	b = num2.N1();
	a += b;
	c += (a < b) ? 1 : 0;
	this->N1() = a;

	a = num1.N2() + c;
	c = (a < c) ? 1 : 0;
	b = num2.N2();
	a += b;
	c += (a < b) ? 1 : 0;
	this->N2() = a;

	this->N3() = num1.N3() + num2.N3() + c;
#endif // SOFTASM
}

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

void QUADUINT::Add(const QUADUINT& num1, const DBLUINT& num2)
{
#ifdef SOFTASM
	AddU256E_x64(this, &num1, &num2);
#else // SOFTASM
	BASEUINT a = num1.N0();
	BASEUINT b = num2.N0();
	a += b;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a;

	a = num1.N1() + c;
	c = (a < c) ? 1 : 0;
	b = num2.N1();
	a += b;
	c += (a < b) ? 1 : 0;
	this->N1() = a;

	a = num1.N2() + c;
	c = (a < c) ? 1 : 0;
	this->N2() = a;

	this->N3() = num1.N3() + c;
#endif // SOFTASM
}

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

void QUADUINT::Add(const QUADUINT& num1, BASEUINT num2)
{
#ifdef SOFTASM
	AddU256F_x64(this, &num1, num2);
#else // SOFTASM
	BASEUINT a = num1.N0();
	BASEUINT b = num2;
	a += b;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a;

	a = num1.N1() + c;
	c = (a < c) ? 1 : 0;
	this->N1() = a;

	a = num1.N2() + c;
	c = (a < c) ? 1 : 0;
	this->N2() = a;

	this->N3() = num1.N3() + c;
#endif // SOFTASM
}

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

void QUADUINT::Sub(const QUADUINT& num)
{
#ifdef SOFTASM
	SubU256A_x64(this, &num);
#else // SOFTASM
	BASEUINT b = num.N0();
	BASEUINT a = this->N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = num.N1() + c;
	c = (b < c) ? 1 : 0;
	a = this->N1();
	c += (a < b) ? 1 : 0;
	this->N1() = a - b;

	b = num.N2() + c;
	c = (b < c) ? 1 : 0;
	a = this->N2();
	c += (a < b) ? 1 : 0;
	this->N2() = a - b;

	this->N3() -= num.N3() + c;
#endif // SOFTASM
}

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

void QUADUINT::Sub(const DBLUINT& num)
{
#ifdef SOFTASM
	SubU256B_x64(this, &num);
#else // SOFTASM
	BASEUINT b = num.N0();
	BASEUINT a = this->N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = num.N1() + c;
	c = (b < c) ? 1 : 0;
	a = this->N1();
	c += (a < b) ? 1 : 0;
	this->N1() = a - b;

	if (c)
	{
		a = this->N2();
		c = (a < 1) ? 1 : 0;
		a--;
		this->N2() = a;

		if (c) this->N3()--;
	}
#endif // SOFTASM
}

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

void QUADUINT::Sub(BASEUINT num)
{
#ifdef SOFTASM
	SubU256C_x64(this, num);
#else // SOFTASM
	BASEUINT b = num;
	BASEUINT a = this->N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	if (c)
	{
		a = this->N1();
		c = (a < 1) ? 1 : 0;
		a--;
		this->N1() = a;

		if (c)
		{
			a = this->N2();
			c = (a < 1) ? 1 : 0;
			a--;
			this->N2() = a;

			if (c) this->N3()--;
		}
	}
#endif // SOFTASM
}

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

void QUADUINT::Sub(const QUADUINT& num1, const QUADUINT& num2)
{
#ifdef SOFTASM
	SubU256D_x64(this, &num1, &num2);
#else // SOFTASM
	BASEUINT b = num2.N0();
	BASEUINT a = num1.N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = num2.N1() + c;
	c = (b < c) ? 1 : 0;
	a = num1.N1();
	c += (a < b) ? 1 : 0;
	this->N1() = a - b;

	b = num2.N2() + c;
	c = (b < c) ? 1 : 0;
	a = num1.N2();
	c += (a < b) ? 1 : 0;
	this->N2() = a - b;

	this->N3() = num1.N3() - num2.N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::Sub(const QUADUINT& num1, const DBLUINT& num2)
{
#ifdef SOFTASM
	SubU256E_x64(this, &num1, &num2);
#else // SOFTASM
	BASEUINT b = num2.N0();
	BASEUINT a = num1.N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = num2.N1() + c;
	c = (b < c) ? 1 : 0;
	a = num1.N1();
	c += (a < b) ? 1 : 0;
	this->N1() = a - b;

	b = c;
	a = num1.N2();
	c = (a < b) ? 1 : 0;
	this->N2() = a - b;

	this->N3() = num1.N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::Sub(const DBLUINT& num1, const QUADUINT& num2)
{
#ifdef SOFTASM
	SubU256F_x64(this, &num1, &num2);
#else // SOFTASM
	BASEUINT b = num2.N0();
	BASEUINT a = num1.N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = num2.N1() + c;
	c = (b < c) ? 1 : 0;
	a = num1.N1();
	c += (a < b) ? 1 : 0;
	this->N1() = a - b;

	b = num2.N2() + c;
	c = (b < c) ? 1 : 0;
	c += (0 < b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - num2.N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::Sub(const QUADUINT& num1, BASEUINT num2)
{
#ifdef SOFTASM
	SubU256G_x64(this, &num1, num2);
#else // SOFTASM
	BASEUINT b = num2;
	BASEUINT a = num1.N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = c;
	a = num1.N1();
	c = (a < b) ? 1 : 0;
	this->N1() = a - b;

	b = c;
	a = num1.N2();
	c = (a < b) ? 1 : 0;
	this->N2() = a - b;

	this->N3() = num1.N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::Sub(BASEUINT num1, const QUADUINT& num2)
{
#ifdef SOFTASM
	SubU256H_x64(this, num1, &num2);
#else // SOFTASM
	BASEUINT b = num2.N0();
	BASEUINT a = num1;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = num2.N1() + c;
	c = (b < c) ? 1 : 0;
	c += (0 < b) ? 1 : 0;
	this->N1() = 0 - b;

	b = num2.N2() + c;
	c = (b < c) ? 1 : 0;
	c += (0 < b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - num2.N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::InvSub(const QUADUINT& num)
{
#ifdef SOFTASM
	SubU256I_x64(this, &num);
#else // SOFTASM
	BASEUINT b = this->N0();
	BASEUINT a = num.N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = this->N1() + c;
	c = (b < c) ? 1 : 0;
	a = num.N1();
	c += (a < b) ? 1 : 0;
	this->N1() = a - b;

	b = this->N2() + c;
	c = (b < c) ? 1 : 0;
	a = num.N2();
	c += (a < b) ? 1 : 0;
	this->N2() = a - b;

	this->N3() = num.N3() - this->N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::InvSub(const DBLUINT& num)
{
#ifdef SOFTASM
	SubU256J_x64(this, &num);
#else // SOFTASM
	BASEUINT b = this->N0();
	BASEUINT a = num.N0();
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = this->N1() + c;
	c = (b < c) ? 1 : 0;
	a = num.N1();
	c += (a < b) ? 1 : 0;
	this->N1() = a - b;

	b = this->N2() + c;
	c = (b < c) ? 1 : 0;
	c += (0 < b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - this->N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::InvSub(BASEUINT num)
{
#ifdef SOFTASM
	SubU256K_x64(this, num);
#else // SOFTASM
	BASEUINT b = this->N0();
	BASEUINT a = num;
	BASEUINT c = (a < b) ? 1 : 0;
	this->N0() = a - b;

	b = this->N1() + c;
	c = (b < c) ? 1 : 0;
	c += (0 < b) ? 1 : 0;
	this->N1() = 0 - b;

	b = this->N2() + c;
	c = (b < c) ? 1 : 0;
	c += (0 < b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - this->N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::Inc()
{
#ifdef SOFTASM
	IncU256A_x64(this);
#else // SOFTASM
	BASEUINT r = this->N0() + 1;
	BASEUINT c = (r < 1) ? 1 : 0;
	this->N0() = r;

	if (c)
	{
		r = this->N1() + c;
		c = (r < c) ? 1 : 0;
		this->N1() = r;

		if (c)
		{
			r = this->N2() + c;
			c = (r < c) ? 1 : 0;
			this->N2() = r;

			if (c) this->N3()++;
		}
	}
#endif // SOFTASM
}

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

void QUADUINT::Inc(const QUADUINT& num)
{
#ifdef SOFTASM
	IncU256B_x64(this, &num);
#else // SOFTASM
	BASEUINT r = num.N0() + 1;
	BASEUINT c = (r < 1) ? 1 : 0;
	this->N0() = r;

	r = num.N1() + c;
	c = (r < c) ? 1 : 0;
	this->N1() = r;

	r = num.N2() + c;
	c = (r < c) ? 1 : 0;
	this->N2() = r;

	this->N3() = num.N3() + c;
#endif // SOFTASM
}

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

void QUADUINT::Dec()
{
#ifdef SOFTASM
	DecU256A_x64(this);
#else // SOFTASM
	BASEUINT r = this->N0();
	BASEUINT c = (r < 1) ? 1 : 0;
	this->N0() = r - 1;

	if (c)
	{
		r = this->N1();
		c = (r < 1) ? 1 : 0;
		this->N1() = r - 1;

		if (c)
		{
			r = this->N2();
			c = (r < 1) ? 1 : 0;
			this->N2() = r - 1;

			if (c) this->N3()--;
		}
	}
#endif // SOFTASM
}

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

void QUADUINT::Dec(const QUADUINT& num)
{
#ifdef SOFTASM
	DecU256B_x64(this, &num);
#else // SOFTASM
	BASEUINT r = num.N0();
	BASEUINT c = (r < 1) ? 1 : 0;
	this->N0() = r - 1;

	r = num.N1();
	BASEUINT c2 = (r < c) ? 1 : 0;
	this->N1() = r - c;

	r = num.N2();
	c = (r < c2) ? 1 : 0;
	this->N2() = r - c2;

	this->N3() = num.N3() - c;
#endif // SOFTASM
}

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

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

	// small number
	if ((num1.N2() == 0) && (num1.N3() == 0))
	{
		this->Mul(num2, *(const DBLUINT*)&num1);
		return;
	}

	if ((num2.N2() == 0) && (num2.N3() == 0))
	{
		this->Mul(num1, *(const DBLUINT*)&num2);
		return;
	}

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

	// prepare result accumulator
	int i, j;
	BASEUINT r[4];
	BASEUINT* d;
	BASEUINT t, carry, r0, r1;
	HALFUINT b0, b1;

	// multiply b0 (store to result)
	d = &r[0]; // destination
	t = b[0];
	b0 = (HALFUINT)t;
	b1 = (HALFUINT)(t >> HALFBITS);
	carry = 0;

	if (b1 == 0)
	{
		for (j = 0; j < 4; j++) // a[j] loop
		{
			this->MulBH(&r0, &r1, a[j], b0);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			*d++ = r0;
		}
	}
	else
	{
		for (j = 0; j < 4; j++) // a[j] loop
		{
			this->MulBB(&r0, &r1, a[j], b0, b1);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			*d++ = r0;
		}
	}

	// multiply b1,... (add to result)
	for (i = 1; i < 4; i++) // b[i] loop
	{
		d = &r[i]; // destination
		t = b[i];
		b0 = (HALFUINT)t;
		b1 = (HALFUINT)(t >> HALFBITS);
		carry = 0;

		if (b1 == 0)
		{
			for (j = 0; j < 4 - i; j++) // a[j] loop
			{
				this->MulBH(&r0, &r1, a[j], b0);
				r0 += carry;
				carry = (r0 < carry) ? 1 : 0;
				carry += r1;
				t = *d;
				r0 += t;
				carry += (r0 < t) ? 1 : 0;
				*d++ = r0;
			}
		}
		else
		{
			for (j = 0; j < 4 - i; j++) // a[j] loop
			{
				this->MulBB(&r0, &r1, a[j], b0, b1);
				r0 += carry;
				carry = (r0 < carry) ? 1 : 0;
				carry += r1;
				t = *d;
				r0 += t;
				carry += (r0 < t) ? 1 : 0;
				*d++ = r0;
			}
		}
	}

	// save result
	this->N0() = r[0];
	this->N1() = r[1];
	this->N2() = r[2];
	this->N3() = r[3];
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * 0:num2

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

	// small number
	if ((num1.N2() == 0) && (num1.N3() == 0))
	{
		this->Mul(*(const DBLUINT*)&num1, num2);
		return;
	}

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

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

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

	// multiply b0 (store to result)
	d = &r[0]; // destination
	t = b[0];
	b0 = (HALFUINT)t;
	b1 = (HALFUINT)(t >> HALFBITS);
	carry = 0;

	if (b1 == 0)
	{
		for (j = 0; j < 4; j++) // a[j] loop
		{
			this->MulBH(&r0, &r1, a[j], b0);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			*d++ = r0;
		}
	}
	else
	{
		for (j = 0; j < 4; j++) // a[j] loop
		{
			this->MulBB(&r0, &r1, a[j], b0, b1);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			*d++ = r0;
		}
	}

	// multiply b1 (add to result)
	d = &r[1]; // destination
	t = b[1];
	b0 = (HALFUINT)t;
	b1 = (HALFUINT)(t >> HALFBITS);
	carry = 0;

	if (b1 == 0)
	{
		for (j = 0; j < 4 - 1; j++) // a[j] loop
		{
			this->MulBH(&r0, &r1, a[j], b0);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			t = *d;
			r0 += t;
			carry += (r0 < t) ? 1 : 0;
			*d++ = r0;
		}
	}
	else
	{
		for (j = 0; j < 4 - 1; j++) // a[j] loop
		{
			this->MulBB(&r0, &r1, a[j], b0, b1);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			t = *d;
			r0 += t;
			carry += (r0 < t) ? 1 : 0;
			*d++ = r0;
		}
	}

	// save result
	this->N0() = r[0];
	this->N1() = r[1];
	this->N2() = r[2];
	this->N3() = r[3];
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = num1 * 0:0:0:num2

void QUADUINT::Mul(const QUADUINT& num1, BASEUINT num2)
{
#ifdef USEASM
	MulU256C_x64(this, &num1, num2);
#else // USEASM

	// small number
	if ((num1.N2() == 0) && (num1.N3() == 0))
	{
		this->Mul(*(const DBLUINT*)&num1, num2);
		return;
	}

	// prepare operands (r = a * b)
	const BASEUINT* a = num1.m_N;

	// prepare result accumulator
	BASEUINT carry, r0, r1;
	HALFUINT b0, b1;

	// multiply b0 (store to result)
	b0 = (HALFUINT)num2;
	b1 = (HALFUINT)(num2 >> HALFBITS);

	if (b1 == 0)
	{
		this->MulBH(&this->N0(), &carry, a[0], b0);

		this->MulBH(&r0, &r1, a[1], b0);
		r0 += carry;
		carry = (r0 < carry) ? 1 : 0;
		carry += r1;
		this->N1() = r0;

		this->MulBH(&r0, &r1, a[2], b0);
		r0 += carry;
		carry = (r0 < carry) ? 1 : 0;
		carry += r1;
		this->N2() = r0;

		this->MulBH(&r0, &r1, a[3], b0);
		r0 += carry;
		this->N3() = r0;
	}
	else
	{
		this->MulBB(&this->N0(), &carry, a[0], b0, b1);

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

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

		this->MulBB(&r0, &r1, a[3], b0, b1);
		r0 += carry;
		this->N3() = r0;
	}
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num1 * 0:num2

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

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

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

	// 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];
	this->N2() = r0;
	this->N3() = carry;
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num1 * 0:0:0:num2

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

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

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

	// prepare result accumulator
	BASEUINT carry, r0, r1;
	HALFUINT b0, b1;

	// multiply b0 (store to result)
	b0 = (HALFUINT)num2;
	b1 = (HALFUINT)(num2 >> HALFBITS);

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

	// save result
	this->N1() = r0;
	this->N2() = carry;
	this->N3() = 0;
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:0:0:num1 * 0:0:0:num2

void QUADUINT::Mul(BASEUINT num1, BASEUINT num2)
{
#ifdef USEASM
	MulU256F_x64(this, num1, num2);
#else // USEASM
	this->MulBB(&this->N0(), &this->N1(), num1, (HALFUINT)num2, (HALFUINT)(num2 >> HALFBITS));
	this->N2() = 0;
	this->N3() = 0;
#endif // USEASM
}

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

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

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

	// prepare result accumulator
	int i, j;
	BASEUINT r[8];
	BASEUINT* d;
	BASEUINT t, carry, r0, r1;
	HALFUINT b0, b1;

	// multiply b0 (store to result)
	d = &r[0]; // destination
	t = b[0];
	b0 = (HALFUINT)t;
	b1 = (HALFUINT)(t >> HALFBITS);
	carry = 0;

	if (b1 == 0)
	{
		for (j = 0; j < 4; j++) // a[j] loop
		{
			this->MulBH(&r0, &r1, a[j], b0);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			*d++ = r0;
		}
	}
	else
	{
		for (j = 0; j < 4; j++) // a[j] loop
		{
			this->MulBB(&r0, &r1, a[j], b0, b1);
			r0 += carry;
			carry = (r0 < carry) ? 1 : 0;
			carry += r1;
			*d++ = r0;
		}
	}
	r[4] = carry;

	// multiply b1,... (add to result)
	for (i = 1; i < 4; i++) // b[i] loop
	{
		d = &r[i]; // destination
		t = b[i];
		b0 = (HALFUINT)t;
		b1 = (HALFUINT)(t >> HALFBITS);
		carry = 0;

		if (b1 == 0)
		{
			for (j = 0; j < 4; j++) // a[j] loop
			{
				this->MulBH(&r0, &r1, a[j], b0);
				r0 += carry;
				carry = (r0 < carry) ? 1 : 0;
				carry += r1;
				t = *d;
				r0 += t;
				carry += (r0 < t) ? 1 : 0;
				*d++ = r0;
			}
		}
		else
		{
			for (j = 0; j < 4; j++) // a[j] loop
			{
				this->MulBB(&r0, &r1, a[j], b0, b1);
				r0 += carry;
				carry = (r0 < carry) ? 1 : 0;
				carry += r1;
				t = *d;
				r0 += t;
				carry += (r0 < t) ? 1 : 0;
				*d++ = r0;
			}
		}
		*d = carry;
	}

	// save result
	this->N0() = r[0];
	this->N1() = r[1];
	this->N2() = r[2];
	this->N3() = r[3];

	high->N0() = r[4];
	high->N1() = r[5];
	high->N2() = r[6];
	high->N3() = r[7];
#endif // USEASM
}

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

int QUADUINT::Mul10(const QUADUINT& num)
{
#ifdef USEASM
	return MulU256H_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);

	BASEUINT r4 = a[4] * (BASEUINT)10;
	r4 += (r3 >> HALFBITS);
	BASEUINT r5 = a[5] * (BASEUINT)10;
	r5 += (r4 >> HALFBITS);
	this->N2() = (r4 & HALFMASK) | (r5 << HALFBITS);

	BASEUINT r6 = a[6] * (BASEUINT)10;
	r6 += (r5 >> HALFBITS);
	BASEUINT r7 = a[7] * (BASEUINT)10;
	r7 += (r6 >> HALFBITS);
	this->N3() = (r6 & HALFMASK) | (r7 << HALFBITS);

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

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

void QUADUINT::Sqr(const QUADUINT& num)
{
#ifdef USEASM
	SqrU256A_x64(this, &num);
#else // USEASM
	this->Mul(num, num);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:num * 0:num

void QUADUINT::Sqr(const DBLUINT& num)
{
#ifdef USEASM
	SqrU256B_x64(this, &num);
#else // USEASM
	this->Mul(num, num);
#endif // USEASM
}

///////////////////////////////////////////////////////////////////////////////
// this = 0:0:0:num * 0:0:0:num

void QUADUINT::Sqr(BASEUINT num)
{
#ifdef USEASM
	SqrU256C_x64(this, num);
#else // USEASM
	this->MulBB(&this->N0(), &this->N1(), num, (HALFUINT)num, (HALFUINT)(num >> HALFBITS));
	this->N2() = 0;
	this->N3() = 0;
#endif // USEASM
}

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

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

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

	// prepare dividend
	BASEUINT d[8];
	d[0] = num1.N0();
	d[1] = num1.N1();
	d[2] = num1.N2();
	d[3] = num1.N3();
	d[4] = 0;
	d[5] = 0;
	d[6] = 0;
	d[7] = 0;

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

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

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

	// divide
	int i;
	BASEUINT c, b;
	for (i = 4*BASEBITS; i > 0; i--)
	{
		// shift dividend
		d[7] = (d[7] << 1) | (d[6] >> (BASEBITS-1));
		d[6] = (d[6] << 1) | (d[5] >> (BASEBITS-1));
		d[5] = (d[5] << 1) | (d[4] >> (BASEBITS-1));
		d[4] = (d[4] << 1) | (d[3] >> (BASEBITS-1));
		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[7] > n[3]) || ((d[7] == n[3]) &&
			((d[6] > n[2]) || ((d[6] == n[2]) &&
			((d[5] > n[1]) || ((d[5] == n[1]) &&
			(d[4] >= n[0])))))))
		{
			// sub divisor from dividend
			c = (d[4] < n[0]) ? 1 : 0;
			d[4] -= n[0];

			b = n[1] + c;
			c = (b < c) ? 1 : 0;
			c += (d[5] < b) ? 1 : 0;
			d[5] -= b;

			b = n[2] + c;
			c = (b < c) ? 1 : 0;
			c += (d[6] < b) ? 1 : 0;
			d[6] -= b;

			d[7] -= n[3] + c;

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

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

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

	// save remainder
	if (rem != NULL)
	{
		rem->N0() = d[4];
		rem->N1() = d[5];
		rem->N2() = d[6];
		rem->N3() = d[7];
	}
#endif // USEASM
}

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

void QUADUINT::Div(const QUADUINT& num1, const DBLUINT& num2, QUADUINT* 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
	DivU256B_x64(this, &num1, &num2, rem);
#else // USEASM

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

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

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

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

	// divide
	int i;
	BASEUINT c;
	for (i = 4*BASEBITS; i > 0; i--)
	{
		// shift dividend
		c = d[5];
		d[5] = (d[5] << 1) | (d[4] >> (BASEBITS-1));
		d[4] = (d[4] << 1) | (d[3] >> (BASEBITS-1));
		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[5] > n[1]) || ((d[5] == n[1]) && (d[4] >= n[0])))
		{
			// sub divisor from dividend
			c = (d[4] < n[0]) ? 1 : 0;
			d[4] -= n[0];
			d[5] = d[5] - n[1] - c;

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

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

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

	// save remainder
	if (rem != NULL)
	{
		rem->N0() = d[4];
		rem->N1() = d[5];
		rem->N2() = 0;
		rem->N3() = 0;
	}
#endif // USEASM
}

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

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

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

	// divide N3
	BASEUINT d = num1.N3(); // 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->N3() = a;

	// divide N2
	d = num1.N2(); // 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->N2() = a;

	// divide N1
	d = num1.N1(); // dividend
	a = 0; // accumulator
	m = BASELAST; // mask

	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->N1() = a;

	// divide N0
	d = num1.N0(); // dividend
	a = 0; // accumulator
	m = BASELAST; // mask

	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:0:0:0:0:num2, returns remainder

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

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

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

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

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

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

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

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

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

	// r2
	k0 = (n & HALFMASK) | (k1 << HALFBITS);
	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 QUADUINT::Div10(const QUADUINT& num)
{
#ifdef USEASM
	return DivU256F_x64(this, &num);
#else // USEASM

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

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

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

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

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

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

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

	// r2
	k0 = (n & HALFMASK) | (k1 << HALFBITS);
	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 QUADUINT::DivHigh(const QUADUINT& num1L, const QUADUINT& num1H, const QUADUINT& num2, QUADUINT* rem /* = NULL */)
{
	if (num1H >= num2)
	{
		this->SetMax();
		if (rem != NULL) rem->Set0();
		return;
	}

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

	// prepare dividend
	BASEUINT d[8];
	d[0] = num1L.N0();
	d[1] = num1L.N1();
	d[2] = num1L.N2();
	d[3] = num1L.N3();
	d[4] = num1H.N0();
	d[5] = num1H.N1();
	d[6] = num1H.N2();
	d[7] = num1H.N3();

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

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

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

	// divide
	int i;
	BASEUINT c, b;
	for (i = 4*BASEBITS; i > 0; i--)
	{
		// shift dividend
		c = d[7];
		d[7] = (d[7] << 1) | (d[6] >> (BASEBITS-1));
		d[6] = (d[6] << 1) | (d[5] >> (BASEBITS-1));
		d[5] = (d[5] << 1) | (d[4] >> (BASEBITS-1));
		d[4] = (d[4] << 1) | (d[3] >> (BASEBITS-1));
		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[7] > n[3]) || ((d[7] == n[3]) &&
			((d[6] > n[2]) || ((d[6] == n[2]) &&
			((d[5] > n[1]) || ((d[5] == n[1]) &&
			(d[4] >= n[0])))))))
		{
			// sub divisor from dividend
			c = (d[4] < n[0]) ? 1 : 0;
			d[4] -= n[0];

			b = n[1] + c;
			c = (b < c) ? 1 : 0;
			c += (d[5] < b) ? 1 : 0;
			d[5] -= b;

			b = n[2] + c;
			c = (b < c) ? 1 : 0;
			c += (d[6] < b) ? 1 : 0;
			d[6] -= b;

			d[7] -= n[3] + c;

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

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

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

	// save remainder
	if (rem != NULL)
	{
		rem->N0() = d[4];
		rem->N1() = d[5];
		rem->N2() = d[6];
		rem->N3() = d[7];
	}

#endif // USEASM
}

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

int QUADUINT::PreFastDiv(const QUADUINT& 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*4-1)==1)
	{
		this->Set(div);
		return FASTDIV_HIGH;
	}

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

	// prepare dividend
	QUADUINT 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*4-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(4*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 QUADUINT::FastDiv(const QUADUINT& num, const QUADUINT& 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:
		{
			QUADUINT n;
			n.MulHigh(num, mul, this);
			this->RShift(flags >> FASTDIV_BIT);
		}
		return;

	// full precission
	//case FASTDIV_FULL:
	default:
		{
			// multiply this = num * mul
			QUADUINT 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(4*BASEBITS - shift);
		}
		return;
	}
}

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

void QUADUINT::Neg()
{
#ifdef SOFTASM
	NegU256A_x64(this);
#else // SOFTASM
	BASEUINT b = this->N0();
	BASEUINT c = (0 != b) ? 1 : 0;
	this->N0() = 0 - b;

	b = this->N1() + c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N1() = 0 - b;

	b = this->N2() + c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - this->N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::Neg(const QUADUINT& num)
{
#ifdef SOFTASM
	NegU256B_x64(this, &num);
#else // SOFTASM
	BASEUINT b = num.N0();
	BASEUINT c = (0 != b) ? 1 : 0;
	this->N0() = 0 - b;

	b = num.N1() + c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N1() = 0 - b;

	b = num.N2() + c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - num.N3() - c;
#endif // SOFTASM
}

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

void QUADUINT::Neg(const DBLUINT& num)
{
#ifdef SOFTASM
	NegU256C_x64(this, &num);
#else // SOFTASM
	BASEUINT b = num.N0();
	BASEUINT c = (0 != b) ? 1 : 0;
	this->N0() = 0 - b;

	b = num.N1() + c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N1() = 0 - b;

	b = c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - c;
#endif // SOFTASM
}

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

void QUADUINT::Neg(BASEUINT num)
{
#ifdef SOFTASM
	NegU256D_x64(this, num);
#else // SOFTASM
	BASEUINT b = num;
	BASEUINT c = (0 != b) ? 1 : 0;
	this->N0() = 0 - b;

	b = c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N1() = 0 - b;

	b = c;
	c = (b < c) ? 1 : 0;
	c += (0 != b) ? 1 : 0;
	this->N2() = 0 - b;

	this->N3() = 0 - c;
#endif // SOFTASM
}

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

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

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

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

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

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

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

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

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

int QUADUINT::Bits() const
{
#ifdef USEASM
	return (int)BitsU256_x64(this);
#else // USEASM
	if (this->N3() == 0)
		if (this->N2() == 0)
			if (this->N1() == 0)
				return ::Bits64(this->N0());
			else
				return ::Bits64(this->N1()) + BASEBITS;
		else
			return ::Bits64(this->N2()) + 2*BASEBITS;
	else
		return ::Bits64(this->N3()) + 3*BASEBITS;
#endif // USEASM
}

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

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

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

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

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

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

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

	// 1 or 2 highest words
	if (shift >= 2*BASEBITS)
	{
		// 1 highest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or highest word
			if (shift >= 4*BASEBITS)
				this->N3() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '<<' < BASEBITS
				this->N3() = this->N0() << (shift - 3*BASEBITS);

			this->N2() = 0;
			this->N1() = 0;
			this->N0() = 0;
		}
		else
		{
			// 2 higher words
			if (shift == 2*BASEBITS)
			{
				this->N3() = this->N1();
				this->N2() = this->N0();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' < BASEBITS
				this->N3() = (this->N1() << (shift - 2*BASEBITS)) | (this->N0() >> (3*BASEBITS - shift));
				this->N2() = this->N0() << (shift - 2*BASEBITS);
			}
			this->N1() = 0;
			this->N0() = 0;
		}
	}
	else
	{
		// 3 higher words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N3() = this->N2();
				this->N2() = this->N1();
				this->N1() = this->N0();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' BASEBITS
				this->N3() = (this->N2() << (shift - BASEBITS)) | (this->N1() >> (2*BASEBITS - shift));
				this->N2() = (this->N1() << (shift - BASEBITS)) | (this->N0() >> (2*BASEBITS - shift));
				this->N1() = this->N0() << (shift - BASEBITS);
			}
			this->N0() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
			this->N3() = (this->N3() << shift) | (this->N2() >> (BASEBITS - shift));
			this->N2() = (this->N2() << shift) | (this->N1() >> (BASEBITS - shift));
			this->N1() = (this->N1() << shift) | (this->N0() >> (BASEBITS - shift));
			this->N0() <<= shift;
		}
	}
}

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

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

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

void QUADUINT::LShift(const QUADUINT& num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (4*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;
	}

	// 1 or 2 highest words
	if (shift >= 2*BASEBITS)
	{
		// 1 highest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or highest word
			if (shift >= BASEBITS*4)
				this->N3() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '<<' < BASEBITS
				this->N3() = num.N0() << (shift - 3*BASEBITS);

			this->N2() = 0;
			this->N1() = 0;
			this->N0() = 0;
		}
		else
		{
			// 2 higher words
			if (shift == 2*BASEBITS)
			{
				this->N3() = num.N1();
				this->N2() = num.N0();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' < BASEBITS
				this->N3() = (num.N1() << (shift - 2*BASEBITS)) | (num.N0() >> (3*BASEBITS - shift));
				this->N2() = num.N0() << (shift - 2*BASEBITS);
			}
			this->N1() = 0;
			this->N0() = 0;
		}
	}
	else
	{
		// 3 higher words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N3() = num.N2();
				this->N2() = num.N1();
				this->N1() = num.N0();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' BASEBITS
				this->N3() = (num.N2() << (shift - BASEBITS)) | (num.N1() >> (2*BASEBITS - shift));
				this->N2() = (num.N1() << (shift - BASEBITS)) | (num.N0() >> (2*BASEBITS - shift));
				this->N1() = num.N0() << (shift - BASEBITS);
			}
			this->N0() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
			this->N3() = (num.N3() << shift) | (num.N2() >> (BASEBITS - shift));
			this->N2() = (num.N2() << shift) | (num.N1() >> (BASEBITS - shift));
			this->N1() = (num.N1() << shift) | (num.N0() >> (BASEBITS - shift));
			this->N0() = num.N0() << shift;
		}
	}
}

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

void QUADUINT::LShift1(BASEUINT num)
{
	this->N3() = 0;
	this->N2() = 0;
	this->N1() = num >> (BASEBITS - 1);
	this->N0() = num << 1;
}

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

void QUADUINT::LShift(BASEUINT num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (4*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;
	}

	// 1 or 2 highest words
	if (shift >= 2*BASEBITS)
	{
		// 1 highest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or highest word
			if (shift >= BASEBITS*4)
				this->N3() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '<<' < BASEBITS
				this->N3() = num << (shift - 3*BASEBITS);

			this->N2() = 0;
			this->N1() = 0;
			this->N0() = 0;
		}
		else
		{
			// 2 higher words
			if (shift == 2*BASEBITS)
			{
				this->N3() = 0;
				this->N2() = num;
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' < BASEBITS
				this->N3() = num >> (3*BASEBITS - shift);
				this->N2() = num << (shift - 2*BASEBITS);
			}
			this->N1() = 0;
			this->N0() = 0;
		}
	}
	else
	{
		// 3 higher words
		if (shift >= BASEBITS)
		{
			this->N3() = 0;

			if (shift == BASEBITS)
			{
				this->N2() = 0;
				this->N1() = num;
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' BASEBITS
				this->N2() = num >> (2*BASEBITS - shift);
				this->N1() = num << (shift - BASEBITS);
			}
			this->N0() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
			this->N3() = 0;
			this->N2() = 0;
			this->N1() = num >> (BASEBITS - shift);
			this->N0() = num << shift;
		}
	}
}

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

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

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

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

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

	// 1 or 2 smallest words
	if (shift >= 2*BASEBITS)
	{
		// 1 smallest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or smallest word
			if (shift >= BASEBITS*4)
				this->N0() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '>>' < BASEBITS
				this->N0() = this->N3() >> (shift - 3*BASEBITS);

			this->N1() = 0;
			this->N2() = 0;
			this->N3() = 0;
		}
		else
		{
			// 2 lower words
			if (shift == 2*BASEBITS)
			{
				this->N0() = this->N2();
				this->N1() = this->N3();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' < BASEBITS
				this->N0() = (this->N2() >> (shift - 2*BASEBITS)) | (this->N3() << (3*BASEBITS - shift));
				this->N1() = this->N3() >> (shift - 2*BASEBITS);
			}
			this->N2() = 0;
			this->N3() = 0;
		}
	}
	else
	{
		// 3 lower words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N0() = this->N1();
				this->N1() = this->N2();
				this->N2() = this->N3();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' BASEBITS
				this->N0() = (this->N1() >> (shift - BASEBITS)) | (this->N2() << (2*BASEBITS - shift));
				this->N1() = (this->N2() >> (shift - BASEBITS)) | (this->N3() << (2*BASEBITS - shift));
				this->N2() = this->N3() >> (shift - BASEBITS);
			}
			this->N3() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
			this->N0() = (this->N0() >> shift) | (this->N1() << (BASEBITS - shift));
			this->N1() = (this->N1() >> shift) | (this->N2() << (BASEBITS - shift));
			this->N2() = (this->N2() >> shift) | (this->N3() << (BASEBITS - shift));
			this->N3() >>= shift;
		}
	}
}

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

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

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

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

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

	// 1 or 2 smallest words
	if (shift >= 2*BASEBITS)
	{
		// 1 smallest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or smallest word
			if (shift >= BASEBITS*4)
				this->N0() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '>>' < BASEBITS
				this->N0() = num.N3() >> (shift - 3*BASEBITS);

			this->N1() = 0;
			this->N2() = 0;
			this->N3() = 0;
		}
		else
		{
			// 2 lower words
			if (shift == 2*BASEBITS)
			{
				this->N0() = num.N2();
				this->N1() = num.N3();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' < BASEBITS
				this->N0() = (num.N2() >> (shift - 2*BASEBITS)) | (num.N3() << (3*BASEBITS - shift));
				this->N1() = num.N3() >> (shift - 2*BASEBITS);
			}
			this->N2() = 0;
			this->N3() = 0;
		}
	}
	else
	{
		// 3 lower words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N0() = num.N1();
				this->N1() = num.N2();
				this->N2() = num.N3();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' BASEBITS
				this->N0() = (num.N1() >> (shift - BASEBITS)) | (num.N2() << (2*BASEBITS - shift));
				this->N1() = (num.N2() >> (shift - BASEBITS)) | (num.N3() << (2*BASEBITS - shift));
				this->N2() = num.N3() >> (shift - BASEBITS);
			}
			this->N3() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
			this->N0() = (num.N0() >> shift) | (num.N1() << (BASEBITS - shift));
			this->N1() = (num.N1() >> shift) | (num.N2() << (BASEBITS - shift));
			this->N2() = (num.N2() >> shift) | (num.N3() << (BASEBITS - shift));
			this->N3() = num.N3() >> shift;
		}
	}
}

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

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

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

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

	// shift to right 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;
	else
		// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
		this->N0() = num >> shift;

	this->N1() = 0;
	this->N2() = 0;
	this->N3() = 0;
}

///////////////////////////////////////////////////////////////////////////////
// convert to ASCIIZ text into buffer of size approx. 1 + 2.41*size in bytes
// (uint256: 78+1, 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 QUADUINT::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
	QUADUINT 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 QUADUINT::ToText(bool sign /* = false */) const
{
	CText txt;
	int len = (int)(BASEBITS/8*4*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 QUADUINT::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 QUADUINT::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 QUADUINT::Sqrt(const QUADUINT& num, QUADUINT* 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
	QUADUINT 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 QUADSINT::FactInit()
{
	// already initialized
	if (QUADSINT::Fact_ != NULL) return true;

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

	// loop
	QUADSINT 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
		QUADSINT::Fact_ = (QUADSINT*)realloc(QUADSINT::Fact_, (size_t)((i+1)*sizeof(QUADSINT)));
		if (QUADSINT::Fact_ == NULL)
		{
			QUADSINT::FactMax = 0;
			return false;
		}

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

	return true;
}

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

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

void QUADSINT::Fact()
{
	if ((this->N1() == 0) && (this->N2() == 0) && (this->N3() == 0))
		this->Fact(this->N0());
	else
		if (this->N3() >= 0)
			this->SetMax();
		else
			this->Set0();
}

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

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

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

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

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

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

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

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

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

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

void QUADSINT::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);
			QUADUINT::Mul(n1, n2);
		}
		else
		{
			QUADUINT::Mul(n1, *(const DBLUINT*)&num2);
			this->Neg();
		}
	}
	else
	{
		if (num2.IsNeg())
		{
			DBLUINT n2;
			n2.Neg(*(const DBLUINT*)&num2);
			QUADUINT::Mul(*(const DBLUINT*)&num1, n2);
			this->Neg();
		}
		else
		{
			QUADUINT::Mul(*(const DBLUINT*)&num1, *(const DBLUINT*)&num2);
		}
	}
}

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

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

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

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

void QUADSINT::Mul(BASESINT num1, BASESINT num2)
{
	if (num1 < 0)
	{
		if (num2 < 0)
		{
			QUADUINT::Mul((BASEUINT)-num1, (BASEUINT)-num2);
		}
		else
		{
			QUADUINT::Mul((BASEUINT)-num1, (BASEUINT)num2);
			this->Neg();
		}
	}
	else
	{
		if (num2 < 0)
		{
			QUADUINT::Mul((BASEUINT)num1, (BASEUINT)-num2);
			this->Neg();
		}
		else
		{
			QUADUINT::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 QUADSINT::MulHigh(const QUADSINT& num1, const QUADSINT& num2, QUADSINT* high)
{
	if (num1.IsNeg())
	{
		QUADUINT n1;
		n1.Neg(*(const QUADUINT*)&num1);

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

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

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

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

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

///////////////////////////////////////////////////////////////////////////////
// this = s:num * s:num

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	// 1 or 2 highest words
	if (shift >= 2*BASEBITS)
	{
		// 1 highest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or highest word
			if (shift >= 4*BASEBITS)
				this->N3() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '<<' < BASEBITS
				this->N3() = this->N0() << (shift - 3*BASEBITS);

			this->N2() = 0;
			this->N1() = 0;
			this->N0() = 0;
		}
		else
		{
			// 2 higher words
			if (shift == 2*BASEBITS)
			{
				this->N3() = this->N1();
				this->N2() = this->N0();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' < BASEBITS
				this->N3() = (this->N1() << (shift - 2*BASEBITS)) | (this->N0() >> (3*BASEBITS - shift));
				this->N2() = this->N0() << (shift - 2*BASEBITS);
			}
			this->N1() = 0;
			this->N0() = 0;
		}
	}
	else
	{
		// 3 higher words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N3() = this->N2();
				this->N2() = this->N1();
				this->N1() = this->N0();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' BASEBITS
				this->N3() = (this->N2() << (shift - BASEBITS)) | (this->N1() >> (2*BASEBITS - shift));
				this->N2() = (this->N1() << (shift - BASEBITS)) | (this->N0() >> (2*BASEBITS - shift));
				this->N1() = this->N0() << (shift - BASEBITS);
			}
			this->N0() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
			this->N3() = (this->N3() << shift) | (this->N2() >> (BASEBITS - shift));
			this->N2() = (this->N2() << shift) | (this->N1() >> (BASEBITS - shift));
			this->N1() = (this->N1() << shift) | (this->N0() >> (BASEBITS - shift));
			this->N0() <<= shift;
		}
	}
}

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

void QUADSINT::LShift(const QUADSINT& num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (4*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;
	}

	// 1 or 2 highest words
	if (shift >= 2*BASEBITS)
	{
		// 1 highest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or highest word
			if (shift >= BASEBITS*4)
				this->N3() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '<<' < BASEBITS
				this->N3() = num.N0() << (shift - 3*BASEBITS);

			this->N2() = 0;
			this->N1() = 0;
			this->N0() = 0;
		}
		else
		{
			// 2 higher words
			if (shift == 2*BASEBITS)
			{
				this->N3() = num.N1();
				this->N2() = num.N0();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' < BASEBITS
				this->N3() = (num.N1() << (shift - 2*BASEBITS)) | (num.N0() >> (3*BASEBITS - shift));
				this->N2() = num.N0() << (shift - 2*BASEBITS);
			}
			this->N1() = 0;
			this->N0() = 0;
		}
	}
	else
	{
		// 3 higher words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N3() = num.N2();
				this->N2() = num.N1();
				this->N1() = num.N0();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' BASEBITS
				this->N3() = (num.N2() << (shift - BASEBITS)) | (num.N1() >> (2*BASEBITS - shift));
				this->N2() = (num.N1() << (shift - BASEBITS)) | (num.N0() >> (2*BASEBITS - shift));
				this->N1() = num.N0() << (shift - BASEBITS);
			}
			this->N0() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
			this->N3() = (num.N3() << shift) | (num.N2() >> (BASEBITS - shift));
			this->N2() = (num.N2() << shift) | (num.N1() >> (BASEBITS - shift));
			this->N1() = (num.N1() << shift) | (num.N0() >> (BASEBITS - shift));
			this->N0() = num.N0() << shift;
		}
	}
}

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

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

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

void QUADSINT::LShift(BASESINT num, int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (4*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;
	}

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

	// 1 or 2 highest words
	if (shift >= 2*BASEBITS)
	{
		// 1 highest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or highest word
			if (shift >= BASEBITS*4)
				this->N3() = 0;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '<<' < BASEBITS
				this->N3() = num << (shift - 3*BASEBITS);

			this->N2() = 0;
			this->N1() = 0;
			this->N0() = 0;
		}
		else
		{
			// 2 higher words
			if (shift == 2*BASEBITS)
			{
				this->N3() = k;
				this->N2() = num;
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' < BASEBITS
				this->N3() = ((BASEUINT)num >> (3*BASEBITS - shift)) | (k << (shift - 2*BASEBITS));
				this->N2() = num << (shift - 2*BASEBITS);
			}
			this->N1() = 0;
			this->N0() = 0;
		}
	}
	else
	{
		// 3 higher words
		if (shift >= BASEBITS)
		{
			this->N3() = k;

			if (shift == BASEBITS)
			{
				this->N2() = k;
				this->N1() = num;
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '<<' < BASEBITS, 0 < '>>' BASEBITS
				this->N2() = ((BASEUINT)num >> (2*BASEBITS - shift)) | (k << (shift - BASEBITS));
				this->N1() = num << (shift - BASEBITS);
			}
			this->N0() = 0;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '>>' < BASEBITS)
			this->N3() = k;
			this->N2() = k;
			this->N1() = ((BASEUINT)num >> (BASEBITS - shift)) | (k << shift);
			this->N0() = num << shift;
		}
	}
}

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

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

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

void QUADSINT::RShift(int shift)
{
#ifdef HWSHIFT
	// hardware compatibility (mask to number of bits)
	shift = (unsigned int)shift & (4*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;

	// 1 or 2 smallest words
	if (shift >= 2*BASEBITS)
	{
		// 1 smallest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or smallest word
			if (shift >= BASEBITS*4)
				this->N0() = k;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '>>' < BASEBITS
				if (shift == BASEBITS*3)
					this->N0() = (BASEUINT)this->N3();
				else
					this->N0() = ((BASEUINT)this->N3() >> (shift - 3*BASEBITS)) | (k << (4*BASEBITS - shift));

			this->N1() = k;
			this->N2() = k;
			this->N3() = k;
		}
		else
		{
			// 2 lower words
			if (shift == 2*BASEBITS)
			{
				this->N0() = this->N2();
				this->N1() = this->N3();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' < BASEBITS
				this->N0() = (this->N2() >> (shift - 2*BASEBITS)) | (this->N3() << (3*BASEBITS - shift));
				this->N1() = ((BASEUINT)this->N3() >> (shift - 2*BASEBITS)) | (k << (3*BASEBITS - shift));
			}
			this->N2() = k;
			this->N3() = k;
		}
	}
	else
	{
		// 3 lower words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N0() = this->N1();
				this->N1() = this->N2();
				this->N2() = this->N3();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' BASEBITS
				this->N0() = (this->N1() >> (shift - BASEBITS)) | (this->N2() << (2*BASEBITS - shift));
				this->N1() = (this->N2() >> (shift - BASEBITS)) | (this->N3() << (2*BASEBITS - shift));
				this->N2() = ((BASEUINT)this->N3() >> (shift - BASEBITS)) | (k << (2*BASEBITS - shift)) ;
			}
			this->N3() = k;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
			this->N0() = (this->N0() >> shift) | (this->N1() << (BASEBITS - shift));
			this->N1() = (this->N1() >> shift) | (this->N2() << (BASEBITS - shift));
			this->N2() = (this->N2() >> shift) | (this->N3() << (BASEBITS - shift));
			this->N3() = ((BASEUINT)this->N3() >> shift) | (k << (BASEBITS - shift));
		}
	}
}

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

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

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

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

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

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

	// 1 or 2 smallest words
	if (shift >= 2*BASEBITS)
	{
		// 1 smallest word or overload
		if (shift >= 3*BASEBITS)
		{
			// overload or smallest word
			if (shift >= BASEBITS*4)
				this->N0() = k;
			else
				// here: 3*BASEBITS <= shift < 4*BASEBITS
				//       0 <= '>>' < BASEBITS
				if (shift == 3*BASEBITS)
					this->N0() = (BASEUINT)num.N3();
				else
					this->N0() = ((BASEUINT)num.N3() >> (shift - 3*BASEBITS)) | (k << (4*BASEBITS - shift));

			this->N1() = k;
			this->N2() = k;
			this->N3() = k;
		}
		else
		{
			// 2 lower words
			if (shift == 2*BASEBITS)
			{
				this->N0() = num.N2();
				this->N1() = num.N3();
			}
			else
			{
				// here: 2*BASEBITS < shift < 3*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' < BASEBITS
				this->N0() = (num.N2() >> (shift - 2*BASEBITS)) | (num.N3() << (3*BASEBITS - shift));
				this->N1() = ((BASEUINT)num.N3() >> (shift - 2*BASEBITS)) | (k << (3*BASEBITS - shift));
			}
			this->N2() = k;
			this->N3() = k;
		}
	}
	else
	{
		// 3 lower words
		if (shift >= BASEBITS)
		{
			if (shift == BASEBITS)
			{
				this->N0() = num.N1();
				this->N1() = num.N2();
				this->N2() = num.N3();
			}
			else
			{
				// here: BASEBITS < shift < 2*BASEBITS
				//       0 < '>>' < BASEBITS, 0 < '<<' BASEBITS
				this->N0() = (num.N1() >> (shift - BASEBITS)) | (num.N2() << (2*BASEBITS - shift));
				this->N1() = (num.N2() >> (shift - BASEBITS)) | (num.N3() << (2*BASEBITS - shift));
				this->N2() = ((BASEUINT)num.N3() >> (shift - BASEBITS)) | (k << (2*BASEBITS - shift));
			}
			this->N3() = k;
		}
		else
		{
			// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
			this->N0() = (num.N0() >> shift) | (num.N1() << (BASEBITS - shift));
			this->N1() = (num.N1() >> shift) | (num.N2() << (BASEBITS - shift));
			this->N2() = (num.N2() >> shift) | (num.N3() << (BASEBITS - shift));
			this->N3() = ((BASEUINT)num.N3() >> shift) | (k << (BASEBITS - shift));
		}
	}
}

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

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

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

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

	// shift to right 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;
	else
		// small shift (here: 0 < shift < BASEBITS, 0 < '<<' < BASEBITS)
		this->N0() = ((BASEUINT)num >> shift) | (k << (BASEBITS - shift));

	this->N1() = k;
	this->N2() = k;
	this->N3() = k;
}

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

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

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

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

	// loop
	QUADFIX 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
		QUADFIX::Fact_ = (QUADFIX*)realloc(QUADFIX::Fact_, (size_t)((i+1)*sizeof(QUADFIX)));
		QUADFIX::FactRec_ = (QUADFIX*)realloc(QUADFIX::FactRec_, (size_t)((i+1)*sizeof(QUADFIX)));
		if ((QUADFIX::Fact_ == NULL) || (QUADFIX::FactRec_ == NULL))
		{
			QUADFIX::FactMax = 0;
			return false;
		}

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

	return true;
}

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

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

void QUADFIX::Fact()
{
	if (this->N3() == 0)
		this->Fact(this->N2());
	else
		if (this->N3() >= 0)
			this->SetMax();
		else
			this->Set0();
}

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

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

void QUADFIX::FactRec()
{
	if (this->N3() == 0)
		this->Fact(this->N2());
	else
		if (this->N3() >= 0)
			this->Set0();
		else
			this->SetMax();
}

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

double QUADFIX::Double()
{
	return	(double)this->N3()*BASELAST*2 +
			(double)this->N2() +
			(double)this->N1()/BASELAST/2 +
			(double)this->N0()/BASELAST/2/BASELAST/2;
}

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

float QUADFIX::Float()
{
	return	(float)this->N3()*BASELAST*2 +
			(float)this->N2() +
			(float)this->N1()/BASELAST/2 +
			(float)this->N0()/BASELAST/2/BASELAST/2;
}

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

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

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

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

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

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

	// N3
	num *= BASELAST;
	num *= 2;
	k = (BASEUINT)num;
	num -= k;
	k2 = (BASEUINT)num;
	num -= k2;
	this->N1() = k + k2;

	// N4
	num *= BASELAST;
	num *= 2;
	k = (BASEUINT)num;
	num -= k;
	k2 = (BASEUINT)num;
	this->N0() = k + k2;

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

///////////////////////////////////////////////////////////////////////////////
// this = double, with modulo overflow

void QUADFIX::SetMod(double num)
{
	num = fmod(num + (double)BASELAST*2, (double)BASELAST*2) - (double)BASELAST*2;
	this->Set(num);
}

///////////////////////////////////////////////////////////////////////////////
// this = dblfix

void QUADFIX::Set(const DBLFIX& num)
{
	this->N0() = 0;
	this->N1() = num.N0();
	this->N2() = num.N1();
	this->N3() = (num.IsNeg() ? -1 : 0);
}

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

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

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

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

void QUADFIX::Div(const QUADFIX& num1, const QUADFIX& num2)
{
	QUADSINT n1, n2;
	BASESINT k = num1.IsNeg() ? -1 : 0;
	n2.N3() = k;
	n2.N2() = k;
	n2.N1() = num1.N3();
	n2.N0() = num1.N2();
	n1.N3() = num1.N1();
	n1.N2() = num1.N0();
	n1.N1() = 0;
	n1.N0() = 0;
	QUADSINT::DivHigh(n1, n2, num2);
}

///////////////////////////////////////////////////////////////////////////////
// convert to ASCIIZ text into buffer of size approx. 3 + 2.41*size in bytes
// (fix256: 1+38+1+39+1=80, fix128: 1+19+1+20+1=43, fix64: 1+9+1+10+1=23, returns number of characters
// dec = max. number of decimal places (max. QUADFIX_DEC+1)
//   default precision: fix256: 38, fix128: 19, fix64: 9
//   max. precision: fix256: 39, fix128: 20, fix64: 10 (last digit can be incorrect, previous one can be incorrectly rounded)

int QUADFIX::ToTextBuf(char* buf, int bufsize, int dec /* = QUADFIX_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 > QUADFIX_DEC+1) dec = QUADFIX_DEC+1;

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

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

	// prepare integer part
	DBLUINT s2;
	s2.N0() = s.N2();
	s2.N1() = s.N3();

	// decode integer part of the number (to end of buffer)
	bufsize--; // reserve space for terminating zero
	int n = 0; // counter of digits
	char* d = &buf[bufsize];
	while (n < bufsize)
	{
		d--;
		*d = s2.Div10() + '0';
		n++;
		if (s2.Equ0()) 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
	s.N3() = 0;
	while (n3 < dec)
	{
		s.N2() = 0;
		s.Mul10();
		*d++ = (char)(s.N2() + '0');
		n3++;
	}

	// round number up
	if ((BASESINT)s.N1() < 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 QUADFIX::ToText(int dec /* = QUADFIX_DEC */) const
{
	CText txt;
	int len = (int)(BASEBITS/8*4*2.5+3);
	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 QUADFIX::FromTextBuf(const char* text, int len /* = -1 */)
{
	// clear accumulator
	this->Set0();
	char ch;
	u32 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
	DBLUINT s;
	s.Set0();
	while ((len > 0) && ((ch = *text++) != 0))
	{
		if ((ch < '0') || (ch > '9')) break;
		len--;
		s.Mul10();
		s.Add((BASEUINT)(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->N2() = ch - '0';
			this->Div10();
		}
	}

	// store integer part
	this->N2() = s.N0();
	this->N3() = s.N1();

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

	return len0 - len;
}

// convert number from text (returns new position)
int QUADFIX::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 QUADFIX::ConstInit()
{
	BASESINT i;
	QUADFIX x, y, z;

	// Euler's constant e = 1 + 1/1! + 1/2! + 1/3! + 1/4! + ...
	// e = 2.718281828 4590452353 6028747135266249775 724709369995957496696762772407663035354759...
	// https://www.math.utah.edu/~pa/math/e.html
#define SHFT (BASEBITS/2) // higher precision
	x = (BASESINT)((BASESINT)2 << SHFT);
	y = (BASESINT)((BASESINT)1 << SHFT);
	i = 1;
	for (;;) // loops fix64: 15, fix128: 25, fix256: 41
	{
		i++;
		y.Div(i);
		if (y.IsZero()) break;
		x.Add(y);
	}

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

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

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

	QUADFIX::EulTxt = QUADFIX::Eul.ToText(); // decimal places fix64: 9, fix128: 19, fix256: 38
	// fix64: 2.718281828
	// fix128: 2.7182818284590452353
	// fix256: 2.71828182845904523536028747135266249776
	QUADFIX::Eul2Txt = QUADFIX::Eul2.ToText();
	QUADFIX::Eul4Txt = QUADFIX::Eul4.ToText();
	QUADFIX::Eul8Txt = QUADFIX::Eul8.ToText();
	QUADFIX::Eul1Txt = QUADFIX::Eul1.ToText();
	QUADFIX::Eul12Txt = QUADFIX::Eul12.ToText();
	QUADFIX::Eul14Txt = QUADFIX::Eul14.ToText();
	QUADFIX::Eul18Txt = QUADFIX::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.141592653 5897932384 6264338327950288419 716939937510582097494459230781640628620899...
	// 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 fix64: 35, fix128: 73, fix256: 153
	{
		// 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
	QUADFIX::Pi180.Div(x, (BASESINT)180);
	QUADFIX::Pi180.RShift(SHFT);
	QUADFIX::Pi180Txt = QUADFIX::Pi180.ToText();

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

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

	// Pi
	QUADFIX::Pi.RShift(x, SHFT);
	QUADFIX::PiTxt = QUADFIX::Pi.ToText(); // decimal places fix64: 9, fix128: 19, fix256: 37 (38)
	// fix64: 3.141592653
	// fix128: 3.1415926535897932385
	// fix256: 3.1415926535897932384626433832795028842(0)

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

	// Pi/4
	QUADFIX::Pi14.RShift(x, SHFT+2);
	QUADFIX::Pi14Txt = QUADFIX::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";
	QUADFIX::Ln10.FromText(t);
	QUADFIX::Ln10Txt = QUADFIX::Ln10.ToText();

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

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

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

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

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

	return true;
}

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

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

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

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

	// find root (usually around 6 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 QUADFIX::NormRad(const QUADFIX& num)
{
	// get number
	this->Set(num);

	// small overflow (-10*Pi .. 10*Pi)
	if (((this->N3() == -1) || ((this->N3() == 0) && (this->N2() < 3*11))) &&
		((this->N3() == 0) || ((this->N3() == -1) && (this->N2() > -3*11))))
	{
		while (this->IsNeg()) this->Add(QUADFIX::Pi2);
		while (this->GreaterEqu(QUADFIX::Pi2)) this->Sub(QUADFIX::Pi2);
	}
	else
	{
		this->Mod(QUADFIX::Pi2);
		if (this->IsNeg()) this->Add(QUADFIX::Pi2);
	}
}

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

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

	// small overflow (-10*360 .. 10*360)
	if (((this->N3() == -1) || ((this->N3() == 0) && (this->N2() < 360*10))) &&
		((this->N3() == 0) || ((this->N3() == -1) && (this->N2() > -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 QUADFIX::Sin(const QUADFIX& num)
{
	// normalize angle
	QUADFIX n;
	n.NormRad(num);

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

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

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

	// loop (1 to 12 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 QUADFIX::ASin0(const QUADFIX& num, bool arccos)
{
	// get number
	QUADFIX 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
	QUADFIX sqr;
	sqr.Sqr(n);

	// too many loops, convert to similar function y = sqrt(1 - x*x)
	bool sim = false;
	if (n.Greater(QUADFIX::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 220 loops)
	BASESINT i;
	QUADFIX 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(QUADFIX::Pi12, *this);

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

	i=i;
}

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

void QUADFIX::ATan0(const QUADFIX& num, bool cotg)
{
	// get number
	QUADFIX 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
	QUADFIX sqr;
	sqr.Sqr(n);
	sqr.Neg();

	// loop
	BASESINT i;
	QUADFIX 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 QUADFIX::Log(const QUADFIX& num)
{
	// zero or negative
	QUADFIX 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() - 2*BASEBITS;
	if (i < 0) k.LShift(-i);
	if (i > 0) k.RShift(i);

	// prepare coefficient
	QUADFIX 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 QUADFIX::Exp(const QUADFIX& num)
{
	// save input variable
	QUADFIX x;
	x.Set(num);

	// prepare increment
	QUADFIX 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() || ((QUADSINT*)&k)->EquM1()) break;

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

#undef QUADUINT
#undef DBLUINT
#undef BASEUINT
#undef HALFUINT

#undef QUADSINT
#undef DBLSINT
#undef BASESINT
#undef HALFSINT

#undef QUADFIX
#undef DBLFIX

#undef USEASM
#undef SOFTASM

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

#undef QUADFIX_MAX
#undef QUADFIX_MIN
#undef QUADFIX_DEC
