// Generic Geometry Library - Numeric Adaptor
//
// Copyright Barend Gehrels 1995-2009, Geodan Holding B.V. Amsterdam, the Netherlands.
// Copyright Bruno Lalande 2008, 2009
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#include "gmp.h"

#include <boost/numeric/conversion/cast.hpp>

#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>



//-------------------------------------------------------------------------------------------------
// Numeric adaptor, adapts for numbers (adaption implementation in POLICY)
//-------------------------------------------------------------------------------------------------
template <typename POLICY>
struct numeric_adaptor
{
	inline numeric_adaptor()
	{
		POLICY::init(value);
	}

	// Copy constructor
	inline numeric_adaptor(const numeric_adaptor<POLICY>& v)
	{
		POLICY::init(value);
		POLICY::copy(v.value, value);
	}

	// Constructor with a normal IEEE type 
	template <typename T>
	inline numeric_adaptor(const T& v)
	{
		POLICY::init(value);
		POLICY::template set<T>(value, v);
	}

	virtual ~numeric_adaptor()
	{
		POLICY::destruct(value);
	}

	// Assignment from other value
	inline numeric_adaptor<POLICY> operator=(const numeric_adaptor<POLICY>& v)
	{
		POLICY::copy(v.value, this->value);
		return *this;
	}

	// Assignment from normal IEEE type
	template <typename T>
	inline numeric_adaptor<POLICY> operator=(const T& v)
	{
		POLICY::template set<T>(this->value, v);
		return *this;
	}

	// Cast to normal IEEE type
	template <typename T>
	inline operator T() const
	{
		return POLICY::template big_numeric_cast<T>(value);
	}


	// Comparisons
	inline bool operator<(const numeric_adaptor<POLICY>& other) const
	{
		return POLICY::compare(value, other.value) < 0;
	}

	inline bool operator>(const numeric_adaptor<POLICY>& other) const
	{
		return POLICY::compare(value, other.value) > 0;
	}

	inline bool operator==(const numeric_adaptor<POLICY>& other) const
	{
		return POLICY::compare(value, other.value) == 0;
	}

	// Operators
	friend inline numeric_adaptor<POLICY> operator+(
			const numeric_adaptor<POLICY>& a, 
			const numeric_adaptor<POLICY>& b)
	{
		typename POLICY::type r;
		POLICY::init(r);
		POLICY::add(r, a.value, b.value);
		return numeric_adaptor<POLICY>(r, true);
	}

	friend inline numeric_adaptor<POLICY> operator*(
			const numeric_adaptor<POLICY>& a, 
			const numeric_adaptor<POLICY>& b)
	{
		typename POLICY::type r;
		POLICY::init(r);
		POLICY::multiply(r, a.value, b.value);
		return numeric_adaptor<POLICY>(r, true);
	}

	friend inline numeric_adaptor<POLICY> operator-(
			const numeric_adaptor<POLICY>& a, 
			const numeric_adaptor<POLICY>& b)
	{
		typename POLICY::type r;
		POLICY::init(r);
		POLICY::subtract(r, a.value, b.value);
		return numeric_adaptor<POLICY>(r, true);
	}

	friend inline numeric_adaptor<POLICY> operator/(
			const numeric_adaptor<POLICY>& a, 
			const numeric_adaptor<POLICY>& b)
	{
		typename POLICY::type r;
		POLICY::init(r);
		POLICY::divide(r, a.value, b.value);
		return numeric_adaptor<POLICY>(r, true);
	}

	// Functions
	static inline numeric_adaptor<POLICY> sqrt(const numeric_adaptor<POLICY>& v)
	{
		typename POLICY::type r;
		POLICY::init(r);
		POLICY::sqrt(r, v.value);
		return numeric_adaptor<POLICY>(r, true);
	}

	private :
		typename POLICY::type value;

		// Constructor with a type. Bool (or any other signature changing parameter)
		// is necessary for cases where type == CT
		inline numeric_adaptor<POLICY>(const typename POLICY::type& v, bool)
		{
			POLICY::init(value);
			POLICY::copy(v, value);
		}

};


//-------------------------------------------------------------------------------------------------
// Policy for GMP big number library
//-------------------------------------------------------------------------------------------------
struct gmp_policy
{
	typedef mpf_t type;

	static inline void init(type& value)
	{
		mpf_init(value);
	}

	static inline void destruct(type& value)
	{
		mpf_clear(value);
	}

	template <typename CT>
	static inline void set(type& value, const CT& v)
	{
		mpf_set_d(value, v);
	}

	static inline void copy(const type& source, type& dest)
	{
		mpf_set(dest, source);
	}

	static inline void add(type& r, const type& a, const type& b)
	{
		mpf_add(r, a, b);
	}

	static inline void subtract(type& r, const type& a, const type& b)
	{
		mpf_sub(r, a, b);
	}

	static inline void multiply(type& r, const type& a, const type& b)
	{
		mpf_mul(r, a, b);
	}

	static inline void divide(type& r, const type& a, const type& b)
	{
		mpf_div(r, a, b);
	}

	static inline void sqrt(type& r, const type& a)
	{
		mpf_sqrt(r, a);
	}

	template <typename CT>
	static inline CT big_numeric_cast(const type& b)
	{
		return mpf_get_d(b);
	}

	static inline int compare(const type& a, const type& b)
	{
		return mpf_cmp(a, b);
	}
};


//-------------------------------------------------------------------------------------------------
// Policy for normal IEEE floating point arithmetic
//-------------------------------------------------------------------------------------------------
template <typename T>
struct ieee_policy
{
	typedef T type;

	static inline void init(type& value) {}
	static inline void destruct(type& value) {}

	template <typename CT>
	static inline void set (type& value, const CT& v)	{ value = boost::numeric_cast<T>(v); }

	static inline void copy(const type& source, type& dest) { dest = source; }

	static inline void add(type& r, const type& a, const type& b) { r = a + b; }
	static inline void subtract(type& r, const type& a, const type& b) { r = a - b; }
	static inline void multiply(type& r, const type& a, const type& b) { r = a * b; }
	static inline void divide(type& r, const type& a, const type& b) { r = a / b; }

	static inline void sqrt(type& r, const type& a) { r = ::sqrt(a); }

	template <typename CT>
	static inline CT big_numeric_cast(const type& v)
	{
		return boost::numeric_cast<CT>(v);
	}

	static inline int compare(const type& a, const type& b)
	{
		return a < b ? -1 : a > b ? 1 : 0;
	}
};

//-------------------------------------------------------------------------------------------------
// Policy for CLN big number library - not worked out
//-------------------------------------------------------------------------------------------------
struct cln_policy
{
};

//-------------------------------------------------------------------------------------------------
// Policy for mpfr big number library - not worked out
//-------------------------------------------------------------------------------------------------
struct mpfr_policy
{
};


//-------------------------------------------------------------------------------------------------
// sample algorithm, using NUM as numeric adaptor, T as normal IEEE type
//-------------------------------------------------------------------------------------------------
template <typename T, typename NUM>
void sample(const std::string& header, T c1, T c2, T c4, T ta, T tb, T tc)
{
	std::cout << std::endl << "---" << header << std::endl;
	NUM v1 = c1;
	NUM v2 = c2;
	NUM v3 = v1 + v2;
	std::cout << T(v3) << std::endl;

	NUM v4 = c4;
	NUM v5 = (v1 + v2) * v4;
	std::cout << T(v5) << std::endl;
	v5 = v1 + v2 * v4;
	std::cout << T(v5) << std::endl;

	NUM v6 = NUM::sqrt(v3);
	std::cout << T(v6) << std::endl;

	v6 = 4;
	std::cout << T(v6) << std::endl;

	v6 = v2;
	std::cout << T(v6) << std::endl;
	std::cout << std::endl;

	if (v1 > v2)
	{
		std::cout << "v1 > v2" << std::endl;
	}
	if (v1 < v2)
	{
		std::cout << "v1 < v2" << std::endl;
	}
	if (v1 == v2)
	{
		std::cout << "v1 == v2" << std::endl;
	}

	// Test Heron formule
	{
		NUM a = ta;
		NUM b = tb;
		NUM c = tc;
		NUM s((a + b + c) / NUM(2.0));
		NUM area = NUM::sqrt(s * (s - a) * (s - b) * (s - c));
		std::cout << "area: " << T(area) << std::endl;
	}

}



int main()
{
	std::cout << std::setprecision(12);
	double a = 31622.77662;
	double b = 0.000023;
	double c = 31622.77661;

	sample<float, numeric_adaptor<ieee_policy<float> > >("use float, calculate with float", 2.0, 3.0, 5.0, a, b, c);
	sample<float, numeric_adaptor<ieee_policy<double> > >("use float, calculate with double", 2.0, 3.0, 5.0, a, b, c);
	sample<double, numeric_adaptor<ieee_policy<long double> > >("use double, calculate with long double", 2.0, 3.0, 5.0, a, b, c);
	sample<double, numeric_adaptor<gmp_policy> >("use double, calculate with gmp", 2.0, 3.0, 5.0, a, b, c);
	return 0;
}