#include "mpfr_real.h"

//MpfrReal::PrecisionType &MpfrReal::CurrPrecision = MpfrReal::PrecisionType(__gmp_default_fp_bit_precision);
MpfrReal::PrecisionType &MpfrReal::CurrPrecision = __gmp_default_fp_bit_precision;
RoundingModes &MpfrReal::CurrRndMode = __gmp_default_rounding_mode;



//--------------------------------------------------------------
//
// Constructors and destructors
//
//--------------------------------------------------------------

// inline MpfrReal::MpfrReal ()
MpfrReal::MpfrReal ()
  {
  mpfr_init2(mpfr_rep, CurrPrecision);
  nbref = new RefCounter(0);
  }

// inline MpfrReal::MpfrReal (double d, 
MpfrReal::MpfrReal (double d, 
                    RoundingModes rnd, 
                    PrecisionType prec)
  {
  mpfr_init2(mpfr_rep, prec);
  mpfr_set_d(mpfr_rep, d, rnd);
  nbref=new RefCounter(1);
  }

// inline MpfrReal::MpfrReal (int i, 
MpfrReal::MpfrReal (int i, 
                    RoundingModes rnd, 
                    PrecisionType prec)
  {
  mpfr_init2(mpfr_rep, prec);
  mpfr_set_si(mpfr_rep, long(i), rnd);
  nbref=new RefCounter(1);
  }

// inline MpfrReal::MpfrReal (long int i, 
MpfrReal::MpfrReal (long int i, 
                    RoundingModes rnd, 
                    PrecisionType prec)
  {
  mpfr_init2(mpfr_rep, prec);
  mpfr_set_si(mpfr_rep, i, rnd);
  nbref=new RefCounter(1);
  }

// inline MpfrReal::MpfrReal (unsigned long int i,
MpfrReal::MpfrReal (unsigned long int i,
                    RoundingModes rnd, 
                    PrecisionType prec)
  {
  mpfr_init2(mpfr_rep, prec);
  mpfr_set_ui(mpfr_rep, i, rnd);
  nbref=new RefCounter(1);
  }

MpfrReal::MpfrReal(string s)
  {
  mpfr_init2(mpfr_rep,MpfrReal::CurrPrecision);
  if ( mpfr_set_str(mpfr_rep, const_cast<char *>(s.c_str()), 10, MpfrReal::CurrRndMode) )
    {
    cerr << "Pb while reading " << s << endl ;
    // maybe throw an exception...
    }
  nbref=new RefCounter(1);
  }

// inline MpfrReal::MpfrReal (const MpfrReal& r)
MpfrReal::MpfrReal (const MpfrReal& r)
  {
  mpfr_rep[0]=r.mpfr_rep[0];
  nbref=r.nbref;
  nbref->incr();
  }
  
// inline MpfrReal::~MpfrReal ()
MpfrReal::~MpfrReal ()
  {
  if (nbref->decr() <= 0)
    mpfr_clear(mpfr_rep);
  }


//--------------------------------------------------------------
//
// Assignment and physical copy
//
//--------------------------------------------------------------

// Assignment is a logical copy
// inline MpfrReal& MpfrReal::operator = (const MpfrReal& r)
MpfrReal& MpfrReal::operator = (const MpfrReal& r)
  {
  if (this == &r)
    return *this;
  else if (nbref->decr() <= 0)
    {
    mpfr_clear(mpfr_rep);
    }
  mpfr_rep[0]=r.mpfr_rep[0];
  nbref=r.nbref;
  nbref->incr();
  return *this;
  }
  
// Physical copy
MpfrReal& MpfrReal::copy (const MpfrReal& r, RoundingModes rnd, PrecisionType prec)
  {
  if (nbref->decr() <= 0)
    {
    mpfr_clear(mpfr_rep);
    }
  // If the desired precision is different from the current one
  if ( prec != GetDefaultPrecision() )
    {
    mpfr_init2(mpfr_rep, prec);
    mpfr_set(mpfr_rep, r.mpfr_rep, rnd);
    }
  else
    mpfr_init_set(mpfr_rep, r.mpfr_rep, rnd);

  nbref=new RefCounter(1);
  return *this;
  }
  

//--------------------------------------------------------------
//
// Precision and rounding: consultation, modification
//
//--------------------------------------------------------------

// inline void MpfrReal::SetDefaultPrecision (MpfrReal::PrecisionType newprec)
void MpfrReal::SetDefaultPrecision (MpfrReal::PrecisionType newprec)
  {
  mpfr_set_default_prec(newprec);
  }

// inline void MpfrReal::SetPrecision (MpfrReal::PrecisionType newprec)
void MpfrReal::SetPrecision (MpfrReal::PrecisionType newprec)
  {
  mpfr_set_prec(mpfr_rep, newprec);
  }

// inline const MpfrReal::PrecisionType MpfrReal::GetDefaultPrecision ()
const MpfrReal::PrecisionType MpfrReal::GetDefaultPrecision ()
  {
  return CurrPrecision;
  }

// inline const MpfrReal::PrecisionType MpfrReal::GetPrecision ()
const MpfrReal::PrecisionType MpfrReal::GetPrecision ()
  {
  return mpfr_get_prec(mpfr_rep);
  }

// inline void MpfrReal::SetDefaultRndMode (RoundingModes newrndmode)
void MpfrReal::SetDefaultRndMode (RoundingModes newrndmode)
  {
  mpfr_set_default_rounding_mode(newrndmode);
  }

// inline const RoundingModes MpfrReal::GetDefaultRndMode ()
const RoundingModes MpfrReal::GetDefaultRndMode ()
  {
  return CurrRndMode;
  }

// inline void MpfrReal::Round (RoundingModes rnd, PrecisionType prec)
void MpfrReal::Round (RoundingModes rnd, PrecisionType prec)
  {
  mpfr_round(mpfr_rep, rnd, prec);
  }

//--------------------------------------------------------------
//
// Constants: Pi and log(2)
//
//--------------------------------------------------------------

MpfrReal& MpfrReal::Pi (RoundingModes rnd, PrecisionType prec) 
  {
  MpfrReal& res= * new MpfrReal;
  mpfr_init2(res.mpfr_rep, prec);
  mpfr_const_pi(res.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
  }

MpfrReal& MpfrReal::Log2 (RoundingModes rnd, PrecisionType prec) 
  {
  MpfrReal& res= * new MpfrReal;
  mpfr_init2(res.mpfr_rep, prec);
  mpfr_const_log2(res.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
}
    
  
//--------------------------------------------------------------
//
// Input Output
//
//--------------------------------------------------------------

ostream& MpfrReal::put(ostream& o, RoundingModes rnd, PrecisionType prec) const
  {
  if (MPFR_IS_NAN(mpfr_rep))
    { o << "NaN"; return o; }
  if ( iszero(*this) )
    { o << "0"; return o; }

  // Mimicking mpfr_get_str, especially for the computation of the allocated length of s
  char *s, *fraction, first_digit;
  mp_exp_t e;
  int size_s;
  s = mpfr_get_str(NULL, &e, 10, int(prec), mpfr_rep,rnd);
  size_s = int( double(MPFR_PREC(mpfr_rep))*log(2.0)/log(10.0) ) + 2;
  // to print the floating point
  if (s[0] == '-')
    {
    o << s[0];
    first_digit = s[1];
    fraction = &s[2];
    size_s++;
    }
  else
    {
    first_digit = s[0];
    fraction = &s[1];
    }
  o << first_digit << "." << fraction;
  (*_mp_free_func)(s, size_s);
  if (--e)
    {
    o << "E" << e ;
    }
  return o;
  }
  
ostream& operator << (ostream& o, const MpfrReal& r)
  {
  return r.put(o);
  }

istream& operator >> (istream& i, MpfrReal& r)
// Mimicking mpfr_set_str and Integer::>> of Givaro
  {
  if (!i.good())
    { 
    cerr << " Pb reading on the input stream " << endl;
    // maybe throw an exception...
    }
  if (!i)
    return i;

  // read white spaces
  i >> ws;

  // reading the string, the rounding is done by mpfr_set_str
  int noend = 1, nread=0;
  size_t alloc_size=100;
  char c, *str=(char *)(*_mp_allocate_func)(alloc_size);

  // read the characters on i until a white space is read
  while (noend)
    {
    // read until alloc_size char are read, or a white space is encountered
    while ( (noend) && (nread < alloc_size) )
      {
      i.get(c);
      if (i.eof())
        { noend=0; }
      else if (isspace(c))
        { noend=0; i.putback(c); }
      else
        { str[nread++] = c; }
      }

    if (nread >= alloc_size)
      {
      size_t old_alloc_size = alloc_size;
      alloc_size = alloc_size * 3 / 2;
      str = (char *)(*_mp_reallocate_func)(str, old_alloc_size, alloc_size);
      }
    }

  str[nread] = '\0';

  if ( mpfr_set_str(r.mpfr_rep, str, 10, MpfrReal::CurrRndMode) )
    {
    cerr << " Pb reading on the input stream " << endl;
    // maybe throw an exception...
    }

  return i;
  }


//--------------------------------------------------------------
//
// Comparisons
//
//--------------------------------------------------------------

int compare (const MpfrReal& r1, const MpfrReal& r2)
  { return mpfr_cmp(r1.mpfr_rep, r2.mpfr_rep); }

// inline int compare (const MpfrReal& r1, const int r2)
int compare (const MpfrReal& r1, const int r2)
  { return mpfr_cmp_si(r1.mpfr_rep, long(r2) ); }

// inline int compare (const MpfrReal& r1, const unsigned int r2)
int compare (const MpfrReal& r1, const unsigned int r2)
  { return mpfr_cmp_ui(r1.mpfr_rep, (unsigned long) r2) ; }

// inline int compare (const MpfrReal& r1, const long int r2)
int compare (const MpfrReal& r1, const long int r2)
  { return mpfr_cmp_si(r1.mpfr_rep, r2); }

// inline int compare (const MpfrReal& r1, const unsigned long r2)
int compare (const MpfrReal& r1, const unsigned long r2)
  { return mpfr_cmp_ui(r1.mpfr_rep, r2); }

// inline int operator == (const MpfrReal& r1, const MpfrReal& r2)
int operator == (const MpfrReal& r1, const MpfrReal& r2)
  { return (compare(r1,r2) == 0); }

// inline int operator == (const int r1, const MpfrReal& r2)
int operator == (const int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator == (const unsigned int r1, const MpfrReal& r2)
int operator == (const unsigned int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator == (const long int r1, const MpfrReal& r2)
int operator == (const long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator == (const unsigned long int r1, const MpfrReal& r2)
int operator == (const unsigned long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator == (const MpfrReal& r1, const int r2)
int operator == (const MpfrReal& r1, const int r2)
  { return (compare(r1,r2) == 0); }

// inline int operator == (const MpfrReal& r1, const unsigned int r2)
int operator == (const MpfrReal& r1, const unsigned int r2)
  { return (compare(r1,r2) == 0); }

// inline int operator == (const MpfrReal& r1, const long int r2)
int operator == (const MpfrReal& r1, const long int r2)
  { return (compare(r1,r2) == 0); }

// inline int operator == (const MpfrReal& r1, const unsigned long int r2)
int operator == (const MpfrReal& r1, const unsigned long int r2)
  { return (compare(r1,r2) == 0); }

// inline int operator != (const MpfrReal& r1, const MpfrReal& r2)
int operator != (const MpfrReal& r1, const MpfrReal& r2)
  { return (compare(r1,r2) != 0); }

// inline int operator != (const int r1, const MpfrReal& r2)
int operator != (const int r1, const MpfrReal& r2)
  { return (compare(r2,r1) != 0); }

// inline int operator != (const unsigned int r1, const MpfrReal& r2)
int operator != (const unsigned int r1, const MpfrReal& r2)
  { return (compare(r2,r1) != 0); }

// inline int operator != (const long int r1, const MpfrReal& r2)
int operator != (const long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) != 0); }

// inline int operator != (const unsigned long int r1, const MpfrReal& r2)
int operator != (const unsigned long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) != 0); }

// inline int operator != (const MpfrReal& r1, const int r2)
int operator != (const MpfrReal& r1, const int r2)
  { return (compare(r1,r2) != 0); }

// inline int operator != (const MpfrReal& r1, const unsigned int r2)
int operator != (const MpfrReal& r1, const unsigned int r2)
  { return (compare(r1,r2) != 0); }

// inline int operator != (const MpfrReal& r1, const long int r2)
int operator != (const MpfrReal& r1, const long int r2)
  { return (compare(r1,r2) != 0); }

// inline int operator != (const MpfrReal& r1, const unsigned long int r2)
int operator != (const MpfrReal& r1, const unsigned long int r2)
  { return (compare(r1,r2) != 0); }

// inline int operator <  (const MpfrReal& r1, const MpfrReal& r2)
int operator <  (const MpfrReal& r1, const MpfrReal& r2)
  { return (compare(r1,r2) < 0); }

// inline int operator <  (const int r1, const MpfrReal& r2)
int operator <  (const int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >  0); }

// inline int operator <  (const unsigned int r1, const MpfrReal& r2)
int operator <  (const unsigned int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >  0); }

// inline int operator <  (const long int r1, const MpfrReal& r2)
int operator <  (const long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >  0); }

// inline int operator <  (const unsigned long int r1, const MpfrReal& r2)
int operator <  (const unsigned long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >  0); }

// inline int operator <  (const MpfrReal& r1, const int r2)
int operator <  (const MpfrReal& r1, const int r2)
  { return (compare(r1,r2) <  0); }

// inline int operator <  (const MpfrReal& r1, const unsigned int r2)
int operator <  (const MpfrReal& r1, const unsigned int r2)
  { return (compare(r1,r2) <  0); }

// inline int operator <  (const MpfrReal& r1, const long int r2)
int operator <  (const MpfrReal& r1, const long int r2)
  { return (compare(r1,r2) <  0); }

// inline int operator <  (const MpfrReal& r1, const unsigned long int r2)
int operator <  (const MpfrReal& r1, const unsigned long int r2)
  { return (compare(r1,r2) <  0); }

// inline int operator <= (const MpfrReal& r1, const MpfrReal& r2)
int operator <= (const MpfrReal& r1, const MpfrReal& r2)
  { return (compare(r1,r2) <= 0); }

// inline int operator <= (const int r1, const MpfrReal& r2)
int operator <= (const int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >= 0); }

// inline int operator <= (const unsigned int r1, const MpfrReal& r2)
int operator <= (const unsigned int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >= 0); }

// inline int operator <= (const long int r1, const MpfrReal& r2)
int operator <= (const long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >= 0); }

// inline int operator <= (const unsigned long int r1, const MpfrReal& r2)
int operator <= (const unsigned long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) >= 0); }

// inline int operator <= (const MpfrReal& r1, const int r2)
int operator <= (const MpfrReal& r1, const int r2)
  { return (compare(r1,r2) <= 0); }

// inline int operator <= (const MpfrReal& r1, const unsigned int r2)
int operator <= (const MpfrReal& r1, const unsigned int r2)
  { return (compare(r1,r2) <= 0); }

// inline int operator <= (const MpfrReal& r1, const long int r2)
int operator <= (const MpfrReal& r1, const long int r2)
  { return (compare(r1,r2) <= 0); }

// inline int operator <= (const MpfrReal& r1, const unsigned long int r2)
int operator <= (const MpfrReal& r1, const unsigned long int r2)
  { return (compare(r1,r2) <= 0); }

// inline int operator >  (const MpfrReal& r1, const MpfrReal& r2)
int operator >  (const MpfrReal& r1, const MpfrReal& r2)
  { return (compare(r1,r2) > 0); }

// inline int operator >  (const int r1, const MpfrReal& r2)
int operator >  (const int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator >  (const unsigned int r1, const MpfrReal& r2)
int operator >  (const unsigned int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator >  (const long int r1, const MpfrReal& r2)
int operator >  (const long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator >  (const unsigned long int r1, const MpfrReal& r2)
int operator >  (const unsigned long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) == 0); }

// inline int operator >  (const MpfrReal& r1, const int r2)
int operator >  (const MpfrReal& r1, const int r2)
  { return (compare(r1,r2) >  0); }

// inline int operator >  (const MpfrReal& r1, const unsigned int r2)
int operator >  (const MpfrReal& r1, const unsigned int r2)
  { return (compare(r1,r2) >  0); }

// inline int operator >  (const MpfrReal& r1, const long int r2)
int operator >  (const MpfrReal& r1, const long int r2)
  { return (compare(r1,r2) >  0); }

// inline int operator >  (const MpfrReal& r1, const unsigned long int r2)
int operator >  (const MpfrReal& r1, const unsigned long int r2)
  { return (compare(r1,r2) >  0); }

// inline int operator >= (const MpfrReal& r1, const MpfrReal& r2)
int operator >= (const MpfrReal& r1, const MpfrReal& r2)
  { return (compare(r1,r2) >= 0); }

// inline int operator >= (const int r1, const MpfrReal& r2)
int operator >= (const int r1, const MpfrReal& r2)
  { return (compare(r2,r1) <= 0); }

// inline int operator >= (const unsigned int r1, const MpfrReal& r2)
int operator >= (const unsigned int r1, const MpfrReal& r2)
  { return (compare(r2,r1) <= 0); }

// inline int operator >= (const long int r1, const MpfrReal& r2)
int operator >= (const long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) <= 0); }

// inline int operator >= (const unsigned long int r1, const MpfrReal& r2)
int operator >= (const unsigned long int r1, const MpfrReal& r2)
  { return (compare(r2,r1) <= 0); }

// inline int operator >= (const MpfrReal& r1, const int r2)
int operator >= (const MpfrReal& r1, const int r2)
  { return (compare(r1,r2) >= 0); }

// inline int operator >= (const MpfrReal& r1, const unsigned int r2)
int operator >= (const MpfrReal& r1, const unsigned int r2)
  { return (compare(r1,r2) >= 0); }

// inline int operator >= (const MpfrReal& r1, const long int r2)
int operator >= (const MpfrReal& r1, const long int r2)
  { return (compare(r1,r2) >= 0); }

// inline int operator >= (const MpfrReal& r1, const unsigned long int r2)
int operator >= (const MpfrReal& r1, const unsigned long int r2)
  { return (compare(r1,r2) >= 0); }

// inline const MpfrReal& min (const MpfrReal& a, MpfrReal& b)
const MpfrReal& min (const MpfrReal& a, MpfrReal& b)
  {
  if (a <= b)
    return a;
  else
    return b;
  }

// inline const MpfrReal& max (const MpfrReal& a, MpfrReal& b)
const MpfrReal& max (const MpfrReal& a, MpfrReal& b)
  {
  if (a >= b)
    return a;
  else
    return b;
  }



//--------------------------------------------------------------
//
// Arithmetic operations
//
//--------------------------------------------------------------

// Addition-----------------------------------------------------
MpfrReal& MpfrReal::operator + (const MpfrReal& b) const 
  {
  MpfrReal& res = * new MpfrReal;

  if ( iszero(*this) )
   { res = b; return res; }
  else if ( iszero(b) )
    { res = *this; return res;}
  else
    {
    res = *this;
    res +=b;
    return res;
    }
  }

/*
MpfrReal& operator + (const MpfrReal& a, const MpfrReal& b) 
  {
  MpfrReal& res = * new MpfrReal;

  if ( iszero(b) )
    res = a;
  else if ( iszero(a) )
    res = b;
  else
    {
    res = a;
    res +=b;
    }
  return res;
  }
*/

// inline MpfrReal& operator + (const MpfrReal& a, const double b) 
MpfrReal& operator + (const MpfrReal& a, const double b) 
  {
  return a + MpfrReal(b);
  }

// inline MpfrReal& operator + (const MpfrReal& a, const int b) 
MpfrReal& operator + (const MpfrReal& a, const int b) 
  {
  MpfrReal& res = * new MpfrReal;
  if (b == 0)
    res = a;
  else if ( iszero(a) )
    res = b;
  else
    { res = a; res +=b; }
  return res;
  }

// inline MpfrReal& operator + (const MpfrReal& a, const unsigned int b) 
MpfrReal& operator + (const MpfrReal& a, const unsigned int b) 
  {
  MpfrReal& res = * new MpfrReal;
  if (b == 0)
    res = a;
  else if ( iszero(a) )
    res = MpfrReal((unsigned long int) b);
  else
    { res = a; res +=b; }
  return res;
  }

// inline MpfrReal& operator + (const MpfrReal& a, const long int b) 
MpfrReal& operator + (const MpfrReal& a, const long int b) 
  {
  MpfrReal& res = * new MpfrReal;
  if (b == 0)
    res = a;
  else if ( iszero(a) )
    res = MpfrReal(b);
  else
    { res = a; res +=b; }
  return res;
  }

// inline MpfrReal& operator + (const MpfrReal& a, const unsigned long int b) 
MpfrReal& operator + (const MpfrReal& a, const unsigned long int b) 
  {
  MpfrReal& res = * new MpfrReal;
  if (b == 0)
    res = a;
  else if ( iszero(a) )
    res = MpfrReal(b);
  else
    { res = a; res +=b; }
  return res;
  }

// inline MpfrReal& operator + (const double a, const MpfrReal& b)
MpfrReal& operator + (const double a, const MpfrReal& b)
  { return MpfrReal(a)+b; }

// inline MpfrReal& operator + (const int a, const MpfrReal& b)
MpfrReal& operator + (const int a, const MpfrReal& b)
  { return b+a; }

// inline MpfrReal& operator + (const unsigned int a, const MpfrReal& b)
MpfrReal& operator + (const unsigned int a, const MpfrReal& b)
  { return b+a; }

// inline MpfrReal& operator + (const long int a, const MpfrReal& b)
MpfrReal& operator + (const long int a, const MpfrReal& b)
  { return b+a; }

// inline MpfrReal& operator + (const unsigned long int a, const MpfrReal& b)
MpfrReal& operator + (const unsigned long int a, const MpfrReal& b)
  { return b+a; }

MpfrReal& MpfrReal::operator += (const MpfrReal& b) 
  {
  if ( iszero(b) )
    return *this;
  if ( iszero(*this) )
    {
    if ( nbref->decr() <= 0 )			// no other ref. on the value of *this
      mpfr_clear(mpfr_rep);
    mpfr_rep[0] = b.mpfr_rep[0];
    nbref = b.nbref;
    nbref->incr();
    return *this;
    }
  if ( nbref->decr() <= 0 )			// no other ref. on the value of *this
    {						// the memory can be reused
    mpfr_add(mpfr_rep, mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else						// the value and memory used must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    mpfr_add(mpfr_rep, tmp.mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator += (const double b) 
MpfrReal& MpfrReal::operator += (const double b) 
  {
  return (*this) += MpfrReal(b);
  }

MpfrReal& MpfrReal::operator += (const long int b) 	// using the more efficient mpfr_add_ui
  {
  if ( b == 0)
    return *this;
  if ( iszero(*this) )
    {
    *this = MpfrReal(b);
    }
  if ( nbref->decr() <= 0 )			// the memory can be reused
    {
    if (b > 0)
      mpfr_add_ui(mpfr_rep, mpfr_rep, (unsigned long int)b, CurrRndMode);
    else
      mpfr_sub_ui(mpfr_rep, mpfr_rep, (unsigned long int)(-b), CurrRndMode);
    nbref->refvalue() = 1;
    }
  else						// the previous value and memory must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    if (b > 0)
      mpfr_add_ui(mpfr_rep, tmp.mpfr_rep, (unsigned long int)b, CurrRndMode);
    else
      mpfr_sub_ui(mpfr_rep, tmp.mpfr_rep, (unsigned long int)(-b), CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

MpfrReal& MpfrReal::operator += (const unsigned long int b) 	// using mpfr_add_ui
  {
  if ( b == 0)
    return *this;
  if ( iszero(*this) )
    {
    *this = MpfrReal(b);
    }
  if ( nbref->decr() <= 0 )			// the memory can be reused
    {
    mpfr_add_ui(mpfr_rep, mpfr_rep, (unsigned long int)b, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else
    {
    MpfrReal tmp = *this;			// the previous value and memory must be preserved
    mpfr_init(mpfr_rep);
    mpfr_add_ui(mpfr_rep, tmp.mpfr_rep, (unsigned long int)b, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator += (const int b) 	// using the efficient mpfr_add_ui
MpfrReal& MpfrReal::operator += (const int b) 	// using the efficient mpfr_add_ui
  {
  return (*this)+= long(b);
  }

// inline MpfrReal& MpfrReal::operator += (const unsigned int b) // using mpfr_add_ui
MpfrReal& MpfrReal::operator += (const unsigned int b) // using mpfr_add_ui
  {
  return (*this)+= (unsigned long int)(b);
  }

MpfrReal& MpfrReal::add (MpfrReal& res, const MpfrReal& r1, const MpfrReal& r2, RoundingModes rnd)
  {
  if ( iszero(r1) )
    return res = r2;
  else if ( iszero(r2) )
    {
    return res = r1;
    }
      
  if ( res.nbref->decr() <= 0 )			// the memory can be reused
    {
    mpfr_add(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref->refvalue() = 1;
    }
  else						// the previous memory must be preserved
    {
    mpfr_init(res.mpfr_rep);
    mpfr_add(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref = new RefCounter(1);
    }
  return res;
  }
  
// Subtraction --------------------------------------------------------
MpfrReal& MpfrReal::operator - (const MpfrReal& b) const
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);                     
  mpfr_sub(res.mpfr_rep, mpfr_rep, b.mpfr_rep, CurrRndMode);  
  res.nbref = new RefCounter(1);
  return res;
  }

/*
MpfrReal& operator - (const MpfrReal& a, const MpfrReal& b) 
  {
  MpfrReal& res = * new MpfrReal;
  if ( iszero(b) )
    res = a;
  else if ( iszero(a) )
    res = -b;
  else
    { res = a; res -= b; }
  return res;
  }
*/

// inline MpfrReal& operator - (const MpfrReal& a, const double b) 
MpfrReal& operator - (const MpfrReal& a, const double b) 
  {
  return a - MpfrReal(b);
  }

// inline MpfrReal& operator - (const MpfrReal &a, const int b)
MpfrReal& operator - (const MpfrReal &a, const int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res -= b;
  return res;
  }

// inline MpfrReal& operator - (const MpfrReal &a, const unsigned int b)
MpfrReal& operator - (const MpfrReal &a, const unsigned int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res -= b;
  return res;
  }

// inline MpfrReal& operator - (const MpfrReal &a, const long int b)
MpfrReal& operator - (const MpfrReal &a, const long int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res -= b;
  return res;
  }

// inline MpfrReal& operator - (const MpfrReal &a, const unsigned long int b)
MpfrReal& operator - (const MpfrReal &a, const unsigned long int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res -= b;
  return res;
  }

// inline MpfrReal& operator - (const double a, const MpfrReal&  b)
MpfrReal& operator - (const double a, const MpfrReal&  b)
  {
  return MpfrReal(a) - b;
  }

// inline MpfrReal& operator - (const int a, const MpfrReal&  b)
MpfrReal& operator - (const int a, const MpfrReal&  b)
  {
  return -(b-a);
  }

// inline MpfrReal& operator - (const unsigned int a, const MpfrReal&  b)
MpfrReal& operator - (const unsigned int a, const MpfrReal&  b)
  {
  return -(b-a);
  }

// inline MpfrReal& operator - (const long int a, const MpfrReal&  b)
MpfrReal& operator - (const long int a, const MpfrReal&  b)
  {
  return -(b-a);
  }

// inline MpfrReal& operator - (const unsigned long int a, const MpfrReal&  b)
MpfrReal& operator - (const unsigned long int a, const MpfrReal&  b)
  {
  return -(b-a);
  }

MpfrReal& MpfrReal::operator - () const
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);                     
  mpfr_neg(res.mpfr_rep, mpfr_rep, CurrRndMode);  
  res.nbref = new RefCounter(1);
  return res;
  }

MpfrReal& MpfrReal::operator -= (const MpfrReal& b) 
  {
  if ( iszero(b) )
    return *this;
  if ( nbref->decr() <= 0 )			// the memory can be reused
    {
    if ( iszero(*this) )
      mpfr_neg(mpfr_rep, b.mpfr_rep, CurrRndMode);  
    else
      mpfr_sub(mpfr_rep, mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else						// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    if ( iszero(tmp) )
      mpfr_neg(mpfr_rep, b.mpfr_rep, CurrRndMode);  
    else
      mpfr_sub(mpfr_rep, tmp.mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator -= (const double b) 
MpfrReal& MpfrReal::operator -= (const double b) 
  {
  return (*this) -= MpfrReal(b);
  }

// inline MpfrReal& MpfrReal::operator -= (const int b) 
MpfrReal& MpfrReal::operator -= (const int b) 
  {
  (*this) += -b;
  return *this;
  }

// inline MpfrReal& MpfrReal::operator -= (const long int b) 
MpfrReal& MpfrReal::operator -= (const long int b) 
  {
  (*this) += -b;
  return *this;
  }

MpfrReal& MpfrReal::operator -= (const unsigned long int b) 	// using the more efficient mpfr_sub_ui
  {
  if ( b == 0 )
    return *this;
  if ( nbref->decr() <= 0 )			// the memory can be reused
    {
    mpfr_sub_ui(mpfr_rep, mpfr_rep, b, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else						// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    mpfr_sub_ui(mpfr_rep, tmp.mpfr_rep, b, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator -= (const unsigned int b) 
MpfrReal& MpfrReal::operator -= (const unsigned int b) 
  {
  (*this) -= (unsigned long int)b;
  return *this;
  }

MpfrReal& MpfrReal::sub (MpfrReal& res, const MpfrReal& r1, const MpfrReal& r2, RoundingModes rnd)
  {
  if ( iszero(r2) )
    return res = r1;
  else if ( iszero(r1) )
    return res = - r2;

  if ( res.nbref->decr() <= 0 )			// the memory can be reused
    {
    mpfr_sub(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref->refvalue() = 1;
    }
  else						// the previous value must be preserved
    {
    mpfr_init(res.mpfr_rep);
    mpfr_sub(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref = new RefCounter(1);
    }
  return res;
  }

// Multiplication-----------------------------------------------
MpfrReal& MpfrReal::operator * (const MpfrReal& b) const
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);                     
  mpfr_mul(res.mpfr_rep, mpfr_rep, b.mpfr_rep, CurrRndMode);  
  res.nbref = new RefCounter(1);
  return res;
  }

// inline MpfrReal& operator * (const MpfrReal& a, const double b) 
MpfrReal& operator * (const MpfrReal& a, const double b) 
  {
  return a * MpfrReal(b);
  }

// inline MpfrReal& operator * (const MpfrReal& a, const int b)
MpfrReal& operator * (const MpfrReal& a, const int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res *= b;
  return res;
  }

// inline MpfrReal& operator * (const MpfrReal& a, const unsigned int b)
MpfrReal& operator * (const MpfrReal& a, const unsigned int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res *= b;
  return res;
  }

// inline MpfrReal& operator * (const MpfrReal& a, const long int b)
MpfrReal& operator * (const MpfrReal& a, const long int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res *= b;
  return res;
  }

// inline MpfrReal& operator * (const MpfrReal& a, const unsigned long int b)
MpfrReal& operator * (const MpfrReal& a, const unsigned long int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res *= b;
  return res;
  }

// inline MpfrReal& operator * (const double a, const MpfrReal& b)
MpfrReal& operator * (const double a, const MpfrReal& b)
  {
  MpfrReal& res = * new MpfrReal;
  res = b;
  res *= MpfrReal(a);
  return res;
  }

// inline MpfrReal& operator * (const int a, const MpfrReal& b)
MpfrReal& operator * (const int a, const MpfrReal& b)
  {
  MpfrReal& res = * new MpfrReal;
  res = b;
  res *= a;
  return res;
  }

// inline MpfrReal& operator * (const unsigned int a, const MpfrReal& b)
MpfrReal& operator * (const unsigned int a, const MpfrReal& b)
  {
  MpfrReal& res = * new MpfrReal;
  res = b;
  res *= a;
  return res;
  }

// inline MpfrReal& operator * (const long int a, const MpfrReal& b)
MpfrReal& operator * (const long int a, const MpfrReal& b)
  {
  MpfrReal& res = * new MpfrReal;
  res = b;
  res *= a;
  return res;
  }

// inline MpfrReal& operator * (const unsigned long int a, const MpfrReal& b)
MpfrReal& operator * (const unsigned long int a, const MpfrReal& b)
  {
  MpfrReal& res = * new MpfrReal;
  res = b;
  res *= a;
  return res;
  }

MpfrReal& MpfrReal::operator *= (const MpfrReal& b) 
  {
  if ( nbref->decr() <= 0 )			// the memory can be reused
    {
    mpfr_mul(mpfr_rep, mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else						// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    mpfr_mul(mpfr_rep, tmp.mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator *= (const double b) 
MpfrReal& MpfrReal::operator *= (const double b) 
  {
  return (*this) *= MpfrReal(b);
  }

MpfrReal& MpfrReal::operator *= (const long int b) 	// using the more efficient mpfr_mul_ui
  {
  if ( nbref->decr() <= 0 )			// the memory can be reused
    {
    if (b >= 0)
      mpfr_mul_ui(mpfr_rep, mpfr_rep, (unsigned long int)(b), CurrRndMode);
    else
      {
      mpfr_mul_ui(mpfr_rep, mpfr_rep, (unsigned long int)(-b), CurrRndMode);
      mpfr_neg(mpfr_rep, mpfr_rep, CurrRndMode);
      }
    nbref->refvalue() = 1;
    }
  else						// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    if (b >= 0)
      mpfr_mul_ui(mpfr_rep, tmp.mpfr_rep, (unsigned long int)(b), CurrRndMode);
    else
      {
      mpfr_mul_ui(mpfr_rep, tmp.mpfr_rep, (unsigned long int)(-b), CurrRndMode);
      mpfr_neg(mpfr_rep, mpfr_rep, CurrRndMode);
      }
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator *= (const int b) 
MpfrReal& MpfrReal::operator *= (const int b) 
  {
  (*this)*= (long int) (b);
  return *this;
  }

MpfrReal& MpfrReal::operator *= (const unsigned long int b) 	// the efficient mpfr_mul_ui is used
  {
  if ( nbref->decr() <= 0 )			// the memory can be reused
    {
    mpfr_mul_ui(mpfr_rep, mpfr_rep, b, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else						// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    mpfr_mul_ui(mpfr_rep, tmp.mpfr_rep, b, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator *= (const unsigned int b) 
MpfrReal& MpfrReal::operator *= (const unsigned int b) 
  {
  (*this)*= (unsigned long int) (b);
  return *this;
  }

MpfrReal& MpfrReal::mul (MpfrReal& res, const MpfrReal& r1, const MpfrReal& r2, RoundingModes rnd)
  {
  if ( iszero(r1) )
    return res = r1;
  else if ( iszero(r2) )
    return res = r2;

  if ( res.nbref->decr() <= 0 )			// the memory can be reused
    {
    mpfr_mul(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref->refvalue() = 1;
    }
  else						// the previous value must be preserved
    {
    mpfr_init(res.mpfr_rep);
    mpfr_mul(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref = new RefCounter(1);
    }
  return res;
  }

// Division-----------------------------------------------------
MpfrReal& MpfrReal::operator / (const MpfrReal& b) const
  {

  //if ( iszero(b) )
    //{
    //cerr << "Division par 0 dans /= ... " << endl; // throw an exception?
    //NaN or infinity handled by MPFR
    //}

  MpfrReal& res = * new MpfrReal;
  mpfr_init(res.mpfr_rep);                     
  mpfr_div(res.mpfr_rep, mpfr_rep, b.mpfr_rep, CurrRndMode);  
  res.nbref = new RefCounter(1);
  return res;
  }

// inline MpfrReal& operator / (const MpfrReal& a, const double b) 
MpfrReal& operator / (const MpfrReal& a, const double b) 
  {
  return a / MpfrReal(b);
  }

// inline MpfrReal& operator / (const MpfrReal& a, const int b)
MpfrReal& operator / (const MpfrReal& a, const int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res /= b;
  return res;
  }
  
// inline MpfrReal& operator / (const MpfrReal& a, const unsigned int b)
MpfrReal& operator / (const MpfrReal& a, const unsigned int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res /= b;
  return res;
  }
  
// inline MpfrReal& operator / (const MpfrReal& a, const long int b)
MpfrReal& operator / (const MpfrReal& a, const long int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res /= b;
  return res;
  }
  
// inline MpfrReal& operator / (const MpfrReal& a, const unsigned long int b)
MpfrReal& operator / (const MpfrReal& a, const unsigned long int b)
  {
  MpfrReal& res = * new MpfrReal;
  res = a;
  res /= b;
  return res;
  }
  
// inline MpfrReal& operator / (const double a, const MpfrReal& b)
MpfrReal& operator / (const double a, const MpfrReal& b)
  {
  return MpfrReal(a)/b;
  }

// inline MpfrReal& operator / (const long int a, const MpfrReal& b)	// using mpfr_ui_div
MpfrReal& operator / (const long int a, const MpfrReal& b)	// using mpfr_ui_div
  {
  MpfrReal& res = * new MpfrReal;
  mpfr_init(res.mpfr_rep);
  //if ( iszero(b) )
    //{
    //cerr << "Division by 0..." << endl;
    //NaN or infinity handled by MPFR
    //}
  if (a >= 0)
    mpfr_ui_div(res.mpfr_rep, (unsigned long int) a, b.mpfr_rep, MpfrReal::CurrRndMode);
  else
    {
    mpfr_ui_div(res.mpfr_rep, (unsigned long int) (-a), b.mpfr_rep, MpfrReal::CurrRndMode);
    mpfr_neg(res.mpfr_rep, res.mpfr_rep, MpfrReal::CurrRndMode);
    }
  res.nbref = new RefCounter(1);
  return res;
  }

// inline MpfrReal& operator / (const int a, const MpfrReal& b)
MpfrReal& operator / (const int a, const MpfrReal& b)
  {
  return (long) a / b;
  }

// inline MpfrReal& operator / (const unsigned long int a, const MpfrReal& b)	// using mpfr_ui_div
MpfrReal& operator / (const unsigned long int a, const MpfrReal& b)	// using mpfr_ui_div
  {
  MpfrReal& res = * new MpfrReal;
  mpfr_init(res.mpfr_rep);
  //if ( iszero(b) )
    //{
    //cerr << "Division by 0..." << endl;
    //NaN or infinity handled by MPFR
    //}
  mpfr_ui_div(res.mpfr_rep, a, b.mpfr_rep, MpfrReal::CurrRndMode);
  res.nbref = new RefCounter(1);
  return res;
  }

// inline MpfrReal& operator / (const unsigned int a, const MpfrReal& b)
MpfrReal& operator / (const unsigned int a, const MpfrReal& b)
  {
  return (unsigned long int) a / b;
  }

MpfrReal& MpfrReal::operator /= (const MpfrReal& b) 
  {
  //if ( iszero(b) )
    //{
    //cerr << "Division par 0 dans /= ... " << endl; // throw an exception?
    //NaN or infinity handled by MPFR
    //}
  if ( nbref->decr() <= 0 )				// the memory can be reused
    {
    mpfr_div(mpfr_rep, mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else							// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    mpfr_div(mpfr_rep, tmp.mpfr_rep, b.mpfr_rep, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator /= (const double b) 
MpfrReal& MpfrReal::operator /= (const double b) 
  {
  return (*this) /= MpfrReal(b);
  }

MpfrReal& MpfrReal::operator /= (const long int b) 	// using the more efficient mpfr_div_ui
  {
  //if ( b == 0 )
    //{
    //cerr << "Division par 0 dans /= ... " << endl; // throw an exception?
    //NaN or infinity handled by MPFR
    //}
  if ( nbref->decr() <= 0 )				// the memory can be reused
    {
    if (b > 0)
      mpfr_div_ui(mpfr_rep, mpfr_rep, (unsigned long int)b, CurrRndMode);
    else
      {
      mpfr_div_ui(mpfr_rep, mpfr_rep, (unsigned long int)(-b), CurrRndMode);
      mpfr_neg(mpfr_rep, mpfr_rep, CurrRndMode);
      }
    nbref->refvalue() = 1;
    }
  else							// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    if (b > 0)
      mpfr_div_ui(mpfr_rep, tmp.mpfr_rep, (unsigned long int)b, CurrRndMode);
    else
      {
      mpfr_div_ui(mpfr_rep, tmp.mpfr_rep, (unsigned long int)(-b), CurrRndMode);
      mpfr_neg(mpfr_rep, mpfr_rep, CurrRndMode);
      }
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator /= (const int b)
MpfrReal& MpfrReal::operator /= (const int b)
  {
  (*this) /= (long int) b;
  return *this;
  }

MpfrReal& MpfrReal::operator /= (const unsigned long int b) 	// using mpfr_div_ui
  {
  //if ( b == 0 )
    //{
    //cerr << "Division par 0 dans /= ... " << endl; // throw an exception?
    //NaN or infinity handled by MPFR
    //}
  if ( nbref->decr() <= 0 )				// the memory can be reused
    {
    mpfr_div_ui(mpfr_rep, mpfr_rep, b, CurrRndMode);
    nbref->refvalue() = 1;
    }
  else							// the previous value must be preserved
    {
    MpfrReal tmp = *this;
    mpfr_init(mpfr_rep);
    mpfr_div_ui(mpfr_rep, tmp.mpfr_rep, b, CurrRndMode);
    nbref = new RefCounter(1);
    }
  return *this;
  }

// inline MpfrReal& MpfrReal::operator /= (const unsigned int b)
MpfrReal& MpfrReal::operator /= (const unsigned int b)
  {
  (*this) /= (unsigned long int) b;
  return *this;
  }

MpfrReal& MpfrReal::div (MpfrReal& res, const MpfrReal& r1, const MpfrReal& r2, RoundingModes rnd)
  {
  if ( iszero(r1) )
    return res = r1;
  //else if ( iszero(r2) )
    //{
    //cerr << "Division par 0 dans div..." << endl;
    //NaN or infinity handled by MPFR
    //}

  if ( res.nbref->decr() <= 0 )				// the memory can be reused
    {
    mpfr_div(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref->refvalue() = 1;
    }
  else							// the previous value must be preserved
    {
    mpfr_init(res.mpfr_rep);
    mpfr_div(res.mpfr_rep, r1.mpfr_rep, r2.mpfr_rep, CurrRndMode);
    res.nbref = new RefCounter(1);
    }
  return res;
  }


//--------------------------------------------------------------
//
// Mathematical and miscellaneous functions
//
//--------------------------------------------------------------
// NaN of infinity handled by MPFR 
// => no test on the value (sign...) of the operand
void MpfrReal::random (PrecisionType prec)		// member to avoid conflict with GMP random
  {
  if (nbref->decr() <= 0)				// the memory can be reused
    {
    mpfr_set_prec(mpfr_rep, prec);
    nbref->refvalue() = 1;
    }
  else							// the previous value must be preserved
    {
    mpfr_init2(mpfr_rep, prec);
    nbref = new RefCounter(1);
    }
  mpfr_random(mpfr_rep);
  }
  

MpfrReal& abs (const MpfrReal& r, RoundingModes rnd)
  {
  if ( r >= 0 )
    return r+0;
  else
    return -r;
  }
      
MpfrReal& sqrt (const MpfrReal& r, RoundingModes rnd)
  {
  MpfrReal& res = * new MpfrReal;

  // NaN handled by MPFR in case r < 0
  mpfr_init(res.mpfr_rep);
  mpfr_sqrt(res.mpfr_rep, r.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
  }
  
MpfrReal& exp (const MpfrReal& r, RoundingModes rnd)
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);
  mpfr_exp(res.mpfr_rep, r.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
  }
  
MpfrReal& log (const MpfrReal& r, RoundingModes rnd)
  {
  MpfrReal& res = * new MpfrReal;

  // NaN handled by MPFR in case r < 0
  mpfr_init(res.mpfr_rep);
  mpfr_log(res.mpfr_rep, r.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
  }
  
MpfrReal& sin (const MpfrReal& r, RoundingModes rnd)
  {
  MpfrReal& res = * new MpfrReal, dummy = * new MpfrReal;

  mpfr_init(res.mpfr_rep);
  mpfr_sin_cos(res.mpfr_rep, dummy.mpfr_rep, r.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
  }
  
MpfrReal& cos (const MpfrReal& r, RoundingModes rnd)
  {
  MpfrReal& res = * new MpfrReal, dummy = * new MpfrReal;

  mpfr_init(res.mpfr_rep);
  mpfr_sin_cos(dummy.mpfr_rep, res.mpfr_rep, r.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
  }
  
void sincos (MpfrReal& ressin, MpfrReal& rescos, const MpfrReal& r, RoundingModes rnd)
  {
  if ( ressin.nbref->decr() <= 0 )
    {
    ressin.nbref->refvalue() = 1;
    }
  else
    {
    mpfr_init(ressin.mpfr_rep);
    ressin.nbref = new RefCounter(1);
    }
  if ( rescos.nbref->decr() <= 0 )
    {
    rescos.nbref->refvalue() = 1;
    }
  else
    {
    mpfr_init(rescos.mpfr_rep);
    rescos.nbref = new RefCounter(1);
    }
  mpfr_sin_cos(ressin.mpfr_rep, rescos.mpfr_rep, r.mpfr_rep, rnd);
  }
  
/*
MpfrReal& zeta (const MpfrReal& r, RoundingModes rnd)
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);
  mpfr_zeta(res.mpfr_rep, r.mpfr_rep, rnd);
  res.nbref = new RefCounter(1);
  return res;
  }
*/
  
MpfrReal& floor (const MpfrReal& r)
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);
  mpfr_floor(res.mpfr_rep, r.mpfr_rep);
  res.nbref = new RefCounter(1);
  return res;
  }
  
MpfrReal& trunc (const MpfrReal& r)
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);
  mpfr_trunc(res.mpfr_rep, r.mpfr_rep);
  res.nbref = new RefCounter(1);
  return res;
  }
  
MpfrReal& ceil (const MpfrReal& r)
  {
  MpfrReal& res = * new MpfrReal;

  mpfr_init(res.mpfr_rep);
  mpfr_ceil(res.mpfr_rep, r.mpfr_rep);
  res.nbref = new RefCounter(1);
  return res;
  }

MpfrReal& add_one_ulp (const MpfrReal& r)
  {
  MpfrReal& res = * new MpfrReal;

  res.copy(r);
  mpfr_add_one_ulp(res.mpfr_rep);
  return res;
  }

MpfrReal& sub_one_ulp (const MpfrReal& r)
  {
  MpfrReal& res = * new MpfrReal;

  res.copy(r);
  mpfr_sub_one_ulp(res.mpfr_rep);
  return res;
  }