/****************************************************************************
* NFG - Numbers for Geometry                     					        *
*                                                                           *
* Consiglio Nazionale delle Ricerche                                        *
* Istituto di Matematica Applicata e Tecnologie Informatiche                *
* Sezione di Genova                                                         *
* IMATI-GE / CNR                                                            *
*                                                                           *
* Authors: Marco Attene                                                     *
* Copyright(C) 2019: IMATI-GE / CNR                                         *
* All rights reserved.                                                      *
*                                                                           *
* This program is free software; you can redistribute it and/or modify      *
* it under the terms of the GNU Lesser General Public License as published  *
* by the Free Software Foundation; either version 3 of the License, or (at  *
* your option) any later version.                                           *
*                                                                           *
* This program is distributed in the hope that it will be useful, but       *
* WITHOUT ANY WARRANTY; without even the implied warranty of                *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser  *
* General Public License for more details.                                  *
*                                                                           *
* You should have received a copy of the GNU Lesser General Public License  *
* along with this program.  If not, see http://www.gnu.org/licenses/.       *
*                                                                           *
****************************************************************************/

///////////////////////////////////////////////////////////////////////////////////////////////////////
//
// To compile on MSVC: use /fp:strict
// On GNU GCC: use -frounding-math
//
///////////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef NUMERICS_H
#define NUMERICS_H

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <float.h>
#include <math.h>
#include <fenv.h>
#include <iostream>
#include <climits>
#include "memPool.h"

#ifdef USE_GNU_GMP_CLASSES
#include <gmpxx.h>
#endif

// Call the following function (once per thread) before using these number types
void initFPU();

inline void ip_error(const char* msg)
{
	std::cerr << msg;
	exit(0);
}

#if INTPTR_MAX == INT64_MAX
#define	IS64BITPLATFORM
#endif

#ifdef _MSC_VER
#define	ISVISUALSTUDIO
#endif

#ifdef IS64BITPLATFORM
#ifdef __SSE2__
#define USE_SIMD_INSTRUCTIONS
#endif
#ifdef __AVX2__
#define USE_SIMD_INSTRUCTIONS
#define USE_AVX2_INSTRUCTIONS
#endif
#endif

#ifdef ISVISUALSTUDIO

#pragma fenv_access (on)

inline void setFPUModeToRoundUP() { _controlfp(_RC_UP, _MCW_RC); }
inline void setFPUModeToRoundNEAR() { _controlfp(_RC_NEAR, _MCW_RC); }

#else

#ifdef __AVX2__
#pragma GCC target("fma")
#endif

#pragma STDC FENV_ACCESS ON

inline void setFPUModeToRoundUP() { fesetround(FE_UPWARD); }
inline void setFPUModeToRoundNEAR() { fesetround(FE_TONEAREST); }
#endif

#ifdef USE_SIMD_INSTRUCTIONS

#ifdef USE_AVX2_INSTRUCTIONS
#include <immintrin.h>
#else
#include <emmintrin.h>
#endif

#endif

	/////////////////////////////////////////////////////////////////////
	// 	   
	// 	   I N T E R V A L   A R I T H M E T I C
	// 
	/////////////////////////////////////////////////////////////////////

	// An interval_number is a pair of doubles representing an interval.
	// Operations on interval_number require that the rounding mode is
	// set to +INFINITY. Use setFPUModeToRoundUP().

	class interval_number
	{
#ifdef USE_SIMD_INSTRUCTIONS
		__m128d interval; // interval[1] = min_low, interval[0] = high

		static inline __m128d zero() { return _mm_setzero_pd(); }
		static inline __m128d minus_one() { return _mm_set1_pd(-1.0); }
		static inline __m128d sign_low_mask() { return _mm_castsi128_pd(_mm_set_epi64x(LLONG_MIN, 0)); }
		static inline __m128d sign_high_mask() { return _mm_castsi128_pd(_mm_set_epi64x(0, LLONG_MIN)); }
		static inline __m128d sign_fabs_mask() { return _mm_castsi128_pd(_mm_set_epi64x(~LLONG_MIN, ~LLONG_MIN)); }
		static inline __m128d all_high_mask() { return _mm_castsi128_pd(_mm_set_epi64x(0, -1LL)); }

		__m128d getLowSwitched() const { return _mm_xor_pd(interval, sign_low_mask()); }

	public:
		const double *getInterval() const { return (const double*)&interval; }

		interval_number() { }
		interval_number(const double a) : interval(_mm_set_pd(-a, a)) {}
		interval_number(const double minf, const double sup) : interval(_mm_set_pd(minf, sup)) {}
		interval_number(const __m128d& i) : interval(i) {}
		interval_number(const interval_number& b) : interval(b.interval) {}

		double minus_inf() const { return _mm_cvtsd_f64(_mm_shuffle_pd(interval, interval, 1));	}
		double inf() const { return -minus_inf(); }
		double sup() const { return _mm_cvtsd_f64(interval); }

		interval_number& operator=(const interval_number& b) { interval = b.interval; return *this; }

		interval_number operator+(const interval_number& b) const { return interval_number(_mm_add_pd(interval, b.interval)); }

		interval_number operator-(const interval_number& b) const { return interval_number(_mm_add_pd(interval, _mm_shuffle_pd(b.interval, b.interval, 1))); }

		interval_number operator*(const interval_number& b) const;

		interval_number operator-() const { return interval_number(_mm_shuffle_pd(interval, interval, 1)); }
		interval_number operator+(const double b) const { return interval_number(_mm_add_pd(interval, _mm_set_pd(-b, b))); }
		interval_number operator-(const double b) const { return interval_number(_mm_sub_pd(interval, _mm_set_pd(-b, b))); }
		interval_number operator*(const double b) const;
		interval_number operator/(const double b) const;
		interval_number& operator+=(const interval_number& b) { return operator=(*this + b); }
		interval_number& operator-=(const interval_number& b) { return operator=(*this - b); }
		interval_number& operator*=(const interval_number& b) { return operator=(*this * b); }
		interval_number& operator+=(const double b) { return operator=(*this + b); }
		interval_number& operator-=(const double b) { return operator=(*this - b); }
		interval_number& operator*=(const double b) { return operator=(*this * b); }
		interval_number& operator/=(const double b) { return operator=(*this / b); }

		interval_number abs() const;

		interval_number sqr() const;

		interval_number pow(unsigned int e) const;

		friend inline interval_number min(const interval_number& a, const interval_number& b);
		friend inline interval_number max(const interval_number& a, const interval_number& b);

		bool operator<(const double b) const { return sup() < b; }
		bool operator<=(const double b) const { return sup() <= b; }
		bool operator>(const double b) const { return inf() > b; }
		bool operator>=(const double b) const { return inf() >= b; }
		bool operator==(const double b) const { return sup()==inf() && sup()==b; }

		void negate() { interval = _mm_shuffle_pd(interval, interval, 1); }

		bool isNegative() const { return _mm_comilt_sd(interval, zero()); }
		bool isPositive() const { return _mm_comilt_sd(_mm_shuffle_pd(interval, interval, 1), zero()); }

#ifdef USE_AVX2_INSTRUCTIONS
		int sign() const {
			__m128d m = _mm_cmplt_sd(interval, zero());
			__m128i r = _mm_castpd_si128(_mm_blendv_pd(_mm_castsi128_pd(_mm_set1_epi32(1)), _mm_castsi128_pd(_mm_set1_epi32(-1)), m));
			return _mm_cvtsi128_si32(r);
		} // Zero is not accounted for

		//int sign() const {
		//	__m128d mp = _mm_cmplt_sd(_mm_shuffle_pd(interval, interval, 1), zero());
		//	__m128d rp = _mm_blendv_pd(zero(), _mm_castsi128_pd(_mm_set1_epi32(1)), mp);
		//	__m128d m = _mm_cmplt_sd(interval, zero());
		//	__m128i r = _mm_castpd_si128(_mm_blendv_pd(rp, _mm_castsi128_pd(_mm_set1_epi32(-1)), m));
		//	return _mm_cvtsi128_si32(r);
		//} // Zero is accounted for
#else
		int sign() const { return (isNegative()) ? (-1) : (1); } // Zero is not accounted for
#endif

#else // USE_SIMD_INSTRUCTIONS
		typedef union error_approx_type_t
		{
			double d;
			uint64_t u;

			inline error_approx_type_t() {}
			inline error_approx_type_t(double a) : d(a) {}
			inline uint64_t is_negative() const { return u >> 63; }
		} casted_double;

	public:
		double min_low, high;

		const double* getInterval() const { return (const double*)&min_low; }

		interval_number() { }
		interval_number(double a) : min_low(-a), high(a) {}
		interval_number(double minf, double sup) : min_low(minf), high(sup) {}
		interval_number(const interval_number& b) : min_low(b.min_low), high(b.high) {}

		double inf() const { return -min_low; }
		double sup() const { return high; }

		bool isNegative() const { return (high < 0); }
		bool isPositive() const { return (min_low < 0); }
		void negate() { std::swap(min_low, high); }

		bool operator<(const double b) const { return (high < b); }

		interval_number& operator=(const interval_number& b) { min_low = b.min_low; high = b.high; return *this; }

		interval_number operator+(const interval_number& b) const { return interval_number(min_low + b.min_low, high + b.high); }

		interval_number operator-(const interval_number& b) const { return interval_number(b.high + min_low, high + b.min_low); }

		interval_number operator*(const interval_number& b) const;

		interval_number operator-() const { return interval_number(high, min_low); }
		interval_number operator+(const double b) const { return interval_number(min_low - b, high + b); }
		interval_number operator-(const double b) const { return interval_number(min_low + b, high - b); }
		interval_number operator*(const double b) const;
		interval_number operator/(const double b) const;

		interval_number& operator+=(const interval_number& b) { min_low += b.min_low; high += b.high; return *this; }
		interval_number& operator-=(const interval_number& b) { return operator=(*this - b); }
		interval_number& operator*=(const interval_number& b) { return operator=(*this * b); }
		interval_number& operator+=(const double b) { min_low -= b; high += b; return *this; }
		interval_number& operator-=(const double b) { min_low += b; high -= b; return *this; }
		interval_number& operator*=(const double b) { return operator=(*this * b); }
		interval_number& operator/=(const double b) { return operator=(*this / b); }

		interval_number abs() const;

		interval_number sqr() const;

		interval_number pow(unsigned int e) const;

		friend inline interval_number min(const interval_number& a, const interval_number& b) {
			return interval_number(std::max(a.min_low, b.min_low), std::min(a.high, b.high));
		}

		friend inline interval_number max(const interval_number& a, const interval_number& b) {
			return interval_number(std::min(a.min_low, b.min_low), std::max(a.high, b.high));
		}

		bool operator>(const double b) const { return (min_low < -b); }
		bool operator==(const double b) const { return (high == b && min_low == -b); }

		int sign() const { return (isNegative()) ? (-1) : (1); } // Zero is not accounted for

#endif // USE_SIMD_INSTRUCTIONS
		double width() const { return sup() - inf(); }

		bool signIsReliable() const { return (isNegative() || isPositive()); } // Zero is not accounted for
		bool containsZero() const { return !signIsReliable(); }

		bool isNAN() const { return sup() != sup(); }

		inline double getMid() const { return (inf() + sup()) / 2; }
		inline bool isExact() const { return inf() == sup(); }

		//inline void operator+=(const interval_number& b) { *this = operator+(b); }

		// Can be TRUE only if the intervals are disjoint
		inline bool operator<(const interval_number& b) const { return (sup() < b.inf()); }
		inline bool operator>(const interval_number& b) const { return (inf() > b.sup()); }

		// Can be TRUE only if the interval interiors are disjoint
		inline bool operator<=(const interval_number& b) const { return (sup() <= b.inf()); }
		inline bool operator>=(const interval_number& b) const { return (inf() >= b.sup()); }

		// TRUE if the intervals are identical single values
		inline bool operator==(const interval_number& b) const { return (sup() == inf() && sup() == b.inf() && sup() == b.sup()); }

		// TRUE if the intervals have no common values
		inline bool operator!=(const interval_number& b) const { return operator<(b) || operator>(b); }

		// The inverse of an interval. Returns NAN if the interval contains zero
		interval_number inverse() const;
	};


	// The square root of an interval
	// Returns NAN if the interval contains a negative value
	inline interval_number sqrt(const interval_number& p)
	{
		const double inf = p.inf();
		const double sup = p.sup();
		if (inf < 0 || sup < 0) return interval_number(NAN);
		const double srinf = sqrt(inf);
		const double srsup = sqrt(sup);
		if (srinf * srinf > inf) return interval_number((-nextafter(srinf, 0)), srsup);
		else return interval_number(-srinf, srsup);
	}

	inline std::ostream& operator<<(std::ostream& os, const interval_number& p)
	{
		os << "[ " << p.inf() << ", " << p.sup() << " ]";
		return os;
	}


	/////////////////////////////////////////////////////////////////////
	// 	   
	// 	   E X P A N S I O N   A R I T H M E T I C
	// 
	/////////////////////////////////////////////////////////////////////

	// Allocate extra-memory
	//#define AllocDoubles(n) ((double *)malloc((n) * sizeof(double)))
	//#define FreeDoubles(p) (free(p))
	#define AllocDoubles(n) ((double *)expansionObject::mempool.alloc((n) * sizeof(double)))
	#define FreeDoubles(p) (expansionObject::mempool.release(p))

	// An instance of the following must be created to access functions for expansion arithmetic
	class expansionObject
	{
	public:
		inline static thread_local MultiPool mempool = MultiPool(2048, 64);

		static void Quick_Two_Sum(const double a, const double b, double& x, double& y) { x = a + b; y = b - (x - a); }

		static void Two_Sum(const double a, const double b, double& x, double& y) {
			double _bv;
			x = a + b; _bv = x - a; y = (a - (x - _bv)) + (b - _bv); 
		}

		static void Two_One_Sum(const double a1, const double a0, const double b, double& x2, double& x1, double& x0) {
			double _i;
			Two_Sum(a0, b, _i, x0); Two_Sum(a1, _i, x2, x1);
		}

		static void two_Sum(const double a, const double b, double* xy) { Two_Sum(a, b, xy[1], xy[0]); }

		static void Two_Diff(const double a, const double b, double& x, double& y) {
			double _bv;
			x = a - b; _bv = a - x; y = (a - (x + _bv)) + (_bv - b); 
		}

		static void Two_One_Diff(const double a1, const double a0, const double b, double& x2, double& x1, double& x0) {
			double _i;
			Two_Diff(a0, b, _i, x0); Two_Sum(a1, _i, x2, x1);
		}

		static void two_Diff(const double a, const double b, double* xy) { Two_Diff(a, b, xy[1], xy[0]); }

		// Products
#ifndef USE_AVX2_INSTRUCTIONS 
		static void Split(double a, double& _ah, double& _al) {
			double _c = 1.3421772800000003e+008 * a;
			_ah = _c - (_c - a); _al = a - _ah;
		}
			
		static void Two_Prod_PreSplit(double a, double b, double _bh, double _bl, double& x, double& y) {
			double _ah, _al;
			x = a * b; 
			Split(a, _ah, _al); 
			y = (_al * _bl) - (((x - (_ah * _bh)) - (_al * _bh)) - (_ah * _bl));
		}
		
		static void Two_Product_2Presplit(double a, double _ah, double _al, double b, double _bh, double _bl, double& x, double& y) {
			x = a * b; y = (_al * _bl) - (((x - _ah * _bh) - (_al * _bh)) - (_ah * _bl));
		}
#endif

		// [x,y] = [a]*[b]		 Multiplies two expansions [a] and [b] of length one
		static void Two_Prod(const double a, const double b, double& x, double& y);
		static void Two_Prod(const double a, const double b, double* xy) { Two_Prod(a, b, xy[1], xy[0]); }


		// [x,y] = [a]^2		Squares an expansion of length one
		static void Square(const double a, double& x, double& y);
		static void Square(const double a, double* xy) { Square(a, xy[1], xy[0]); }

		// [x2,x1,x0] = [a1,a0]-[b]		Subtracts an expansion [b] of length one from an expansion [a1,a0] of length two
		static void two_One_Diff(const double a1, const double a0, const double b, double& x2, double& x1, double& x0)
		 { Two_One_Diff(a1, a0, b, x2, x1, x0); }

		static void two_One_Diff(const double* a, const double b, double* x) { two_One_Diff(a[1], a[0], b, x[2], x[1], x[0]); }

		// [x3,x2,x1,x0] = [a1,a0]*[b]		Multiplies an expansion [a1,a0] of length two by an expansion [b] of length one
		static void Two_One_Prod(const double a1, const double a0, const double b, double& x3, double& x2, double& x1, double& x0);
		static void Two_One_Prod(const double* a, const double b, double* x) { Two_One_Prod(a[1], a[0], b, x[3], x[2], x[1], x[0]); }

		// [x3,x2,x1,x0] = [a1,a0]+[b1,b0]		Calculates the sum of two expansions of length two
		static void Two_Two_Sum(const double a1, const double a0, const double b1, const double b0, double& x3, double& x2, double& x1, double& x0) {
			double _j, _0;
			Two_One_Sum(a1, a0, b0, _j, _0, x0); Two_One_Sum(_j, _0, b1, x3, x2, x1);	
		}
		
		static void Two_Two_Sum(const double* a, const double* b, double* xy) { Two_Two_Sum(a[1], a[0], b[1], b[0], xy[3], xy[2], xy[1], xy[0]); }

		// [x3,x2,x1,x0] = [a1,a0]-[b1,b0]		Calculates the difference between two expansions of length two
		static void Two_Two_Diff(const double a1, const double a0, const double b1, const double b0, double& x3, double& x2, double& x1, double& x0) {
			double _j, _0, _u3;
			Two_One_Diff(a1, a0, b0, _j, _0, x0); Two_One_Diff(_j, _0, b1, _u3, x2, x1); x3 = _u3; 
		}

		static void Two_Two_Diff(const double* a, const double* b, double* x) { Two_Two_Diff(a[1], a[0], b[1], b[0], x[3], x[2], x[1], x[0]); }

		// Calculates the second component 'y' of the expansion [x,y] = [a]-[b] when 'x' is known
		static void Two_Diff_Back(const double a, const double b, double& x, double& y) { 
			double _bv;
			_bv = a - x; y = (a - (x + _bv)) + (_bv - b); 
		}

		static void Two_Diff_Back(const double a, const double b, double* xy) { Two_Diff_Back(a, b, xy[1], xy[0]); }

		// [h] = [a1,a0]^2		Squares an expansion of length 2
		// 'h' must be allocated by the caller with 6 components.
		static void Two_Square(const double& a1, const double& a0, double* x);

		// [h7,h6,...,h0] = [a1,a0]*[b1,b0]		Calculates the product of two expansions of length two.
		// 'h' must be allocated by the caller with eight components.
		static void Two_Two_Prod(const double a1, const double a0, const double b1, const double b0, double* h);
		static void Two_Two_Prod(const double* a, const double* b, double* xy) { Two_Two_Prod(a[1], a[0], b[1], b[0], xy); }

		// [e] = -[e]		Inplace inversion
		static void Gen_Invert(const int elen, double* e) { for (int i = 0; i < elen; i++) e[i] = -e[i]; }

		// [h] = [e] + [f]		Sums two expansions and returns number of components of result
		// 'h' must be allocated by the caller with at least elen+flen components.
		static int Gen_Sum(const int elen, const double* e, const int flen, const double* f, double* h);

		// Same as above, but 'h' is allocated internally. The caller must still call 'free' to release the memory.
		static int Gen_Sum_With_Alloc(const int elen, const double* e, const int flen, const double* f, double** h)
		{
			*h = AllocDoubles(elen + flen);
			return Gen_Sum(elen, e, flen, f, *h);
		}

		// [h] = [e] + [f]		Subtracts two expansions and returns number of components of result
		// 'h' must be allocated by the caller with at least elen+flen components.
		static int Gen_Diff(const int elen, const double* e, const int flen, const double* f, double* h);

		// Same as above, but 'h' is allocated internally. The caller must still call 'free' to release the memory.
		static int Gen_Diff_With_Alloc(const int elen, const double* e, const int flen, const double* f, double** h)
		{
			*h = AllocDoubles(elen + flen);
			return Gen_Diff(elen, e, flen, f, *h);
		}

		// [h] = [e] * b		Multiplies an expansion by a scalar
		// 'h' must be allocated by the caller with at least elen*2 components.
		static int Gen_Scale(const int elen, const double* e, const double b, double* h);

		// [h] = [e] * 2		Multiplies an expansion by 2
		// 'h' must be allocated by the caller with at least elen components. This is exact up to overflows.
		static void Double(const int elen, const double* e, double* h) { for (int i = 0; i < elen; i++) h[i] = 2 * e[i]; }
		
		// [h] = [e] * n		Multiplies an expansion by n
		// If 'n' is a power of two, the multiplication is exact
		static void ExactScale(const int elen, double* e, const double n) { for (int i = 0; i < elen; i++) e[i] *= n; }

		// [h] = [a] * [b]
		// 'h' must be allocated by the caller with at least 2*alen*blen components.
		static int Sub_product(const int alen, const double* a, const int blen, const double* b, double* h);

		// [h] = [a] * [b]
		// 'h' must be allocated by the caller with at least MAX(2*alen*blen, 8) components.
		static int Gen_Product(const int alen, const double* a, const int blen, const double* b, double* h);

		// Same as above, but 'h' is allocated internally. The caller must still call 'free' to release the memory.
		static int Gen_Product_With_Alloc(const int alen, const double* a, const int blen, const double* b, double** h);


		// Assume that *h is pre-allocated with hlen doubles.
		// If more elements are required, *h is re-allocated internally.
		// In any case, the function returns the size of the resulting expansion.
		// The caller must verify whether reallocation took place, and possibly call 'free' to release the memory.
		// When reallocation takes place, *h becomes different from its original value.

		static int Double_With_PreAlloc(const int elen, const double* e, double** h, const int hlen);

		static int Gen_Scale_With_PreAlloc(const int elen, const double* e, const double& b, double** h, const int hlen);

		static int Gen_Sum_With_PreAlloc(const int elen, const double* e, const int flen, const double* f, double** h, const int hlen);

		static int Gen_Diff_With_PreAlloc(const int elen, const double* e, const int flen, const double* f, double** h, const int hlen);

		static int Gen_Product_With_PreAlloc(const int alen, const double* a, const int blen, const double* b, double** h, const int hlen);

		// Approximates the expansion to a double
		static double To_Double(const int elen, const double* e);

		static void print(const int elen, const double* e) { for (int i = 0; i < elen; i++) printf("%e ", e[i]); printf("\n");}
	};

#ifdef USE_GNU_GMP_CLASSES
	typedef mpz_class bignatural;
	typedef mpq_class bigfloat;
	typedef mpq_class bigrational;
#else
	/////////////////////////////////////////////////////////////////////
	// 	   
	// 	   B I G   N A T U R A L
	// 
	/////////////////////////////////////////////////////////////////////

	// Preallocates memory for bignaturals having at most 32 limbs.
	// Larger numbers will use the standard heap.
	inline static thread_local MultiPool nfgMemoryPool;

	// A bignatural is an arbitrarily large non-negative integer.
	// It is made of a sequence of digits in base 2^32.
	// Leading zero-digits are not allowed.
	// The value 'zero' is represented by an empty digit sequence.

	class bignatural {
		uint32_t* digits;	// Ptr to the digits
		uint32_t m_size;		// Actual number of digits
		uint32_t m_capacity;	// Current vector capacity

		inline static uint32_t* BN_ALLOC(uint32_t num_bytes) { return (uint32_t*)nfgMemoryPool.alloc(num_bytes); }
		inline static void BN_FREE(uint32_t* ptr) { nfgMemoryPool.release(ptr); }

		void init(const bignatural& m);
		void init(const uint64_t m); 

	public:
		// Creates a 'zero'
		bignatural() : digits(NULL), m_size(0), m_capacity(0) { }

		~bignatural() { BN_FREE(digits); }

		bignatural(const bignatural& m) { init(m); }

		bignatural(bignatural&& m) noexcept : digits(m.digits), m_size(m.m_size), m_capacity(m.m_capacity) { 
			m.digits = nullptr;
		}

		// Creates from uint64_t
		bignatural(uint64_t m) { init(m); }

		// If the number fits a uint64_t convert and return true
		bool toUint64(uint64_t& n) const;

		// If the number fits a uint32_t convert and return true
		bool toUint32(uint32_t& n) const;

		bignatural& operator=(const bignatural& m);

		bignatural& operator=(const uint64_t m);

		inline const uint32_t& back() const { return digits[m_size - 1]; }

		inline const uint32_t& operator[](int i) const { return digits[i]; }

		inline uint32_t size() const { return m_size; }

		inline bool empty() const { return m_size == 0; }

		// Left-shift by n bits and possibly add limbs as necessary
		void operator<<=(uint32_t n);

		// Right-shift by n bits
		void operator>>=(uint32_t n);

		bool operator==(const bignatural& b) const;

		bool operator!=(const bignatural& b) const;

		bool operator>=(const bignatural& b) const;

		bool operator>(const bignatural& b) const;

		bool operator<=(const bignatural& b) const { return b >= *this; }

		bool operator<(const bignatural& b) const { return b > *this; }

		bignatural operator+(const bignatural& b) const;

		// Assume that b is smaller than or equal to this number!
		bignatural operator-(const bignatural& b) const;

		bignatural operator*(const bignatural& b) const;

		// Short division algorithm
		bignatural divide_by(const uint32_t D, uint32_t& remainder) const;

		uint32_t getNumSignificantBits() const;

		bool getBit(uint32_t b) const;

		// Long division
		bignatural divide_by(const bignatural& divisor, bignatural& remainder) const;

		// Greatest common divisor (Euclidean algorithm)
		bignatural GCD(const bignatural& D) const;

		// String representation in decimal form
		std::string get_dec_str() const;

		// String representation in binary form
		std::string get_str() const;

		// Count number of zeroes on the right (least significant binary digits)
		uint32_t countEndingZeroes() const;

	protected:
		inline uint32_t& back() { return digits[m_size - 1]; }

		inline void pop_back() { m_size--; }

		inline uint32_t& operator[](int i) { return digits[i]; }

		inline void push_back(uint32_t b) {
			if (m_size == m_capacity) increaseCapacity((m_capacity|1)<<2);
			digits[m_size++] = b;
		}

		inline void push_bit_back(uint32_t b) {
			if (m_size) {
				operator<<=(1);
				back() |= b;
			}
			else if (b) push_back(1);
		}

		inline void reserve(uint32_t n) {
			if (n > m_capacity) increaseCapacity(n);
		}

		inline void resize(uint32_t n) {
			reserve(n);
			m_size = n;
		}

		inline void fill(uint32_t v) {
			for (uint32_t i = 0; i < m_size; i++) digits[(int)i] = v;
		}

		inline void pop_front() {
			for (uint32_t i = 1; i < m_size; i++) digits[i - 1] = digits[i];
			pop_back();
		}

		// Count number of zeroes on the left (most significant digits)
		uint32_t countLeadingZeroes() const;

		// Count number of zeroes on the right of the last 1 in the least significant limb
		// Assumes that number is not zero and last limb is not zero!
		uint32_t countEndingZeroesLSL() const {
			const uint32_t d = back();
			uint32_t i = 31;
			while (!(d << i)) i--;
			return 31 - i;
		}

		inline void pack() {
			uint32_t i = 0;
			while (i < m_size && digits[i] == 0) i++;
			if (i) {
				uint32_t* dold = digits + i;
				uint32_t* dnew = digits;
				uint32_t* dend = digits + m_size;
				while (dold < dend) *dnew++ = *dold++;
				m_size -= i;
			}
		}

		// a and b must NOT be this number!
		void toSum(const bignatural& a, const bignatural& b);

		// a and b must NOT be this number!
		// Assume that b is smaller or equal than a!
		void toDiff(const bignatural& a, const bignatural& b);

		// a and b must NOT be this number!
		void toProd(const bignatural& a, const bignatural& b);

	private:

		// Multiplies by a single limb, left shift, and add to accumulator. Does not pack!
		void addmul(uint32_t b, uint32_t left_shifts, bignatural& result) const;

		void increaseCapacity(uint32_t new_capacity);

		friend class bigfloat;
	};

	inline std::ostream& operator<<(std::ostream& os, const bignatural& p)
	{
		os << p.get_dec_str();
		return os;
	}

	/////////////////////////////////////////////////////////////////////
	// 	   
	// 	   B I G   F L O A T
	// 
	/////////////////////////////////////////////////////////////////////

	// A bigfloat is a floting point number with arbitrarily large mantissa.
	// In principle, we could have made the exponent arbitrarily large too,
	// but in practice this appears to be useless.
	// Exponents are in the range [-INT32_MAX, INT32_MAX]
	//
	// A bigfloat f evaluates to f = sign * mantissa * 2^exponent
	//
	// mantissa is a bignatural whose least significant bit is 1.
	// Number is zero if mantissa is empty.

	class bigfloat {
		bignatural mantissa; // .back() is less significant. Use 32-bit limbs to avoid overflows using 64-bits
		int32_t exponent; // In principle we might still have under/overflows, but not in practice
		int32_t sign;	// Redundant but keeps alignment

	public:
		// Default constructor creates a zero-valued bigfloat
		bigfloat() : exponent(0), sign(0) {}

		// Lossless conversion from double
		bigfloat(const double d);

		// Truncated approximation
		double get_d() const;
		
		bigfloat operator+(const bigfloat& b) const;

		bigfloat operator-(const bigfloat& b) const;

		bigfloat operator*(const bigfloat& b) const;

		void invert() { sign = -sign; }

		bigfloat inverse() const {
			bigfloat r = *this;
			r.invert();
			return r;
		}

		inline int sgn() const { return sign; }

		std::string get_str() const;

		const bignatural& getMantissa() const { return mantissa; }
		int32_t getExponent() const { return exponent; }

	private:

		// Right-shift as long as the least significant bit is zero
		void pack();

		// Left-shift the mantissa by n bits and reduce the exponent accordingly
		void leftShift(uint32_t n) {
			mantissa <<= n;
			exponent -= n;
		}
	};

	inline int sgn(const bigfloat& f) {
		return f.sgn();
	}

	inline bigfloat operator-(const bigfloat& f) {
		return f.inverse();
	}

	inline bigfloat operator*(double d, const bigfloat& f) {
		return f * bigfloat(d);
	}

	inline std::ostream& operator<<(std::ostream& os, const bigfloat& p)
	{
		if (p.sgn() < 0) os << "-";
		os << p.getMantissa().get_dec_str() << " * 2^" << p.getExponent();
		return os;
	}

/////////////////////////////////////////////////////////////////////
// 	   
// 	   B I G   R A T I O N A L
// 
/////////////////////////////////////////////////////////////////////

// A bigrational is a fraction of two coprime bignaturals with a sign.
// Number is zero if sign is zero

	class bigrational {
		bignatural numerator, denominator;
		int32_t sign;	// Redundant but keeps alignment

		// Iteratively divide both num and den by two as long as they are both even
		void compress();

		// Make numerator and denominator coprime (divide both by GCD)
		void canonicalize();

	public:
		// Create a zero
		bigrational() : sign(0) {}

		// Create from a bigfloat (lossless)
		bigrational(const bigfloat& f);

		// Create from explicit numerator, denominator and sign.
		bigrational(const bignatural& num, const bignatural& den, int32_t s) :
			numerator(num), denominator(den), sign(s) {
			canonicalize();
		}

		// Convert to multiplicative inverse
		void invert() {
			if (!sign) ip_error("bigrational::invert() : inverse zero!\n");
			std::swap(numerator, denominator);
		}

		// Return multiplicative inverse
		bigrational inverse() const { bigrational r = *this; r.invert(); return r; }

		// Invert sign
		void negate() { sign = -sign; }

		// Return additive inverse
		bigrational negation() const { bigrational r = *this; r.negate(); return r; }

		// Standard arithmetic ops
		bigrational operator+(const bigrational& r) const;

		bigrational operator-(const bigrational& r) const { return operator+(r.negation()); }

		bigrational operator*(const bigrational& r) const {
			if (sign == 0 || r.sign == 0) return bigrational();
			else return bigrational(numerator * r.numerator, denominator * r.denominator, sign * r.sign);
		}

		bigrational operator/(const bigrational& r) const {
			if (!r.sign) ip_error("bigrational::operator/ : division by zero!\n");
			return operator*(r.inverse()); 
		}

		// Comparison operators
		bool operator==(const bigrational& r) const {
			return (sign == r.sign && numerator == r.numerator && denominator == r.denominator);
		}

		bool operator!=(const bigrational& r) const {
			return (sign != r.sign || numerator != r.numerator || denominator != r.denominator);
		}

		bool hasGreaterModule(const bigrational& r) const {
			return numerator * r.denominator > r.numerator * denominator;
		}

		bool hasGrtrOrEqModule(const bigrational& r) const {
			return numerator * r.denominator >= r.numerator * denominator;
		}

		bool operator>(const bigrational& r) const {
			return (sign > r.sign || (sign > 0 && r.sign > 0 && hasGreaterModule(r)) || (sign < 0 && r.sign < 0 && r.hasGreaterModule(*this)));
		}

		bool operator>=(const bigrational& r) const {
			return (sign > r.sign || (sign > 0 && r.sign > 0 && hasGrtrOrEqModule(r)) || (sign < 0 && r.sign < 0 && r.hasGrtrOrEqModule(*this)));
		}

		bool operator<(const bigrational& r) const {
			return (sign < r.sign || (sign < 0 && r.sign < 0 && hasGreaterModule(r)) || (sign > 0 && r.sign > 0 && r.hasGreaterModule(*this)));
		}

		bool operator<=(const bigrational& r) const {
			return (sign < r.sign || (sign < 0 && r.sign < 0 && hasGrtrOrEqModule(r)) || (sign > 0 && r.sign > 0 && r.hasGrtrOrEqModule(*this)));
		}

		// Conversion to double (truncated)
		double get_d() const;

		const bignatural& get_num() const { return numerator; }
		const bignatural& get_den() const { return denominator; }
		int32_t sgn() const { return sign; }

		// Return decimal representation
		std::string get_dec_str() const;

		// Return binary representation
		std::string get_str() const;
	};

	inline bigrational operator-(const bigrational& p) { return p.negation(); }

	inline int32_t sgn(const bigrational& p) { return p.sgn(); }

	inline std::ostream& operator<<(std::ostream& os, const bigrational& p)
	{
		os << p.get_dec_str();
		return os;
	}
#endif // USE_GNU_GMP_CLASSES

#include "numerics.hpp"

#endif //NUMERICS_H