Scl_mpf.hpp

Go to the documentation of this file.

/*********************************************************************
template<class T> inline 
Scl<T>& rem (const Scl<T>& dividend, const Scl<T>& divisor) NRVAL(r(dividend))
{
  NRCODE(r.rem(divisor);)
  NONRCODE(Scl<T> r(dividend); r.rem(divisor); return r;)
}
*      This file is part of the source code of BORDERBASIX software. *
*      (C) B. Mourrain, GALAAD, INRIA                                *
**********************************************************************
History: 
$Id: Scl_mpf.hpp,v 1.1.1.1 2006/10/06 08:01:40 trebuche Exp $
**********************************************************************/
#ifndef _SCL_MPF_H_
#define _SCL_MPF_H_
//----------------------------------------------------------------------
#include <iostream>
#include <iomanip>
#include <string>
#include <borderbasix/mdebug.hpp>
#include "Scl.hpp"
//----------------------------------------------------------------------
typedef mpf_t   MPF;
#define _MPF_
unsigned long MPF_PRECISION;
//----------------------------------------------------------------------
//template<> inline shared_object<MPF>::rep::~rep() {mpf_clear(&obj);}
//----------------------------------------------------------------------
//----------------------------------------------------------------------
template<> inline Scl<MPF>::~Scl (){mpf_clear(rep());}
//----------------------------------------------------------------------
template<> inline
void Scl<MPF>::init ( )                   {mpf_init(rep());}
//------------------- Scl<T> (...) -------------------------------------
template<> inline
Scl<MPF>::Scl ( )                    {
    //  cout<<"default construct "<<endl;
    mpf_init(rep());}
//----------------------------------------------------------------------
template<> inline 
Scl<MPF>::operator int ()  const
{
    return (int) mpf_get_d(rep());
}

template<> inline
double Scl<MPF>::castDouble(Scl<MPF> &ref)
{
    return mpf_get_d(ref.rep());
}

template<> inline 
Scl<MPF>::Scl(const Scl<MPF>& rhs) 
{
    //  COUNT<MPF>('c');
    //cout<<"const SCl "<<endl;
    mpf_init_set(data,rhs.rep());
}
//----------------------------------------------------------------------
template<> inline 
Scl<MPF> Scl<MPF>::operator =  (const Scl<MPF>& rhs)
{
    //COUNT<MPF>('=');
    if (this != &rhs) {
        //  mpf_init(rep());
        mpf_set(rep(), rhs.rep());
    }
    return *this;
}
//----------------------------------------------------------------------
template<> inline
Scl<MPF>::Scl (signed long int sl)   {mpf_init_set_si(rep(),sl);}
//----------------------------------------------------------------------
template<> inline
Scl<MPF>::Scl (unsigned long int ul) {mpf_init_set_ui(rep(),ul);}
//----------------------------------------------------------------------
template<> inline
Scl<MPF>::Scl (int si)               {
    // cout<<"init set si "<<endl;
    mpf_init_set_si(rep(), si);}
//----------------------------------------------------------------------
template<> inline
Scl<MPF>::Scl (const char *string, unsigned int base)
{
    if (mpf_init_set_str(rep(), string, base))
        std::cerr << "Scl<MPF>: The string " << string
                  << " is not a valid number in base " << base << std::endl;
}
template<> inline
Scl<MPF>::Scl (double d)
{
    //cout<<"je passe dans construct double "<<endl;
    mpf_init_set_d(rep(),d);
}
// == ------------------------------------------------------------------
template<> inline 
bool Scl<MPF>::operator == (const Scl<MPF>& rhs) const
{
    return mpf_cmp(rep(), rhs.rep()) == 0;
}

template<> inline 
bool Scl<MPF>::operator == (long sl) const {
    return mpf_cmp_si(rep(), sl) == 0;
}

template<> inline 
bool Scl<MPF>::operator == (int si) const {
    return mpf_cmp_si(rep(), (long) si) == 0;
}

template<> inline 
bool Scl<MPF>::operator == (unsigned long ul) const
{
    return mpf_cmp_ui(rep(), ul) == 0;
}
template<> inline 
bool Scl<MPF>::operator == (double ul) const
{
    return mpf_cmp_d(rep(), ul) == 0;
}

// != ------------------------------------------------------------------
template<> inline 
bool Scl<MPF>::operator != (const Scl<MPF>& rhs) const
{
    return mpf_cmp(rep(), rhs.rep()) != 0;
}
template<> inline
bool Scl<MPF>::operator != (double sl) const {
    return mpf_cmp_d(rep(), sl) != 0;
}


template<> inline 
bool Scl<MPF>::operator != (long sl) const {
    return mpf_cmp_si(rep(), sl) != 0;
}

template<> inline
bool Scl<MPF>::operator != (int si) const {
    return mpf_cmp_si(rep(), (long) si) != 0;
}

template<> inline 
bool Scl<MPF>::operator != (unsigned long ul) const
{
    return mpf_cmp_ui(rep(), ul) != 0;
}
// > -------------------------------------------------------------------
template<> inline 
bool Scl<MPF>::operator > (const Scl<MPF>& rhs) const
{
    //)<<"RR > RR";
    return mpf_cmp(rep(), rhs.rep()) > 0;
}

template<> inline 
bool Scl<MPF>::operator > (double sl) const {
    return mpf_cmp_d(rep(), sl) > 0;
}

template<> inline 
bool Scl<MPF>::operator > (long sl) const {
    return mpf_cmp_si(rep(), sl) > 0;
}

template<> inline 
bool Scl<MPF>::operator > (int si) const {
    return mpf_cmp_si(rep(), (long) si) > 0;
}

template<> inline 
bool Scl<MPF>::operator > (unsigned long ul) const
{
    return mpf_cmp_ui(rep(), ul) > 0;
}


// <= --------------------------------------------------------------------
template<> inline 
bool Scl<MPF>::operator >= (const Scl<MPF>& rhs) const
{
    return mpf_cmp(rep(), rhs.rep()) >= 0;
}

template<> inline 
bool Scl<MPF>::operator >= (long sl) const {
    return mpf_cmp_si(rep(), sl) >= 0;
}
template<> inline 
bool Scl<MPF>::operator >= (double sl) const {
    return mpf_cmp_d(rep(), sl) >= 0;
}

template<> inline 
bool Scl<MPF>::operator >= (int si) const
{
    return mpf_cmp_si(rep(), (long) si) >= 0;
}

template<> inline 
bool Scl<MPF>::operator >= (unsigned long ul) const
{
    return mpf_cmp_ui(rep(), ul) >= 0;
}

// < -------------------------------------------------------------------
template<> inline 
bool Scl<MPF>::operator < (const Scl<MPF>& rhs) const
{
    return mpf_cmp(rep(), rhs.rep()) < 0;
}

template<> inline 
bool Scl<MPF>::operator < (long sl) const {
    return mpf_cmp_si(rep(), sl) < 0;
}
template<> inline 
bool Scl<MPF>::operator < (double sl) const {
    return mpf_cmp_d(rep(), sl) < 0;
}

template<> inline 
bool Scl<MPF>::operator < (int si) const
{
    return mpf_cmp_si(rep(), (long) si) < 0;
}


template<> inline 
bool Scl<MPF>::operator < (unsigned long ul) const
{
    return mpf_cmp_ui(rep(), ul) < 0;
}

// <= ------------------------------------------------------------------

template<> inline 
bool Scl<MPF>::operator <= (const Scl<MPF>& rhs) const
{
    return mpf_cmp(rep(), rhs.rep()) <= 0;
}
template<> inline
bool Scl<MPF>::operator <= (double sl) const {
    return mpf_cmp_d(rep(), sl) <= 0;
}


template<> inline 
bool Scl<MPF>::operator <= (long sl) const {
    return mpf_cmp_si(rep(), sl) <= 0;
}

template<> inline 
bool Scl<MPF>::operator <= (int si) const
{
    return mpf_cmp_si(rep(), (long) si) <= 0;
}

template<> inline 
bool Scl<MPF>::operator <= (unsigned long ul) const
{
    return mpf_cmp_ui(rep(), ul) <= 0;
}

// = -------------------------------------------------------------------
template<> inline 
Scl<MPF>  Scl<MPF>::operator =  (unsigned long ul)
{
    mpf_set_ui(rep(), ul); return *this;
}

template<> inline
Scl<MPF>  Scl<MPF>::operator =  (double ul)
{
    mpf_set_d(rep(), ul); return *this;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator =  (long sl)
{
    mpf_set_si(rep(), sl);  return *this;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator = (int ul)
{
    mpf_set_si(rep(), ul);  return *this;
}
// += -------------------------------------------------------------------
template<> inline 
Scl<MPF> Scl<MPF>::operator += (const Scl<MPF>& rhs)
{
    Scl<MPF> tmp;
    //mpf_init(tmp.rep());
    mpf_add(tmp.rep(), rep(), rhs.rep());
    mpf_swap(tmp.rep(),rep());
    return *this;

}

template<> inline 
Scl<MPF> Scl<MPF>::operator += (unsigned long ul)
{
    mpf_add_ui(rep(), rep(), ul); return *this;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator += (long sl)
{
    if (sl >= 0)
        mpf_add_ui(rep(), rep(), (unsigned long) sl);
    else
        mpf_sub_ui(rep(), rep(), ((unsigned long) -sl));
    return *this;
}

template<> inline
Scl<MPF> Scl<MPF>::operator += (int sl)
{
    if (sl >= 0)
        mpf_add_ui(rep(), rep(), (unsigned long) sl);
    else
        mpf_sub_ui(rep(), rep(), ((unsigned long) -sl));
    return *this;
}


template<> inline 
Scl<MPF> Scl<MPF>::operator += (double ul)
{
    *this += Scl<MPF>(ul);  return *this;
}

// -= -------------------------------------------------------------------
template<> inline 
Scl<MPF> Scl<MPF>::operator -= (const Scl<MPF>& rhs)
{
    mpf_sub(rep(), rep(), rhs.rep()); return *this;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator -= (double ul)
{
    Scl<MPF> tmp(ul);
    *this-=tmp; return *this;
}
template<> inline 
Scl<MPF> Scl<MPF>::operator -= (unsigned long ul)
{
    mpf_sub_ui(rep(), rep(), ul); return *this;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator -= (long sl)
{
    if (sl >= 0)
        mpf_sub_ui(rep(), rep(), (unsigned long) sl);
    else
        mpf_add_ui(rep(), rep(), (unsigned long) -sl);
    return *this;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator -= (int sl)
{
    if (sl >= 0)
        mpf_sub_ui(rep(), rep(), (unsigned long) sl);
    else
        mpf_add_ui(rep(), rep(), (unsigned long) -sl);
    return *this;
    //  *this -= Scl<MPF>(ul); return *this;
}
// *= -----------------------------------------------------------------
template<> inline 
Scl<MPF> Scl<MPF>::operator *= (const Scl<MPF>& rhs)
{
    mpf_mul(rep(), rep(), rhs.rep()); return *this;
}
template<> inline 
Scl<MPF> Scl<MPF>::operator *= (unsigned long ul)
{
    mpf_mul_ui(rep(), rep(), ul); return *this;
}
template<> inline 
Scl<MPF> Scl<MPF>::operator *= (long sl)
{
    if (sl >= 0)
        mpf_mul_ui(rep(), rep(), (unsigned long) sl);
    else
    {
        mpf_neg(rep(),rep());
        mpf_mul_ui(rep(), rep(),(unsigned long)(-sl));
    }
    return *this;
}
template<> inline 
Scl<MPF> Scl<MPF>::operator *= (double d)
{
    return *this*=(Scl<MPF>)d;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator *= (int ul)
{
    if (ul >= 0)
        mpf_mul_ui(rep(), rep(), (unsigned long) ul);
    else {
        mpf_mul_ui(rep(), rep(), (unsigned long) (-ul));
        mpf_neg(rep(), rep());
    }
    return *this;
}
// /= ------------------------------------------------------------------
template<> 
Scl<MPF> Scl<MPF>::operator /= (const Scl<MPF>& rhs)
{
    mpf_div(rep(), rep(), rhs.rep());
    //mpf_divexact(&rep(), &rep(), &rhs.rep());
    return *this;
}
template<> inline 
Scl<MPF> Scl<MPF>::operator /= (unsigned long ul)
{
    mpf_div_ui(rep(), rep(), ul); return *this;
}

template<> inline 
Scl<MPF> Scl<MPF>::operator /= (long sl)
{
    if (sl >= 0)
        mpf_div_ui(rep(), rep(), (unsigned long) sl);
    else
    {
        mpf_neg(rep(),rep());
        mpf_div_ui(rep(), rep(), ((unsigned long) -sl));
    }
    return *this;
}
template<> inline 
Scl<MPF> Scl<MPF>::operator /= (int sl)
{
    if (sl >= 0)
        mpf_div_ui(rep(), rep(), (unsigned long) sl);
    else
    {
        mpf_neg(rep(),rep());
        mpf_div_ui(rep(), rep(), ((unsigned long) -sl));
    }
    return *this;
}

//======================================================================
// ARITHMETIC OPERATORS
//======================================================================
#if 0
Scl<MPF> operator -(const Scl<MPF>& a1) 
{
    Scl<MPF> result(a1);
    result*=(-1);
    return result;
}
#endif 
Scl<MPF> operator +(const Scl<MPF>& a1, const Scl<MPF>& a2) 
{
    Scl<MPF> result;
    mpf_add(result.rep(), a1.rep(), a2.rep());
    return result;
}
// //----------------------------------------------------------------------
Scl<MPF> operator -(const Scl<MPF>& a1, const Scl<MPF>& a2) 
{
    Scl<MPF> result;
    mpf_sub(result.rep(), a1.rep(), a2.rep());
    return result;
}
//----------------------------------------------------------------------
Scl<MPF> operator *(const Scl<MPF>& a1, const Scl<MPF>& a2) 
{
    Scl<MPF> result;
    mpf_mul(result.rep(), a1.rep(), a2.rep());
    return result;
}
Scl<MPF> operator *(const Scl<MPF>& a1, double a2)
{
    Scl<MPF> result(a1);
    result*=a2;
    return result;
}

//----------------------------------------------------------------------
Scl<MPF> operator /(const Scl<MPF>& a1, const Scl<MPF>& a2) 
{
    Scl<MPF> result;
    mpf_div(result.rep(), a1.rep(), a2.rep());
    return result;
}
Scl<MPF> operator /(double a2,const Scl<MPF>& a1)
{
    Scl<MPF> result(a2);
    result/=a1;
    return result;
}
Scl<MPF> operator /(const Scl<MPF>& a1, double a2)
{
    Scl<MPF> result=a1;
    Scl<MPF>tmp(a2);
    result/=tmp;
    return result;
}
Scl<MPF> operator /(int a2,const Scl<MPF>& a1)
{
    Scl<MPF> result(a2);
    result/=a1;
    return result;
}
Scl<MPF> operator /(const Scl<MPF>& a1, int a2)
{
    Scl<MPF> result(a1);
    Scl<MPF>tmp(a2);
    result/=tmp;
    return result;
}

Scl<MPF> operator /(long int a2,const Scl<MPF>& a1)
{
    Scl<MPF> result(a2);
    result/=a1;
    return result;
}
//Scl<MPF> operator /(const Scl<MPF>& a1, long int a2)
//{
//  Scl<MPF> result(a1);
//  Scl<MPF>tmp(a2);
//  result/=tmp;
//  return result;
//}


//======================================================================
// ASSIGNEMENTS  
//======================================================================
//  inline void assign(Scl<MPF> & p, 
//          const OP<'+',Scl<MPF>, Scl<MPF> >& M) 
//  {
//    mpf_add(p.rep(),M.op1.rep(),M.op2.rep());
//  }
//  //----------------------------------------------------------------------
//  inline
//  void assign(Scl<MPF> & p, 
//          const OP<'-',Scl<MPF>, Scl<MPF> >& M) 
//  { 
//    mpf_sub(p.rep(),M.op1.rep(),M.op2.rep());
//  }
//  //----------------------------------------------------------------------
//  inline void assign(Scl<MPF> & p, 
//          const OP<'*',Scl<MPF>, Scl<MPF> >& M) 
//  { 
//    mpf_mul(p.rep(),M.op1.rep(),M.op2.rep());
//  }
//  //----------------------------------------------------------------------
//  void assign(Scl<MPF> & p, 
//          const OP<'/',Scl<MPF>,Scl<MPF> >& M) 
//  { 
//    assert(M.op2!=0);
//    mpf_div(p.rep(),M.op1.rep(),M.op2.rep());
//}
//======================================================================
template<> inline 
void Scl<MPF>::operator ++ ( )
{
    mpf_add_ui(rep(), rep(), 1);
}
//----------------------------------------------------------------------
template<> inline 
void Scl<MPF>::operator -- ( )
{
    mpf_sub_ui(rep(), rep(), 1);
}
//======================================================================
inline void convert(Scl<MPF>& n, char *s) 
{
    // printf("%s\n",s);
    mpf_init_set_str(n.rep(), s, 10);
    //cout<<"apres convert "<<mpf_get_d(n.rep())<<endl;;
}
//======================================================================
// INPUT OUTPUT
//======================================================================
std::ostream& operator << (std::ostream& os, const Scl<MPF>& b)
{
    // os<<setprecision(40);
    unsigned sz= mpf_get_prec(b.rep())+1;
    char* tmp = new char [sz];
    int n=os.precision();
    gmp_sprintf(tmp,"%.*Ff",n,b.rep());
    os<< tmp;
    return os;

    mp_exp_t exponent ;
    char* str = mpf_get_str(NULL,&exponent, 10, 0, b.rep());

    if(mpf_cmp_si(b.rep(),0)<0) {
        os<< "-0."<<(++str);
    } else {
        os<< "0."<<str;
    }
    if (exponent != 0)
        os <<"e"<<exponent;
    //mpf_out_str(stdout, 10,16,b.rep());
    return os;
}
//----------------------------------------------------------------------
std::istream& operator >> (std::istream& is, Scl<MPF>& b)
{
    std::string s;is >> s;
    //cout<<"ici string "<<s<<endl;
    mpf_init_set_str(b.rep(),s.c_str(), 10);
    //cout<<"valeur de la converstion "<<b<<endl;
    //mp_exp_t exponent ;
    //is<< mpf_get_str(NULL,&exponent, 10,0, b.rep());
    return is;
}

bool has_positive_sign(const Scl<MPF>& b) {
    return (b>0);
}

//======================================================================
// SPECIAL FUNCTIONS
//======================================================================
void mpf_precision(unsigned long l)
{
    mpf_set_default_prec(l);
}
//----------------------------------------------------------------------
void Precision( Scl<MPF>& b, unsigned long l)
{
    mpf_set_prec(b.rep(),l);
}
//======================================================================
template<> inline 
void Scl<MPF>::Div2Exp (unsigned long exponent_of_2)
{
    mpf_div_2exp(rep(), rep(), exponent_of_2);
}

template<> inline void
Scl<MPF>::abs ( )
{
    if (mpf_sgn(rep()) < 0)
        mpf_neg(rep(),rep());
}

template<> inline 
Scl<MPF>& Scl<MPF>::negate ( )
{
    mpf_neg(rep(),rep());
    //rep()._mp_size = - rep()._mp_size;
    return *this;
}

template<> inline 
void Scl<MPF>::sqrt ( )
{
    mpf_sqrt(rep(), rep());
}

//-------------------------- log(...) -----------------------------------
#include<mpfr.h>
inline Scl<MPF> log (const Scl<MPF>& b)
{
    mpfr_t tmp;
    Scl<MPF> res;
    mpfr_set_default_prec(mpf_get_default_prec()+10);
    mpfr_init(tmp);
    mpfr_set_f(tmp,b.rep(),GMP_RNDN);
    mpfr_log(tmp,tmp,GMP_RNDN);
    mpfr_get_f(res.rep(),tmp,GMP_RNDN);
    mpfr_clear(tmp);
    return res;

    //return mpf_size(b.rep());
}
//-------------------------- sign ---------------------------------------------
inline int sign (const Scl<MPF>& b)
{
    return mpf_sgn(b.rep());
    //  if (b.rep()._mp_size == 0)
    //  return 0;
    //else
    //  return b.rep()._mp_size > 0 ? 1 : -1;
}
//--------------------------------------------------------------------------------
//#ifdef ONEMORETIME


template<> inline
Scl<MPF>& Scl<MPF>:: pow (int exp){

    mpf_pow_ui(rep(),rep(),exp);
    return *this;
}


inline Scl<MPF> pow(const Scl<MPF>& b1, const Scl<MPF>& b2)
{
    //cout<<"puissance "<<b1<<" "<<b2<<endl;
    mpfr_t tmpb1,tmpb2;
    Scl<MPF> res;
    if(1)
    {
        // mpfr_set_default_prec(mpf_get_default_prec()+10);
        mpfr_init(tmpb1);
        mpfr_init(tmpb2);
        mpfr_set_f(tmpb1,b1.rep(),MPFR_RNDN);
        mpfr_set_f(tmpb2,b2.rep(),MPFR_RNDN);
        mpfr_pow(tmpb1,tmpb1,tmpb2,MPFR_RNDN);
        mpfr_get_f(res.rep(),tmpb1,MPFR_RNDN);
        //cout<<"res puissance pos"<<res<<endl;
        mpfr_clear(tmpb1);
        mpfr_clear(tmpb2);
    }
    else
    {
        mpfr_init(tmpb1);
        mpfr_init(tmpb2);
        mpfr_set_f(tmpb1,b1.rep(),MPFR_RNDN);
        mpfr_set_f(tmpb2,b2.rep(),MPFR_RNDN);
        mpfr_neg(tmpb2,tmpb2,MPFR_RNDN);
        //mpfr_printf("le tmpb1 : %Rf\n",tmpb1);
        //mpfr_printf("le tmpb2 : %Rf\n",tmpb2);
        mpfr_pow(tmpb1,tmpb1,tmpb2,MPFR_RNDN);
        //mpfr_printf("le temp : %Rf\n",tmpb1);
        fflush(stdout);
        mpfr_ui_div(tmpb1,1l,tmpb1,MPFR_RNDN);
        mpfr_get_f(res.rep(),tmpb1,MPFR_RNDN);
        //cout<<"res puissance neg "<<res<<endl;
        mpfr_clear(tmpb1);
        mpfr_clear(tmpb2);

    }
    return res;

}

//-------------------------- compare (...) ------------------------------------
inline int compare (const Scl<MPF>& b1, const Scl<MPF>& b2)
{
    return mpf_cmp(b1.rep(), b2.rep());
}

inline int compare (const Scl<MPF>& b, unsigned long ul)
{
    return mpf_cmp_ui(b.rep(), ul);
}

inline int compare (const Scl<MPF>& b, long sl)
{
    return mpf_cmp_si(b.rep(), sl);
}

bool operator > (double ul, const Scl<MPF> &rhs)
{
    //mdebug()<<"double > RR"<<ul<<rhs;
    return mpf_cmp_d(rhs.rep(), ul) < 0;
}
bool operator < (double ul, const Scl<MPF> &rhs)
{
    return mpf_cmp_d(rhs.rep(), ul) > 0;
}
bool operator < (int ul, const Scl<MPF> &rhs)
{
    return mpf_cmp_si(rhs.rep(), ul) > 0;
}
bool operator > (int ul, const Scl<MPF> &rhs)
{
    //mdebug()<<"int > RR"<<ul<<rhs<<(mpf_cmp_si(rhs.rep(), ul) < 0);
    return mpf_cmp_si(rhs.rep(), ul) < 0;
}
//----------------------------------------------------------------------
#endif // //SCL_MPF_H