/Users/mourrain/Devel/mmx/numerix/include/numerix/floating.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : floating.hpp
* DESCRIPTION: Floating point numbers
* COPYRIGHT  : (C) 2003--2006  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_FLOATING_HPP
#define __MMX_FLOATING_HPP
#include <numerix/rational.hpp>
#include <numerix/string_scnot.hpp>
#include <numerix/mmx_mpfr.hpp>
namespace mmx {
#define TMPL_DEF template<typename V=std_floating>
#define TMPL template<typename V>
#define TMPL2 template<typename V,typename V2>
#define Floating floating<V>
#define Floating2 floating<V2>
#define Floating_rep floating_rep<V>
#define Rnd V::rnd()
TMPL class floating;
TMPL class floating_rep;
TMPL Floating copy (const Floating& x);

/******************************************************************************
* Different variants of floating types
******************************************************************************/

class std_floating {
public:
  static inline mp_rnd_t rnd () {
    return mmx_rounding_mode; }
  static inline xnat prec () {
    return mmx_bit_precision; }
  static inline xnat prec (const mpfr_t& x) {
    //(void) x; return DEFAULT_LOW; }
    return min (mmx_bit_precision, (xnat) mpfr_get_prec (x)); }
  static inline xnat prec (const mpfr_t& x, const mpfr_t& y) {
    //(void) x; (void) y; return DEFAULT_LOW; }
    return min (mmx_bit_precision,
                (xnat) max (mpfr_get_prec (x), mpfr_get_prec (y))); }
  static inline xnat disc () {
    return mmx_discrepancy; }
};

class single_floating: public std_floating {
public:
  static inline xnat prec () {
    return DEFAULT_LOW; }
  static inline xnat prec (const mpfr_t& x) {
    (void) x; return DEFAULT_LOW; }
  static inline xnat prec (const mpfr_t& x, const mpfr_t& y) {
    (void) x; (void) y; return DEFAULT_LOW; }
  static inline xnat disc () {
    return DEFAULT_LOW; }
};

template<typename V, mp_rnd_t R>
class rnd_floating: public V {
public:
  static inline mp_rnd_t rnd () {
    return R; }
};

template<typename V, xnat P, xnat D=P>
class fixed_floating: public V {
public:
  static inline xnat prec () {
    return P; }
  static inline xnat prec (const mpfr_t& x) {
    (void) x; return P; }
  static inline xnat prec (const mpfr_t& x, const mpfr_t& y) {
    (void) x; (void) y; return P; }
  static inline xnat disc () {
    return D; }
};

/******************************************************************************
* The floating type
******************************************************************************/

TMPL_DEF
class floating_rep REP_STRUCT {
private:
  mpfr_t x;

public:
  inline floating_rep () { mpfr_init2 (x, V::prec()); }
  inline floating_rep (xnat prec) { mpfr_init2 (x, prec); }
  inline ~floating_rep () { mpfr_clear (x); }
  friend class Floating;
};

TMPL_DEF
class floating {
INDIRECT_PROTO_1 (floating, floating_rep, V)
public:
  // this should really be private, but we do not want to declare
  // the whole set of inline functions below twice
  inline floating (const mpfr_t& x):
    rep (new Floating_rep (V::prec (x))) {}
  inline floating (const mpfr_t& x, const mpfr_t& y):
    rep (new Floating_rep (V::prec (x, y))) {}

  inline const mpfr_t& operator * () const { return rep->x; }
  inline mpfr_t& operator * () { secure(); return rep->x; }

public:
  // constructors
  inline floating ():
    rep (new Floating_rep ()) {}
  inline floating (signed int i):
    rep (new Floating_rep ()) { mpfr_set_si (rep->x, i, Rnd); }
  inline floating (unsigned int i):
    rep (new Floating_rep ()) { mpfr_set_ui (rep->x, i, Rnd); }
  inline floating (signed short int i):
    rep (new Floating_rep ()) { mpfr_set_si (rep->x, i, Rnd); }
  inline floating (unsigned short int i):
    rep (new Floating_rep ()) { mpfr_set_ui (rep->x, i, Rnd); }
  inline floating (signed long int i):
    rep (new Floating_rep ()) { mpfr_set_si (rep->x, i, Rnd); }
  inline floating (unsigned long int i):
    rep (new Floating_rep ()) { mpfr_set_ui (rep->x, i, Rnd); }
  inline floating (float x2):
    rep (new Floating_rep ()) { mpfr_set_d (rep->x, x2, Rnd); }
  inline floating (double x2):
    rep (new Floating_rep ()) { mpfr_set_d (rep->x, x2, Rnd); }
  inline floating (const integer& i):
    rep (new Floating_rep ()) { mpfr_set_z (rep->x, *i, Rnd); }
  inline floating (const rational& x):
    rep (new Floating_rep ()) { mpfr_set_q (rep->x, *x, Rnd); }
  inline floating (const char* s, int base= 0):
    rep (new Floating_rep ()) { mpfr_set_str (rep->x, s, base, Rnd); }
  inline floating (const string& s, int base= 0):
    rep (new Floating_rep ()) {
      char* r= as_charp (s);
      mpfr_set_str (rep->x, r, base, Rnd);
      free_charp (r); }
};
INDIRECT_IMPL_1 (floating, floating_rep, typename V, V)

template<> struct symbolic_type_information<floating<> > {
  static const nat id= SYMBOLIC_FLOATING; };

TMPL
struct fast_helper<Floating > {
  typedef double fast_type;
  static inline fast_type dd (const Floating& x) { return as_double (x); }
  static inline Floating uu (const fast_type& x) { return Floating (x); }
};

STYPE_TO_TYPE(TMPL,default_radius_type,Floating,floating<single_floating> );

/******************************************************************************
* Standard functions
******************************************************************************/

TMPL inline Floating copy (const Floating& x) {
  Floating r; mpfr_set (*r, *x, Rnd); return r; }
TMPL inline Floating duplicate (const Floating& x) {
  return copy (x); }

TMPL nat hash (const Floating& x) {
  mpfr_t y;
  mpfr_init2 (y, 32);
  mpfr_set (y, *x, GMP_RNDN);
  if (!mpfr_number_p (y)) return 12321;
  if (mpfr_zero_p (y)) return 1234321;
  nat m= (nat) y->_mpfr_d[0];
  nat e= (nat) y->_mpfr_exp;
  nat s= (nat) y->_mpfr_sign;
  mpfr_clear (y);
  return m ^ (m >> ((nat) 31)) ^ e ^ (e << ((nat) 3)) ^ s ^ (s << ((nat) 7));
}

/******************************************************************************
* Conversion functions
******************************************************************************/

TMPL string
as_string (const Floating& x) {
  return mpfr_to_string (*x);
}

TMPL syntactic
flatten (const Floating& g) {
  return flatten_number (as_string (g));
}

TMPL syntactic
flatten (const Floating& g, xnat dd) {
  return flatten_number (as_string (g, dd));
}

TMPL
struct binary_helper<Floating >: public void_binary_helper<Floating > {
  static inline string short_type_name () { return "R"; }
  static inline generic full_type_name () { return "Floating"; }
  static inline generic disassemble (const Floating& x) {
    xnat p= precision (x);
    string s= mpfr_to_string (*x, p + 9);
    return gen_vec (as<generic> (p), as<generic> (s)); }
  static inline Floating assemble (const generic& v) {
    xnat old= mmx_bit_precision;
    mmx_bit_precision= as<xnat> (vector_access (v, 0));
    Floating r= Floating (as<string> (vector_access (v, 1)));
    mmx_bit_precision= old;
    return r; }
  static inline void write (const port& out, const Floating& x) {
    mpfr_binary_write (out, *x); }
  static inline Floating read (const port& in) {
    Floating x; mpfr_binary_read (in, *x); return x; }
};

template<typename W, typename V> floating<W>
as_floating (const floating<V>& x) {
  floating<W> r;
  mpfr_set (*r, *x, W::rnd());
  return r;
}

template<typename W,typename V>
struct as_helper<floating<W>,floating<V> > {
  static inline floating<W> cv (const floating<V>& x) {
    floating<W> r;
    mpfr_set (*r, *x, W::rnd());
    return r;
  }
};

template<typename V> floating<V>
as_floating (const rational& x) {
  floating<V> r;
  mpfr_set_q (*r, *x, V::rnd());
  return r;
}

template<typename V>
struct as_helper<Floating,rational> {
  static inline Floating cv (const rational& x) {
    return as_floating<V> (x);
  }
};

template<typename V>
struct as_helper<double,Floating > {
  static inline double cv (const Floating& x) { return as_double (x); }
};

template<typename F, typename V> inline Floating
promote (const F& x, const Floating&) { return as<Floating > (x); }

TMPL inline int as_int (const Floating& x) {
  return mpfr_get_si (*x, Rnd); }
TMPL inline nat as_nat (const Floating& x) {
  return mpfr_get_ui (*x, Rnd); }
TMPL inline double as_double (const Floating& x) {
  return mpfr_get_d (*x, Rnd); }
TMPL inline integer as_integer (const Floating& x) {
  integer r; mpfr_get_z (*r, *x, Rnd); return r; }

/******************************************************************************
* Constants
******************************************************************************/

TMPL inline void set_nan (Floating& r) {
  mpfr_init (*r); }
TMPL inline void set_infinity (Floating& r) {
  mpfr_set_inf (*r, 1); }
TMPL inline void set_maximal (Floating& r) {
  mpfr_set_inf (*r, 1); }
TMPL inline void set_minimal (Floating& r) {
  mpfr_set_inf (*r, -1); }
TMPL inline void set_largest (Floating& r) {
  mpfr_set_inf (*r, 1); mpfr_nextbelow (*r); }
TMPL inline void set_smallest (Floating& r) {
  mpfr_set_si (*r, 0, GMP_RNDU); mpfr_nextabove (*r); }
TMPL inline void set_accuracy (Floating& r) {
  r= next_above (Floating (1)) - Floating (1); }
TMPL inline void set_log2 (Floating& r) {
  mpfr_const_log2 (*r, Rnd); }
TMPL inline void set_pi (Floating& r) {
  mpfr_const_pi (*r, Rnd); }
TMPL inline void set_euler (Floating& r) {
  mpfr_const_euler (*r, Rnd); }
TMPL inline Floating times_infinity (const Floating& x) {
  return x * infinity_cst<Floating > (); }

/******************************************************************************
* Basic arithmetic
******************************************************************************/

TMPL inline Floating operator - (const Floating& x1) {
  Floating r (*x1); mpfr_neg (*r, *x1, Rnd); return r; }
TMPL inline Floating operator + (const Floating& x1, const Floating& x2) {
  Floating r (*x1, *x2); mpfr_add (*r, *x1, *x2, Rnd); return r; }
TMPL inline Floating operator - (const Floating& x1, const Floating& x2) {
  Floating r (*x1, *x2); mpfr_sub (*r, *x1, *x2, Rnd); return r; }
TMPL inline Floating operator * (const Floating& x1, const Floating& x2) {
  Floating r (*x1, *x2); mpfr_mul (*r, *x1, *x2, Rnd); return r; }
TMPL inline Floating operator / (const Floating& x1, const Floating& x2) {
  Floating r (*x1, *x2); mpfr_div (*r, *x1, *x2, Rnd); return r; }
TMPL inline Floating square (const Floating& x) {
  Floating r (*x); mpfr_mul (*r, *x, *x, Rnd); return r; }
TMPL inline Floating invert (const Floating& x) {
  Floating one (1); Floating r (*x); mpfr_div (*r, *one, *x, Rnd); return r; }

// in-place arithmetic
TMPL inline Floating& operator += (Floating& x1, const Floating& x2) {
  x1.secure (); mpfr_add (*x1, *x1, *x2, Rnd); return x1; }
TMPL inline Floating& operator -= (Floating& x1, const Floating& x2) {
  x1.secure (); mpfr_sub (*x1, *x1, *x2, Rnd); return x1; }
TMPL inline Floating& operator *= (Floating& x1, const Floating& x2) {
  x1.secure (); mpfr_mul (*x1, *x1, *x2, Rnd); return x1; }
TMPL inline Floating& operator /= (Floating& x1, const Floating& x2) {
  x1.secure (); mpfr_div (*x1, *x1, *x2, Rnd); return x1; }

TMPL inline void neg (Floating& r) {
  r.secure (); mpfr_neg (*r, *r, Rnd); }
TMPL inline void neg (Floating& r, const Floating& x1) {
  r.secure (); mpfr_neg (*r, *x1, Rnd); }
TMPL inline void add (Floating& r, const Floating& x1, const Floating& x2) {
  r.secure (); mpfr_add (*r, *x1, *x2, Rnd); }
TMPL inline  void sub (Floating& r, const Floating& x1, const Floating& x2) {
  r.secure (); mpfr_sub (*r, *x1, *x2, Rnd); }
TMPL inline  void mul (Floating& r, const Floating& x1, const Floating& x2) {
  r.secure (); mpfr_mul (*r, *x1, *x2, Rnd); }
TMPL inline  void div (Floating& r, const Floating& x1, const Floating& x2) {
  r.secure (); mpfr_div (*r, *x1, *x2, Rnd); }
TMPL inline void mul_add (Floating& r, const Floating& x, const Floating& y) {
  r.secure (); mpfr_fma (*r, *x, *y, *r, Rnd); }

ARITH_INT_SUGAR(TMPL,Floating);

/******************************************************************************
* Order related functions
******************************************************************************/

TMPL inline bool is_zero (const Floating& x) {
  return mpfr_zero_p (*x) != 0; }

TMPL2 inline bool operator == (const Floating& x1, const Floating2& x2) {
  return mpfr_equal_p (*x1, *x2) != 0; }
TMPL2 inline bool operator != (const Floating& x1, const Floating2& x2) {
  return mpfr_equal_p (*x1, *x2) == 0; }
TMPL2 inline bool operator < (const Floating& x1, const Floating2& x2) {
  return mpfr_less_p (*x1, *x2) != 0; }
TMPL2 inline bool operator > (const Floating& x1, const Floating2& x2) {
  return mpfr_greater_p (*x1, *x2) != 0; }
TMPL2 inline bool operator <= (const Floating& x1, const Floating2& x2) {
  return mpfr_lessequal_p (*x1, *x2) != 0; }
TMPL2 inline bool operator >= (const Floating& x1, const Floating2& x2) {
  return mpfr_greaterequal_p (*x1, *x2) != 0; }

TMPL inline bool operator == (const Floating& x1, const Floating& x2) {
  return mpfr_equal_p (*x1, *x2) != 0; }
TMPL inline bool operator != (const Floating& x1, const Floating& x2) {
  return mpfr_equal_p (*x1, *x2) == 0; }
TMPL inline bool operator < (const Floating& x1, const Floating& x2) {
  return mpfr_less_p (*x1, *x2) != 0; }
TMPL inline bool operator > (const Floating& x1, const Floating& x2) {
  return mpfr_greater_p (*x1, *x2) != 0; }
TMPL inline bool operator <= (const Floating& x1, const Floating& x2) {
  return mpfr_lessequal_p (*x1, *x2) != 0; }
TMPL inline bool operator >= (const Floating& x1, const Floating& x2) {
  return mpfr_greaterequal_p (*x1, *x2) != 0; }

TMPL inline int sign (const Floating& x1) {
  int r= mpfr_sgn (*x1); return (r<0? -1: (r==0? 0: 1)); }
TMPL inline Floating abs (const Floating& x1) {
  Floating r (*x1); mpfr_abs (*r, *x1, Rnd); return r; }
TMPL inline Floating min (const Floating& x1, const Floating& x2) {
  return x1 <= x2? x1: x2; }
TMPL inline Floating max (const Floating& x1, const Floating& x2) {
  return x1 >= x2? x1: x2; }

TMPL inline Floating floor (const Floating& x1) {
  Floating r (*x1); mpfr_floor (*r, *x1); return r; }
TMPL inline Floating trunc (const Floating& x1) {
  Floating r (*x1); mpfr_trunc (*r, *x1); return r; }  
TMPL inline Floating ceil (const Floating& x1) {
  Floating r (*x1); mpfr_ceil (*r, *x1); return r; }
TMPL inline Floating round (const Floating& x1) {
  Floating r (*x1); mpfr_round (*r, *x1); return r; }

TRUE_TO_EXACT_IDENTITY_SUGAR(TMPL,Floating)
EQUAL_INT_SUGAR(TMPL,Floating);
COMPARE_INT_SUGAR(TMPL,Floating);

/******************************************************************************
* Elementary functions
******************************************************************************/

TMPL inline Floating sqrt (const Floating& x1) {
  Floating r (*x1); mpfr_sqrt (*r, *x1, Rnd); return r; }
TMPL inline Floating cbrt (const Floating& x1) {
  Floating r (*x1); mpfr_cbrt (*r, *x1, Rnd); return r; }
TMPL inline Floating hypot (const Floating& x1, const Floating& x2) {
  Floating r (*x1, *x2); mpfr_hypot (*r, *x1, *x2, Rnd); return r; }
TMPL inline Floating exp (const Floating& x1) {
  Floating r (*x1); mpfr_exp (*r, *x1, Rnd); return r; }
TMPL inline Floating exp2 (const Floating& x1) {
  Floating r (*x1); mpfr_exp2 (*r, *x1, Rnd); return r; }
TMPL inline Floating log (const Floating& x1) {
  Floating r (*x1); mpfr_log (*r, *x1, Rnd); return r; }
TMPL inline Floating log2 (const Floating& x1) {
  Floating r (*x1); mpfr_log2 (*r, *x1, Rnd); return r; }
TMPL inline Floating log10 (const Floating& x1) {
  Floating r (*x1); mpfr_log10 (*r, *x1, Rnd); return r; }
TMPL inline Floating cos (const Floating& x1) {
  Floating r (*x1); mpfr_cos (*r, *x1, Rnd); return r; }
TMPL inline Floating sin (const Floating& x1) {
  Floating r (*x1); mpfr_sin (*r, *x1, Rnd); return r; }
TMPL inline Floating tan (const Floating& x1) {
  Floating r (*x1); mpfr_tan (*r, *x1, Rnd); return r; }
TMPL inline Floating acos (const Floating& x1) {
  Floating r (*x1); mpfr_acos (*r, *x1, Rnd); return r; }
TMPL inline Floating asin (const Floating& x1) {
  Floating r (*x1); mpfr_asin (*r, *x1, Rnd); return r; }
TMPL inline Floating atan (const Floating& x1) {
  Floating r (*x1); mpfr_atan (*r, *x1, Rnd); return r; }
//TMPL inline Floating atan2 (const Floating& x1, const Floating& x2) {
//Floating r (*x1, *x2); mpfr_atan2 (*r, *x1, *x2, Rnd); return r; }
TMPL inline Floating cosh (const Floating& x1) {
  Floating r (*x1); mpfr_cosh (*r, *x1, Rnd); return r; }
TMPL inline Floating sinh (const Floating& x1) {
  Floating r (*x1); mpfr_sinh (*r, *x1, Rnd); return r; }
TMPL inline Floating tanh (const Floating& x1) {
  Floating r (*x1); mpfr_tanh (*r, *x1, Rnd); return r; }
TMPL inline Floating acosh (const Floating& x1) {
  Floating r (*x1); mpfr_acosh (*r, *x1, Rnd); return r; }
TMPL inline Floating asinh (const Floating& x1) {
  Floating r (*x1); mpfr_asinh (*r, *x1, Rnd); return r; }
TMPL inline Floating atanh (const Floating& x1) {
  Floating r (*x1); mpfr_atanh (*r, *x1, Rnd); return r; }
TMPL inline Floating pow (const Floating& x1, const Floating& x2) {
  Floating r (*x1, *x2); mpfr_pow (*r, *x1, *x2, Rnd); return r; }
TMPL inline Floating pow (const Floating& x1, const int& x2) {
  Floating r (*x1); mpfr_pow_ui (*r, *x1, x2, Rnd); return r; }
TMPL inline Floating pow (const int& x1, const Floating& x2) {
  Floating r (*x2); mpfr_ui_pow (*r, x1, *x2, Rnd); return r; }

ATAN2_SUGAR(TMPL,Floating)
INV_TRIGO_SUGAR(TMPL,Floating)
INV_HYPER_SUGAR(TMPL,Floating)

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Floating& x) {
  return mpfr_number_p (*x) != 0; }
TMPL inline bool is_infinite (const Floating& x) {
  return mpfr_inf_p (*x) != 0; }
TMPL inline bool is_fuzz (const Floating& x) {
  (void) x; return false; }
TMPL inline bool is_nan (const Floating& x) {
  return mpfr_nan_p (*x) != 0; }
TMPL inline bool is_reliable (const Floating& x) {
  (void) x; return false; }

TMPL inline Floating change_precision (const Floating& x, xnat prec) {
  Floating r (new Floating_rep (prec)); mpfr_set (*r, *x, Rnd); return r; }
TMPL inline xnat precision (const Floating& x) {
  return (xnat) mpfr_get_prec (*x); }
TMPL inline Floating next_above (const Floating& x1) {
  Floating r (*x1); mpfr_set (*r, *x1, Rnd);
  mpfr_nextabove (*r); return r; }
TMPL inline Floating next_below (const Floating& x1) {
  Floating r (*x1); mpfr_set (*r, *x1, Rnd);
  mpfr_nextbelow (*r); return r; }
TMPL inline Floating rounding_error (const Floating& x) {
  Floating r (*x);
  mpfr_abs (*r, *x, GMP_RNDN);
  mpfr_mul_2si (*r, *r, 1-mpfr_get_prec (*x), GMP_RNDN);
  mpfr_nextabove (*r);
  return r; }
TMPL inline Floating additive_error (const Floating& x) {
  Floating r (*x);
  mpfr_abs (*r, *x, GMP_RNDN);
  mpfr_mul_2si (*r, *r, 1-mpfr_get_prec (*x), GMP_RNDN);
  mpfr_nextabove (*r);
  return r; }

TMPL inline xint exponent (const Floating& x) {
  return mpfr_safe_get_exp (*x); }
TMPL inline double magnitude (const Floating& x) {
  return mpfr_get_magnitude (*x); }
TMPL inline Floating operator << (const Floating& x1, const xint& shift) {
  Floating r (*x1); mpfr_mul_2si (*r, *x1, shift, Rnd); return r; }
TMPL inline Floating operator >> (const Floating& x1, const xint& shift) {
  Floating r (*x1); mpfr_mul_2si (*r, *x1, -shift, Rnd); return r; }
TMPL inline Floating& operator <<= (Floating& x1, const xint& shift) {
  return x1= x1 << shift; }
TMPL inline Floating& operator >>= (Floating& x1, const xint& shift) {
  return x1= x1 >> shift; }

/******************************************************************************
* Special functions and other functions
******************************************************************************/

TMPL inline Floating gamma (const Floating& x1) {
  Floating r (*x1); mpfr_gamma (*r, *x1, Rnd); return r; }
TMPL inline Floating zeta (const Floating& x1) {
  Floating r (*x1); mpfr_zeta (*r, *x1, Rnd); return r; }
TMPL inline Floating erf (const Floating& x1) {
  Floating r (*x1); mpfr_erf (*r, *x1, Rnd); return r; }

TMPL inline Floating uniform_deviate (const Floating& lo, const Floating& hi) {
  Floating r (*lo, *hi); mpfr_urandomb (*r, gmp_randstate);
  return lo + r * (hi - lo); }

/******************************************************************************
* Extras for glue
******************************************************************************/

#define mmx_floating floating<>
#define make_literal_floating make_literal<floating<> >

#undef TMPL_DEF
#undef TMPL
#undef TMPL2
#undef Floating
#undef Floating2
#undef Floating_rep
#undef Rnd
}
#endif // __MMX_FLOATING_HPP