#ifndef _BIGFLOAT_H
#define _BIGFLOAT_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef 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 General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//
// bigfloat<N> class with N digits precision
// as a remplacement for the double
//
//    class bigfloat<N>
//    class numeric_limits <bigfloat<N> >
//    class numeric_flags  <bigfloat<N> >
//
// author: Pierre.Saramito@imag.fr
//
// date: 25 january 2000
//
// IMPLEMENTATION NOTE:
// this is a wrapper class for cln_F:
// provides an interface comparable to double.
//
// TODO: ostream
//   knows setprecision(n) but 
// does not respond to 
//   scientific & fixed & uppercase
//
// TODO: complete the math lib, as IEEE
// -> see the blitz test suite.
//

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <cmath>
#include <cln/cln.h>
#ifdef _RHEOLEF_HAVE_GINAC
#include <ginac/ginac.h>
#endif // _RHEOLEF_HAVE_GINAC
#include "rheolef/numeric_limits.h"
#include "rheolef/numeric_flags.h"
#include "all_float_io.h"
#include <sstream>

template <int N>
class bigfloat {
public:

// allocators/deallocators

  bigfloat(double x = 0);
  bigfloat(const char* s);
  bigfloat(int x);
  bigfloat(long x);
  bigfloat(unsigned int x);
  bigfloat(unsigned long x);

// type conversions

  operator double() const; // compiler do not send any warning, such as
  operator int() const;    // warning: initialization to `int' from `bigfloat' !!
#ifdef _RHEOLEF_HAVE_GINAC
  operator GiNaC::ex() const;
#endif // _RHEOLEF_HAVE_GINAC

// assignments

  bigfloat<N>& operator = (const bigfloat<N>&);
  bigfloat<N>& operator = (const char*&);
  bigfloat<N>& operator = (const double&);
  bigfloat<N>& operator = (const int&);
  bigfloat<N>& operator = (const long&);
  bigfloat<N>& operator = (const unsigned int&);
  bigfloat<N>& operator = (const unsigned long&);

// operators
  template <int M>
  friend bigfloat<M> operator - (const bigfloat<M>&);

# define bigfloat_binary_operator_basic(op,T1,T2) 		\
  template <int M>						\
  friend bigfloat<M> operator op (const T1&, const T2&);

# define bigfloat_binary_operator(op) 				\
bigfloat_binary_operator_basic(op,bigfloat<M>,bigfloat<M>)	\
bigfloat_binary_operator_basic(op,bigfloat<M>,double)		\
bigfloat_binary_operator_basic(op,bigfloat<M>,int)		\
bigfloat_binary_operator_basic(op,bigfloat<M>,long)		\
bigfloat_binary_operator_basic(op,bigfloat<M>,unsigned int)	\
bigfloat_binary_operator_basic(op,bigfloat<M>,unsigned long)	\
bigfloat_binary_operator_basic(op,double,bigfloat<M>)		\
bigfloat_binary_operator_basic(op,int,bigfloat<M>)		\
bigfloat_binary_operator_basic(op,long,bigfloat<M>)		\
bigfloat_binary_operator_basic(op,unsigned int,bigfloat<M>)	\
bigfloat_binary_operator_basic(op,unsigned long,bigfloat<M>)

bigfloat_binary_operator(+)
bigfloat_binary_operator(-)
bigfloat_binary_operator(*)
bigfloat_binary_operator(/)

# undef bigfloat_binary_operator_basic
# undef bigfloat_binary_operator

// computed assignments

# define bigfloat_computed_assignment(op) 			\
  bigfloat<N>& operator op(const bigfloat<N>&);			\
  bigfloat<N>& operator op(const double&);			\
  bigfloat<N>& operator op(const int&);				\
  bigfloat<N>& operator op(const long&);			\
  bigfloat<N>& operator op(const unsigned int&);		\
  bigfloat<N>& operator op(const unsigned long&);

bigfloat_computed_assignment(+=)
bigfloat_computed_assignment(-=)
bigfloat_computed_assignment(*=)
bigfloat_computed_assignment(/=)

#undef bigfloat_computed_assignment

// comparators

#define bigfloat_comparator_basic(op)					\
  template <int M>							\
  friend bool operator op (const bigfloat<M>&, const bigfloat<M>&);
#define bigfloat_comparator_both(op,T)					\
  template <int M>							\
  friend bool operator op (const bigfloat<M>&, const T&);		\
  template <int M>							\
  friend bool operator op (const T&, const bigfloat<M>&);
#define bigfloat_comparator(op)						\
  bigfloat_comparator_basic(op)						\
  bigfloat_comparator_both(op,double)					\
  bigfloat_comparator_both(op,int)					\
  bigfloat_comparator_both(op,long)					\
  bigfloat_comparator_both(op,unsigned int)				\
  bigfloat_comparator_both(op,unsigned long)

bigfloat_comparator(==)
bigfloat_comparator(!=)
bigfloat_comparator(<)
bigfloat_comparator(>)
bigfloat_comparator(<=)
bigfloat_comparator(>=)

#undef bigfloat_comparator
#undef bigfloat_comparator_basic
#undef bigfloat_comparator_both

// math functions

  template <int M> friend bigfloat<M> fabs(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> sqrt(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> sin(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> cos(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> tan(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> log(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> exp(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> sinh(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> cosh(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> tanh(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> atan(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> acos(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> sqr(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> floor(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> ceil(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> trunc(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> round(const bigfloat<M>& x);

  template <int M> friend bigfloat<M> log10(const bigfloat<M>& x);
  template <int M> friend bigfloat<M> pow(const bigfloat<M>& x, const bigfloat<M>& y);
  template <int M> friend bigfloat<M> pow(const bigfloat<M>& x, int n);

// inputs/outputs

  template <int M> friend std::ostream& operator <<  (std::ostream& s, const bigfloat<M>& x);
  template <int M> friend std::istream& operator >>  (std::istream& s, bigfloat<M>& x);

// constants (may be const! -> core dump..)

#ifdef TODO
  static bigfloat<N> Pi;
  static bigfloat<N> one_on_log10;
#endif // TODO

// numeric limits

// Returns the smallest float format which guarantees at least
// n decimal digits in the mantissa (after the decimal point).
// Methode from: cln-1.0.2/src/float/misc/cl_float_format.cc
// Mindestens 1+n Dezimalstellen (inklusive Vorkommastelle)
// bedeutet mindestens ceiling((1+n)*ln(10)/ln(2)) Binarstellen.
// ln(10)/ln(2) = 3.321928095 = (binar) 11.01010010011010011110000100101111...
//                       = (binar) 100 - 0.10101101100101100001111011010001
// Durch diese Berechnungsmethode wird das Ergebnis sicher >= (1+n)*ln(10)/ln(2)
// sein, evtl. um ein paar Bit zu gross, aber nicht zu klein.

  static const int digits
        = ((1+N) << 2) 
	- ((1+N) >> 1)  - ((1+N) >> 3)  - ((1+N) >> 5)
	- ((1+N) >> 6)  - ((1+N) >> 8)  - ((1+N) >> 9)
	- ((1+N) >> 12) - ((1+N) >> 14) - ((1+N) >> 15);

  static const int digits10 		= N;

  static bigfloat<N> epsilon();
  static bigfloat<N> min();
  static bigfloat<N> max();

// computation of some constants

  static bigfloat<N> compute_pi();

// data

protected:

  cln::cl_F _f;

  // internal constructor
  bigfloat(const cln::cl_F& y) : _f(y) {}

public:
  // for debug
  std::ostream& dump(std::ostream&);
};
//
// TODO: complete the fields...
//
namespace std {
template<int N>
class numeric_limits <bigfloat<N> > {
public:
    static const int  digits 	     = bigfloat<N>::digits;
    static const int  digits10 	     = bigfloat<N>::digits10;

    // have no representation...
    static bigfloat<N> min ()          { return bigfloat<N>::min(); }
    static bigfloat<N> max ()          { return bigfloat<N>::max(); }

    static bigfloat<N> epsilon ()      { return bigfloat<N>::epsilon(); }
    static bigfloat<N> denorm_min ()   { return bigfloat<N>::min(); }

    // CLN does not implement features like NaNs, denormalized numbers
    // and gradual underflow.
    // If the exponent range of some floating-point type is too limited 
    // for your application, choose another floating-point type with
    // larger exponent range. 

#ifdef TODO
    static const int  min_exponent   = INT_MIN;
    static const int  max_exponent   = INT_MAX;
    static const int  min_exponent10 = INT_MIN;
    static const int  max_exponent10 = INT_MAX;

    __NUMERIC_LIMITS_FLOAT(<bigfloat<N>)

#endif // TODO
};
} // namespace std
template <int N>
class numeric_flags<bigfloat<N> > {
public:
    static bool output_as_double()       { return _output_as_double; }
    static void output_as_double(bool f) { _output_as_double = f; }
protected:
    static bool _output_as_double;
};

// flag initialisation

template <int N>
bool numeric_flags<bigfloat<N> >::_output_as_double = false;


// numeric limits

template <int N>
inline
bigfloat<N>
bigfloat<N>::min() 
{
    return bigfloat<N>(most_negative_float(cln::float_format_t(digits)));
}
template <int N>
inline
bigfloat<N>
bigfloat<N>::max() 
{
    return bigfloat<N>(most_positive_float(cln::float_format_t(digits)));
}
template <int N>
inline
bigfloat<N>
bigfloat<N>::epsilon() 
{
    return bigfloat<N>(float_epsilon(cln::float_format_t(digits)));
}

// constants

template <int N>
inline
bigfloat<N> 
bigfloat<N>::compute_pi()
{
	return bigfloat<N>(cln::pi(cln::float_format_t(digits)));
}

// allocators/deallocators

#define bigfloat_cstor(T)				\
template <int N>					\
inline							\
bigfloat<N>::bigfloat(T x)				\
	: _f(cln::cl_float(x, cln::float_format_t(digits)))	\
{							\
}
bigfloat_cstor(double)
bigfloat_cstor(int)
bigfloat_cstor(unsigned int)
#undef  bigfloat_cstor

template <int N>
inline
bigfloat<N>::bigfloat(const char* s)
	: _f(cln::cl_float(cln::cl_F(s), cln::float_format_t(digits)))
{
}

// type conversions

template <int N>
inline
bigfloat<N>::operator double() const
{
	return cln::double_approx(_f);
}
template <int N>
inline
bigfloat<N>::operator int() const
{
	return int(double(*this));
}
#ifdef _RHEOLEF_HAVE_GINAC
template <int N>
inline
bigfloat<N>::operator GiNaC::ex() const
{
	std::stringstream s;
	s << std::setprecision(N) << *this;
	// std::string sx = s.str();
	GiNaC::ex x;
	s >> x;
	return x;
}
#endif // _RHEOLEF_HAVE_GINAC

// assignments

template <int N>
inline
bigfloat<N>&
bigfloat<N>::operator = (const bigfloat<N>& x) 
{
	_f = x._f;
	return *this;
}
#define bigfloat_assignment(T)				\
template <int N>					\
inline							\
bigfloat<N>&						\
bigfloat<N>::operator = (const T& x)			\
{							\
	_f = cln::cl_float(x, cln::float_format_t(digits));	\
	return *this;					\
}
bigfloat_assignment(double)
bigfloat_assignment(int)
bigfloat_assignment(unsigned int)
#undef bigfloat_assignment

template <int N>
inline
bigfloat<N>&
bigfloat<N>::operator = (const char*& x)
{
	_f = cln::cl_float(cln::cl_F(x), cln::float_format_t(digits));
	return *this;
}

// operators

template <int N>
inline
bigfloat<N>
operator - (const bigfloat<N>& x)
{
	return bigfloat<N>(-x._f);
}

# define bigfloat_binary_operator_basic(op)				\
template <int N>							\
inline									\
bigfloat<N>								\
operator op (const bigfloat<N>& x, const bigfloat<N>& y)		\
{									\
	return bigfloat<N>(x._f op y._f);				\
}
# define bigfloat_binary_operator_left(op,T1)				\
template <int N>							\
inline									\
bigfloat<N>								\
operator op (const T1& x, const bigfloat<N>& y)				\
{									\
	return bigfloat<N>(x) op y;					\
}
# define bigfloat_binary_operator_right(op,T2)				\
template <int N>							\
inline									\
bigfloat<N>								\
operator op (const bigfloat<N>& x, const T2& y)				\
{									\
	return x op bigfloat<N>(y);					\
}
# define bigfloat_binary_operator_both(op,T)				\
bigfloat_binary_operator_left(op,T)					\
bigfloat_binary_operator_right(op,T)

# define bigfloat_binary_operator(op) 					\
bigfloat_binary_operator_basic(op)					\
bigfloat_binary_operator_both(op,double)				\
bigfloat_binary_operator_both(op,int)					\
bigfloat_binary_operator_both(op,long)					\
bigfloat_binary_operator_both(op,unsigned int)				\
bigfloat_binary_operator_both(op,unsigned long)

bigfloat_binary_operator(+)
bigfloat_binary_operator(-)
bigfloat_binary_operator(*)
bigfloat_binary_operator(/)

# undef bigfloat_binary_operator_basic
# undef bigfloat_binary_operator

// computed assignments

# define bigfloat_computed_assignment_basic(op,T)		\
template <int N>						\
inline								\
bigfloat<N>&							\
bigfloat<N>::operator op##= (const T& x)			\
{								\
    (*this) = (*this) op x;					\
    return *this;						\
}
# define bigfloat_computed_assignment(op)			\
bigfloat_computed_assignment_basic(op,bigfloat<N>)		\
bigfloat_computed_assignment_basic(op,double)			\
bigfloat_computed_assignment_basic(op,int)			\
bigfloat_computed_assignment_basic(op,long)			\
bigfloat_computed_assignment_basic(op,unsigned int)		\
bigfloat_computed_assignment_basic(op,unsigned long)

bigfloat_computed_assignment(+)
bigfloat_computed_assignment(-)
bigfloat_computed_assignment(*)
bigfloat_computed_assignment(/)

# undef bigfloat_computed_assignment_basic
# undef bigfloat_computed_assignment

// comparators

#define bigfloat_comparator_basic(op)				\
template <int N>						\
inline								\
bool								\
operator op (const bigfloat<N>& x, const bigfloat<N>& y)	\
{								\
	return (cln::compare(x._f,y._f) op 0);			\
}
#define bigfloat_comparator_both(op,T)				\
template <int N>						\
inline								\
bool								\
operator op (const bigfloat<N>& x, const T& y)			\
{								\
	return (x op bigfloat<N>(y));				\
}								\
template <int N>						\
inline								\
bool								\
operator op (const T& x, const bigfloat<N>& y)			\
{								\
	return (bigfloat<N>(x) op y);				\
}
#define bigfloat_comparator(op)					\
  bigfloat_comparator_basic(op)					\
  bigfloat_comparator_both(op,double)				\
  bigfloat_comparator_both(op,int)				\
  bigfloat_comparator_both(op,long)				\
  bigfloat_comparator_both(op,unsigned int)			\
  bigfloat_comparator_both(op,unsigned long)

bigfloat_comparator(==)						\
bigfloat_comparator(!=)						\
bigfloat_comparator(<)						\
bigfloat_comparator(>)						\
bigfloat_comparator(<=)						\
bigfloat_comparator(>=)

#undef bigfloat_comparator_basic
#undef bigfloat_comparator_both
#undef bigfloat_comparator

// math functions

#define bigfloat_math_function(fct, cln_fct)		\
template <int N>					\
inline							\
bigfloat<N>						\
fct (const bigfloat<N>& x)				\
{							\
	return cln_fct (x._f);				\
}
bigfloat_math_function(fabs, abs)
bigfloat_math_function(sqrt, sqrt)
bigfloat_math_function(sin, sin)
bigfloat_math_function(cos, cos)
bigfloat_math_function(tan, tan)
bigfloat_math_function(log, ln)
bigfloat_math_function(exp, exp)
bigfloat_math_function(sinh, sinh)
bigfloat_math_function(cosh, cosh)
bigfloat_math_function(tanh, tanh)
bigfloat_math_function(atan, atan)
bigfloat_math_function(sqr, square)
bigfloat_math_function(floor, ffloor)
bigfloat_math_function(ceil, fceiling)
bigfloat_math_function(trunc, ftruncate)
bigfloat_math_function(round, fround)
#undef bigfloat_math_function

template <int N>
inline
bigfloat<N>
log10 (const bigfloat<N>& x)
{
#ifdef TODO
	// How to store one_on_log10 one time for all ?
	return log(x) * bigfloat<N>::one_on_log10;
#else
	return log(x) / log(bigfloat<N>(10));
#endif // TODO
}
template <int N>
inline
bigfloat<N>
pow (const bigfloat<N>& x, const bigfloat<N>& y)
{
	if (y == trunc(y)) {
	    return pow (x, int(y));
	}
	if (x == bigfloat<N>(0)) return bigfloat<N>(1);
	if (x <  bigfloat<N>(0)) {
	    fatal_macro ("argument x="<<x<<" of pow(x,y) may be >= 0"); 
	}
	return exp(y*log(x));
}
template <int N>
inline
bigfloat<N>
pow (const bigfloat<N>& x, int n)
{
	if (x == bigfloat<N>(0)) {
	    if (n == 0) return bigfloat<N>(1);
	    if (n <  0) return bigfloat<N>(1) / bigfloat<N>(0);
	    else        return bigfloat<N>(0);
	}
	if (x <  bigfloat<N>(0)) {
	    if (n % 2 == 1) {
	        return - exp(n*log(-x));
	    } else {
	        return exp(n*log(-x));
	    }
	}
	return exp(n*log(x));
}
template <int N>
inline
bigfloat<N>
acos (const bigfloat<N>& x)
{
	return bigfloat<N>(cln::cl_float(cln::realpart(acos(x._f))));
}

// inputs/outputs

template<int N>
inline
std::istream&
operator >> (std::istream& s, bigfloat<N>& x)
{
  return all_float_read (s, x);
}
template <int N>
std::ostream&
operator << (std::ostream& s, const bigfloat<N>& x)
{
  return all_float_write (s, x);
}
// for debug
template <int N>
inline
std::ostream&
bigfloat<N>::dump(std::ostream& s)
{
    return s << _f;
}
#endif //_BIGFLOAT_H