// C++ Interface for MPFR
// NR 21-11-00 (first steps in C++...)

#ifndef _MPFR_REAL_H_
#define _MPFR_REAL_H_

//#include <stdio.h>

// To be able to use the C++ string
using namespace std;
#ifdef __GNUC__
#include <string>
#define String string
#else
//#include <string_iso_SUNWCC.h>
#endif

#include <iostream.h>
#include <math.h>

extern "C" {
#include "gmp.h"
#include "gmp-impl.h"
#include "mpfr.h"
#include "mpfr-impl.h"
}

#ifndef _GIVARO_REF_COUNTER_H_
#define _GIVARO_REF_COUNTER_H_
// ==================================================================== //
// Definition of the Counter class, Counter
// (c) copyright GIVARO 1994                                            
// author: Th. Gautier, version : 2.7, date: 1995
// This class definition objects to handle reference
// counter for memory allocation (eg array0). 
// ==================================================================== //
#include <stddef.h>

class RefCounter {
public:
   // Cstor and Dstor
inline RefCounter( long l = 0) : counter(l) {} 
//inline RefCounter( const RefCounter& ) : counter(C.counter) {} 
inline ~RefCounter() {}

  //  Return the value
inline long  getvalue() const { return counter ; } 
inline long  val() const { return counter ; }
  // Return a ref to the counter
inline long& refvalue() { return counter ; }
  // Increments the counter and returns the new value 
inline long  incr() { return ++counter ; }
  // Decrements the value and returns the new value 
inline long  decr() { return --counter ; }

protected:
  long counter ;
} ;

#endif

//int my_mpfr_set_str(mpfr_ptr x, char *str, int base, mp_rnd_t rnd_mode);


// Rounding modes
typedef int RoundingModes;
static const RoundingModes RoundUp=int(GMP_RNDU); 
static const RoundingModes RoundDown=int(GMP_RNDD); 
static const RoundingModes RoundNearest=int(GMP_RNDN); 
static const RoundingModes RoundToZero=int(GMP_RNDZ);


class MpfrReal {

  protected:
    mpfr_t mpfr_rep;	// representation of the real in a mpfr format
    RefCounter *nbref;
    // Precision
    //typedef mp_prec_t PrecisionType;
    typedef long unsigned int PrecisionType;
    static PrecisionType &CurrPrecision;
    // Rounding mode
    static RoundingModes& CurrRndMode;

  public:
    // constructors and destructors
    MpfrReal ();
    MpfrReal (double d, 
              RoundingModes rnd = CurrRndMode, 
              PrecisionType prec = CurrPrecision);
    MpfrReal (int i, 
              RoundingModes rnd = CurrRndMode, 
              PrecisionType prec = CurrPrecision);
    MpfrReal (long int i, 
              RoundingModes rnd = CurrRndMode, 
              PrecisionType prec = CurrPrecision);
   MpfrReal (unsigned long int i, 
              RoundingModes rnd = CurrRndMode, 
              PrecisionType prec = CurrPrecision);
    MpfrReal(string s);
    MpfrReal (const MpfrReal& r);
    ~MpfrReal ();
 
    // Assignment and copy operators
    MpfrReal& operator = (const MpfrReal& r);
    MpfrReal& copy (const MpfrReal& r, 
                    RoundingModes rnd = CurrRndMode, 
                    PrecisionType prec = CurrPrecision);

    // Is equal to zero?
    friend int iszero(const MpfrReal& r) { return (!MPFR_IS_NAN(r.mpfr_rep) && !MPFR_NOTZERO(r.mpfr_rep)); }
    friend int isnan(const MpfrReal& r) { return (MPFR_IS_NAN(r.mpfr_rep)); }
    friend int sign(const MpfrReal& r) { return (MPFR_SIGN(r.mpfr_rep)); }

    // Precision and rounding mode
    static void SetDefaultPrecision (PrecisionType newprec);
    void SetPrecision (PrecisionType newprec);
    static const PrecisionType GetDefaultPrecision ();
    const PrecisionType GetPrecision ();
    static void SetDefaultRndMode (RoundingModes newrndmode);
    static const RoundingModes GetDefaultRndMode ();

    // Rounding: should be in place or not? In place
    // => to round not in place, copy and round
    void Round (RoundingModes rnd = CurrRndMode,
                PrecisionType prec = CurrPrecision);

    // "Constants" (but they depend on the precision and the rounding mode)
    static MpfrReal& Pi(RoundingModes rnd = CurrRndMode, 
                        PrecisionType prec = CurrPrecision) ;
    static MpfrReal& Log2(RoundingModes rnd = CurrRndMode, 
                        PrecisionType prec = CurrPrecision) ;

    // Comparison operators
    friend int compare (const MpfrReal& r1, const MpfrReal& r2);
    friend int compare (const MpfrReal& r1, const int r2);
    friend int compare (const MpfrReal& r1, const unsigned int r2);
    friend int compare (const MpfrReal& r1, const long int r2);
    friend int compare (const MpfrReal& r1, const unsigned long int r2);


    // Arithmetic operators
    // Philosophy: are members only the operations between MpfrReal ?
    //static MpfrReal& operator+ (const MpfrReal& r1, const MpfrReal& r2) ;
    MpfrReal& operator+ (const MpfrReal& r) const;
    friend MpfrReal& operator+ (const MpfrReal& r1, const double r2) ;
    friend MpfrReal& operator+ (const MpfrReal& r1, const int r2) ;
    friend MpfrReal& operator+ (const MpfrReal& r1, const unsigned int r2) ;
    friend MpfrReal& operator+ (const MpfrReal& r1, const long int r2) ;
    friend MpfrReal& operator+ (const MpfrReal& r1, const unsigned long int r2) ;
    friend MpfrReal& operator+ (const double r1, const MpfrReal& r2) ;
    friend MpfrReal& operator+ (const int r1, const MpfrReal& r2) ;
    friend MpfrReal& operator+ (const unsigned int r1, const MpfrReal& r2) ;
    friend MpfrReal& operator+ (const long int r1, const MpfrReal& r2) ;
    friend MpfrReal& operator+ (const unsigned long int r1, const MpfrReal& r2) ;
    MpfrReal& operator+= (const MpfrReal& r) ;
    MpfrReal& operator+= (const double r) ;
    MpfrReal& operator+= (const int r) ;
    MpfrReal& operator+= (const unsigned int r) ;
    MpfrReal& operator+= (const long int r) ;
    MpfrReal& operator+= (const unsigned long int r) ;
    static MpfrReal& add (MpfrReal& res, 
                          const MpfrReal& r1, const MpfrReal& r2,
                          RoundingModes rnd = CurrRndMode);

    //static MpfrReal& operator- (const MpfrReal& r1, const MpfrReal& r2);
    MpfrReal& operator- (const MpfrReal& r) const;
    friend MpfrReal& operator- (const MpfrReal& r1, const double r2) ;
    friend MpfrReal& operator- (const MpfrReal& r1, const int r2);
    friend MpfrReal& operator- (const MpfrReal& r1, const unsigned int r2);
    friend MpfrReal& operator- (const MpfrReal& r1, const long int r2);
    friend MpfrReal& operator- (const MpfrReal& r1, const unsigned long int r2);
    friend MpfrReal& operator- (const double r1, const MpfrReal& r2);
    friend MpfrReal& operator- (const int r1, const MpfrReal& r2);
    friend MpfrReal& operator- (const unsigned int r1, const MpfrReal& r2);
    friend MpfrReal& operator- (const long int r1, const MpfrReal& r2);
    friend MpfrReal& operator- (const unsigned long int r1, const MpfrReal& r2);
    MpfrReal& operator- () const;
    MpfrReal& operator-= (const MpfrReal& r) ;
    MpfrReal& operator-= (const double r) ;
    MpfrReal& operator-= (const int r) ;
    MpfrReal& operator-= (const long int r) ;
    MpfrReal& operator-= (const unsigned int r) ;
    MpfrReal& operator-= (const unsigned long int r) ;
    static MpfrReal& sub (MpfrReal& res, 
                          const MpfrReal& r1, const MpfrReal& r2,
                          RoundingModes rnd = CurrRndMode);

    //static MpfrReal& operator* (const MpfrReal& r1, const MpfrReal& r2) ;
    MpfrReal& operator* (const MpfrReal& r) const;
    friend MpfrReal& operator* (const MpfrReal& r1, const double r2) ;
    friend MpfrReal& operator* (const MpfrReal& r1, const int r2) ;
    friend MpfrReal& operator* (const MpfrReal& r1, const unsigned int r2) ;
    friend MpfrReal& operator* (const MpfrReal& r1, const long int r2) ;
    friend MpfrReal& operator* (const MpfrReal& r1, const unsigned long int r2) ;
    friend MpfrReal& operator* (const double r1, const MpfrReal& r2);
    friend MpfrReal& operator* (const int r1, const MpfrReal& r2);
    friend MpfrReal& operator* (const unsigned int r1, const MpfrReal& r2);
    friend MpfrReal& operator* (const long int r1, const MpfrReal& r2);
    friend MpfrReal& operator* (const unsigned long int r1, const MpfrReal& r2);
    MpfrReal& operator*= (const MpfrReal& r) ;
    MpfrReal& operator*= (const double r) ;
    MpfrReal& operator*= (const int r) ;
    MpfrReal& operator*= (const unsigned int r) ;
    MpfrReal& operator*= (const long int r) ;
    MpfrReal& operator*= (const unsigned long int r) ;
    static MpfrReal& mul (MpfrReal& res, 
                          const MpfrReal& r1, const MpfrReal& r2,
                          RoundingModes rnd = CurrRndMode);

    //static MpfrReal& operator/ (const MpfrReal& r1, const MpfrReal& r2) const;
    MpfrReal& operator/ (const MpfrReal& r) const;
    friend MpfrReal& operator/ (const MpfrReal& r1, const double r2) ;
    friend MpfrReal& operator/ (const MpfrReal& r1, const int r2) ;
    friend MpfrReal& operator/ (const MpfrReal& r1, const unsigned int r2) ;
    friend MpfrReal& operator/ (const MpfrReal& r1, const long int r2) ;
    friend MpfrReal& operator/ (const MpfrReal& r1, const unsigned long int r2) ;
    friend MpfrReal& operator/ (const double r1, const MpfrReal& r2) ;
    friend MpfrReal& operator/ (const int r1, const MpfrReal& r2) ;
    friend MpfrReal& operator/ (const unsigned int r1, const MpfrReal& r2) ;
    friend MpfrReal& operator/ (const long int r1, const MpfrReal& r2) ;
    friend MpfrReal& operator/ (const unsigned long int r1, const MpfrReal& r2) ;
    MpfrReal& operator/= (const MpfrReal& r) ;
    MpfrReal& operator/= (const double r) ;
    MpfrReal& operator/= (const int r) ;
    MpfrReal& operator/= (const unsigned int r) ;
    MpfrReal& operator/= (const long int r) ;
    MpfrReal& operator/= (const unsigned long int r) ;
    static MpfrReal& div (MpfrReal& res, 
                          const MpfrReal& r1, const MpfrReal& r2,
                          RoundingModes rnd = CurrRndMode);

    // Input/Output
    ostream& put(ostream& o,
                 RoundingModes rnd = CurrRndMode,
                 PrecisionType prec = CurrPrecision ) const;
                 //PrecisionType prec = PrecisionType(0) ) const;
    friend istream& operator >> (istream &i, MpfrReal& r);
    friend ostream& operator << (ostream &o, const MpfrReal& r);

    // Mathematical functions: exp, sin...
    void random (PrecisionType prec = CurrPrecision);

    friend MpfrReal& abs (const MpfrReal& r, RoundingModes rnd = CurrRndMode);
    friend MpfrReal& sqrt (const MpfrReal& r, RoundingModes rnd = CurrRndMode);
    friend MpfrReal& exp (const MpfrReal& r, RoundingModes rnd = CurrRndMode);
    friend MpfrReal& log (const MpfrReal& r, RoundingModes rnd = CurrRndMode);
    friend MpfrReal& sin (const MpfrReal& r, RoundingModes rnd = CurrRndMode);
    friend MpfrReal& cos (const MpfrReal& r, RoundingModes rnd = CurrRndMode);
    friend void sincos (MpfrReal& res_sin, MpfrReal& res_cos, 
                        const MpfrReal& r, RoundingModes rnd = CurrRndMode);
    //friend MpfrReal& zeta (const MpfrReal& r, RoundingModes rnd = CurrRndMode);

    friend MpfrReal& floor (const MpfrReal& r);
    friend MpfrReal& trunc (const MpfrReal& r);
    friend MpfrReal& ceil (const MpfrReal& r);

    friend MpfrReal& add_one_ulp (const MpfrReal& r);
    friend MpfrReal& sub_one_ulp (const MpfrReal& r);
};


//--------------------------
// Comparison operators
//--------------------------
int operator == (const MpfrReal& r1, const MpfrReal& r2);
int operator == (const MpfrReal& r1, const int r2);
int operator == (const MpfrReal& r1, const unsigned int r2);
int operator == (const MpfrReal& r1, const long int r2);
int operator == (const MpfrReal& r1, const unsigned long int r2);
int operator == (const int r1, const MpfrReal& r2);
int operator == (const unsigned int r1, const MpfrReal& r2);
int operator == (const long int r1, const MpfrReal& r2);
int operator == (const unsigned long int r1, const MpfrReal& r2);

int operator != (const MpfrReal& r1, const MpfrReal& r2);
int operator != (const MpfrReal& r1, const int r2);
int operator != (const MpfrReal& r1, const unsigned int r2);
int operator != (const MpfrReal& r1, const long int r2);
int operator != (const MpfrReal& r1, const unsigned long int r2);
int operator != (const int r1, const MpfrReal& r2);
int operator != (const unsigned int r1, const MpfrReal& r2);
int operator != (const long int r1, const MpfrReal& r2);
int operator != (const unsigned long int r1, const MpfrReal& r2);

int operator <  (const MpfrReal& r1, const MpfrReal& r2);
int operator <  (const MpfrReal& r1, const int r2);
int operator <  (const MpfrReal& r1, const unsigned int r2);
int operator <  (const MpfrReal& r1, const long int r2);
int operator <  (const MpfrReal& r1, const unsigned long int r2);
int operator <  (const int r1, const MpfrReal& r2);
int operator <  (const unsigned int r1, const MpfrReal& r2);
int operator <  (const long int r1, const MpfrReal& r2);
int operator <  (const unsigned long int r1, const MpfrReal& r2);

int operator <= (const MpfrReal& r1, const MpfrReal& r2);
int operator <= (const MpfrReal& r1, const int r2);
int operator <= (const MpfrReal& r1, const unsigned int r2);
int operator <= (const MpfrReal& r1, const long int r2);
int operator <= (const MpfrReal& r1, const unsigned long int r2);
int operator <= (const int r1, const MpfrReal& r2);
int operator <= (const unsigned int r1, const MpfrReal& r2);
int operator <= (const long int r1, const MpfrReal& r2);
int operator <= (const unsigned long int r1, const MpfrReal& r2);

int operator >  (const MpfrReal& r1, const MpfrReal& r2);
int operator >  (const MpfrReal& r1, const int r2);
int operator >  (const MpfrReal& r1, const unsigned int r2);
int operator >  (const MpfrReal& r1, const long int r2);
int operator >  (const MpfrReal& r1, const unsigned long int r2);
int operator >  (const int r1, const MpfrReal& r2);
int operator >  (const unsigned int r1, const MpfrReal& r2);
int operator >  (const long int r1, const MpfrReal& r2);
int operator >  (const unsigned long int r1, const MpfrReal& r2);

int operator >= (const MpfrReal& r1, const MpfrReal& r2);
int operator >= (const MpfrReal& r1, const int r2);
int operator >= (const MpfrReal& r1, const unsigned int r2);
int operator >= (const MpfrReal& r1, const long int r2);
int operator >= (const MpfrReal& r1, const unsigned long int r2);
int operator >= (const int r1, const MpfrReal& r2);
int operator >= (const unsigned int r1, const MpfrReal& r2);
int operator >= (const long int r1, const MpfrReal& r2);
int operator >= (const unsigned long int r1, const MpfrReal& r2);

const MpfrReal& min (const MpfrReal& a, const MpfrReal& b);
const MpfrReal& max (const MpfrReal& a, const MpfrReal& b);

#endif		// _MPFR_REAL_H_