Scl_mpfr.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;)
}
*      (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_MPFR_H_
#define _SCL_MPFR_H_
//----------------------------------------------------------------------
//#include <cmath>
#include <iostream>
#include <iomanip>
#include <string>
#include <borderbasix/mdebug.hpp>
#include "Scl.hpp"
#include <mpfr.h>

//----------------------------------------------------------------------
typedef mpfr_t   MPFR;
#define _MPFR_
unsigned long MPFR_PRECISION;
//----------------------------------------------------------------------

//----------------------------------------------------------------------
template<> inline Scl<MPFR>::~Scl (){
    mpfr_clear(rep());
}
//----------------------------------------------------------------------
template<> inline
void Scl<MPFR>::init ( ) {
    mpfr_init(rep());
}
//----------------------------------------------------------------------
template<> inline
Scl<MPFR>::Scl ( )                    {
    //  cout<<"default construct "<<endl;
    mpfr_init(rep());
}
//----------------------------------------------------------------------
template<> inline 
Scl<MPFR>::operator int ()  const
{
    return (int) mpfr_get_d(rep(),MPFR_RNDN);
}
//----------------------------------------------------------------------
template<> inline
double Scl<MPFR>::castDouble(Scl<MPFR> &ref)
{
    return mpfr_get_d(ref.rep(),MPFR_RNDN);
}
//----------------------------------------------------------------------
template<> inline 
Scl<MPFR>::Scl(const Scl<MPFR>& rhs)
{
    //cout<<"const SCl "<<endl;
    mpfr_init_set(rep(),rhs.rep(),MPFR_RNDN);
}

//----------------------------------------------------------------------
template<> inline
Scl<MPFR>::Scl (signed long int sl)   {
    mpfr_init_set_si(rep(),sl,MPFR_RNDN);
}
//----------------------------------------------------------------------
template<> inline
Scl<MPFR>::Scl (unsigned long int ul) {
    mpfr_init_set_ui(rep(),ul,MPFR_RNDN);
}
//----------------------------------------------------------------------
template<> inline
Scl<MPFR>::Scl (int si) {
    // cout<<"init set si "<<endl;
    mpfr_init_set_si(rep(), si,MPFR_RNDN);
}
//----------------------------------------------------------------------
template<> inline
Scl<MPFR>::Scl (const char *str, unsigned int base)
{
    mpfr_set_default_prec(64);
    if (mpfr_init_set_str(rep(), str, base, MPFR_RNDN))
        std::cerr << "Scl<MPFR>: The string " << str
                  << " is not a valid number in base " << base << std::endl;
}
//----------------------------------------------------------------------
template<> inline
Scl<MPFR>::Scl (double d)
{
    //cout<<"je passe dans construct double "<<endl;
    mpfr_init_set_d(rep(),d,MPFR_RNDN);
}
//----------------------------------------------------------------------
template<> inline 
Scl<MPFR> Scl<MPFR>::operator =  (const Scl<MPFR>& rhs)
{
    //COUNT<MPFR>('=');
    if (this != &rhs) {
        //  mpfr_init(rep());
        mpfr_set(rep(), rhs.rep(),MPFR_RNDN);
    }
    return *this;
}

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

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

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

template<> inline 
bool Scl<MPFR>::operator == (unsigned long ul) const
{
    return mpfr_cmp_ui(rep(), ul) == 0;
}
template<> inline 
bool Scl<MPFR>::operator == (double ul) const
{
    return mpfr_cmp_d(rep(), ul) == 0;
}

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


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

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

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

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

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

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

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


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

template<> inline 
bool Scl<MPFR>::operator >= (long sl) const {
    return mpfr_cmp_si(rep(), sl) >= 0;
}
template<> inline 
bool Scl<MPFR>::operator >= (double sl) const {
    return mpfr_cmp_d(rep(), sl) >= 0;
}

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

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

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

template<> inline 
bool Scl<MPFR>::operator < (long sl) const {
    return mpfr_cmp_si(rep(), sl) < 0;
}
template<> inline 
bool Scl<MPFR>::operator < (double sl) const {
    return mpfr_cmp_d(rep(), sl) < 0;
}

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


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

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

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


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

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

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

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

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

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

template<> inline 
Scl<MPFR> Scl<MPFR>::operator = (int ul)
{
    mpfr_set_si(rep(), ul,MPFR_RNDN);  return *this;
}
// += -------------------------------------------------------------------
template<> inline 
Scl<MPFR> Scl<MPFR>::operator += (const Scl<MPFR>& rhs)
{
    Scl<MPFR> tmp;
    //mpfr_init(tmp.rep());
    mpfr_add(tmp.rep(), rep(), rhs.rep(),MPFR_RNDN);
    mpfr_swap(tmp.rep(),rep());
    return *this;

}

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

template<> inline 
Scl<MPFR> Scl<MPFR>::operator += (long sl)
{
    if (sl >= 0)
        mpfr_add_ui(rep(), rep(), (unsigned long) sl,MPFR_RNDN);
    else
        mpfr_sub_ui(rep(), rep(), ((unsigned long) -sl),MPFR_RNDN);
    return *this;
}

template<> inline
Scl<MPFR> Scl<MPFR>::operator += (int sl)
{
    if (sl >= 0)
        mpfr_add_ui(rep(), rep(), (unsigned long) sl,MPFR_RNDN);
    else
        mpfr_sub_ui(rep(), rep(), ((unsigned long) -sl),MPFR_RNDN);
    return *this;
}


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

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

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

template<> inline 
Scl<MPFR> Scl<MPFR>::operator -= (long sl)
{
    if (sl >= 0)
        mpfr_sub_ui(rep(), rep(), (unsigned long) sl,MPFR_RNDN);
    else
        mpfr_add_ui(rep(), rep(), (unsigned long) -sl,MPFR_RNDN);
    return *this;
}

template<> inline 
Scl<MPFR> Scl<MPFR>::operator -= (int sl)
{
    if (sl >= 0)
        mpfr_sub_ui(rep(), rep(), (unsigned long) sl,MPFR_RNDN);
    else
        mpfr_add_ui(rep(), rep(), (unsigned long) -sl,MPFR_RNDN);
    return *this;
    //  *this -= Scl<MPFR>(ul); return *this;
}
// *= -----------------------------------------------------------------
template<> inline 
Scl<MPFR> Scl<MPFR>::operator *= (const Scl<MPFR>& rhs)
{
    mpfr_mul(rep(), rep(), rhs.rep(),MPFR_RNDN);
    return *this;
}
template<> inline 
Scl<MPFR> Scl<MPFR>::operator *= (unsigned long ul)
{
    mpfr_mul_ui(rep(), rep(), ul,MPFR_RNDN);
    return *this;
}
template<> inline 
Scl<MPFR> Scl<MPFR>::operator *= (long sl)
{
    if (sl >= 0)
        mpfr_mul_ui(rep(), rep(), (unsigned long) sl,MPFR_RNDN);
    else
    {
        mpfr_neg(rep(),rep(),MPFR_RNDN);
        mpfr_mul_ui(rep(), rep(),(unsigned long)(-sl),MPFR_RNDN);
    }
    return *this;
}
template<> inline 
Scl<MPFR> Scl<MPFR>::operator *= (double d)
{
    return *this*=(Scl<MPFR>)d;
}

template<> inline 
Scl<MPFR> Scl<MPFR>::operator *= (int ul)
{
    if (ul >= 0)
        mpfr_mul_ui(rep(), rep(), (unsigned long) ul,MPFR_RNDN);
    else {
        mpfr_mul_ui(rep(), rep(), (unsigned long) (-ul),MPFR_RNDN);
        mpfr_neg(rep(), rep(),MPFR_RNDN);
    }
    return *this;
}
// /= ------------------------------------------------------------------
template<> 
Scl<MPFR> Scl<MPFR>::operator /= (const Scl<MPFR>& rhs)
{
    mpfr_div(rep(), rep(), rhs.rep(),MPFR_RNDN);
    //mpfr_divexact(&rep(), &rep(), &rhs.rep());
    return *this;
}
template<> inline 
Scl<MPFR> Scl<MPFR>::operator /= (unsigned long ul)
{
    mpfr_div_ui(rep(), rep(), ul,MPFR_RNDN);
    return *this;
}

template<> inline 
Scl<MPFR> Scl<MPFR>::operator /= (long sl)
{
    if (sl >= 0)
        mpfr_div_ui(rep(), rep(), (unsigned long) sl,MPFR_RNDN);
    else
    {
        mpfr_neg(rep(),rep(),MPFR_RNDN);
        mpfr_div_ui(rep(), rep(), ((unsigned long) -sl),MPFR_RNDN);
    }
    return *this;
}
template<> inline 
Scl<MPFR> Scl<MPFR>::operator /= (int sl)
{
    if (sl >= 0)
        mpfr_div_ui(rep(), rep(), (unsigned long) sl,MPFR_RNDN);
    else
    {
        mpfr_neg(rep(),rep(),MPFR_RNDN);
        mpfr_div_ui(rep(), rep(), ((unsigned long) -sl),MPFR_RNDN);
    }
    return *this;
}

//======================================================================
// ARITHMETIC OPERATORS
//======================================================================
#if 0
Scl<MPFR> operator -(const Scl<MPFR>& a1)
{
    Scl<MPFR> result(a1);
    result*=(-1);
    return result;
}
#endif 
Scl<MPFR> operator +(const Scl<MPFR>& a1, const Scl<MPFR>& a2)
{
    Scl<MPFR> result;
    mpfr_add(result.rep(), a1.rep(), a2.rep(),MPFR_RNDN);
    return result;
}
// //----------------------------------------------------------------------
Scl<MPFR> operator -(const Scl<MPFR>& a1, const Scl<MPFR>& a2)
{
    Scl<MPFR> result;
    mpfr_sub(result.rep(), a1.rep(), a2.rep(),MPFR_RNDN);
    return result;
}
//----------------------------------------------------------------------
Scl<MPFR> operator *(const Scl<MPFR>& a1, const Scl<MPFR>& a2)
{
    Scl<MPFR> result;
    mpfr_mul(result.rep(), a1.rep(), a2.rep(),MPFR_RNDN);
    return result;
}
Scl<MPFR> operator *(const Scl<MPFR>& a1, double a2)
{
    Scl<MPFR> result(a1);
    result*=a2;
    return result;
}

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

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



template<> inline 
void Scl<MPFR>::operator ++ ( )
{
    mpfr_add_ui(rep(), rep(), 1,MPFR_RNDN);
}
//----------------------------------------------------------------------
template<> inline 
void Scl<MPFR>::operator -- ( )
{
    mpfr_sub_ui(rep(), rep(), 1,MPFR_RNDN);
}
//======================================================================
inline void convert(Scl<MPFR>& n, char *s)
{
    // printf("%s\n",s);
    mpfr_init_set_str(n.rep(), s, 10,MPFR_RNDN);
    //cout<<"apres convert "<<mpfr_get_d(n.rep())<<endl;;
}

//======================================================================
// INPUT OUTPUT
//======================================================================
std::ostream& operator << (std::ostream& os, const Scl<MPFR>& x)
{
    //std::cout<<"\nrrr ";
    //mpfr_out_str(stdout,10,mpfr_get_prec(x.rep()),x.rep(),MPFR_RNDN);
    //std::cout<<"\n";

    //os<<std::setprecision(40);
    unsigned sz= mpfr_get_prec(x.rep())+1;
    //std::cout<<" prec "<<sz<<std::endl;
    char* tmp = new char [sz];

    int n = os.precision();
    mpfr_exp_t e;
    mpfr_get_str (tmp+1, &e, 10, n, x.rep(), MPFR_RNDN);

    int s=(x<0?2:1);
    tmp[0] = tmp[1];
    tmp[s-1]= tmp[s];
    tmp[s] = '.';


    //mpfr_sprintf(tmp,"%.RNf",x.rep());
    os<< tmp;
    if(e!=1)
        os << "e" << e-1;
    //delete[] tmp;
    return os;
}
//----------------------------------------------------------------------
std::istream& operator >> (std::istream& is, Scl<MPFR>& b)
{
    std::string s; is >> s;
    //cout<<"ici string "<<s<<endl;
    mpfr_init_set_str(b.rep(), s.c_str(), 10, MPFR_RNDN);

    return is;
}

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

//======================================================================
// SPECIAL FUNCTIONS
//======================================================================
void Precision(unsigned long l)
{
    mpfr_set_default_prec(l);
}
//----------------------------------------------------------------------
void Precision( Scl<MPFR>& b, unsigned long l)
{
    mpfr_set_prec(b.rep(),l);
}
//======================================================================
template<> inline 
void Scl<MPFR>::Div2Exp (unsigned long exponent_of_2)
{
    mpfr_div_2exp(rep(), rep(), exponent_of_2,MPFR_RNDN);
}

template<> inline void
Scl<MPFR>::abs ( )
{
    if (mpfr_sgn(rep()) < 0)
        mpfr_neg(rep(),rep(),MPFR_RNDN);
}

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

template<> inline 
void Scl<MPFR>::sqrt ( )
{
    mpfr_sqrt(rep(), rep(),MPFR_RNDN);
}

//-------------------------- log(...) -----------------------------------
inline Scl<MPFR> log (const Scl<MPFR>& b)
{
    mpfr_t tmp;
    Scl<MPFR> res;
    mpfr_set_default_prec(mpfr_get_default_prec()+10);
    mpfr_init(tmp);
    mpfr_set(tmp,b.rep(),MPFR_RNDN);
    mpfr_log(tmp,tmp,MPFR_RNDN);
    mpfr_set(res.rep(),tmp,MPFR_RNDN);
    mpfr_clear(tmp);
    return res;

    //return mpfr_size(b.rep());
}
//-------------------------- sign ---------------------------------------------
inline int sign (const Scl<MPFR>& b)
{
    return mpfr_sgn(b.rep());
}
//--------------------------------------------------------------------------------
//#ifdef ONEMORETIME


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

    mpfr_pow_ui(rep(),rep(),exp,MPFR_RNDN);
    return *this;
}


inline Scl<MPFR> pow(const Scl<MPFR>& b1, const Scl<MPFR>& b2)
{
    //cout<<"puissance "<<b1<<" "<<b2<<endl;
    mpfr_t tmpb1,tmpb2;
    Scl<MPFR> res;
    if(1)
    {
        // mpfr_set_default_prec(mpfr_get_default_prec()+10);
        mpfr_init(tmpb1);
        mpfr_init(tmpb2);
        mpfr_set(tmpb1,b1.rep(),MPFR_RNDN);
        mpfr_set(tmpb2,b2.rep(),MPFR_RNDN);
        mpfr_pow(tmpb1,tmpb1,tmpb2,MPFR_RNDN);
        mpfr_set(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(tmpb1,b1.rep(),MPFR_RNDN);
        mpfr_set(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_set(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<MPFR>& b1, const Scl<MPFR>& b2)
{
    return mpfr_cmp(b1.rep(), b2.rep());
}

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

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

bool operator > (double ul, const Scl<MPFR> &rhs)
{
    //mdebug()<<"double > RR"<<ul<<rhs;
    return mpfr_cmp_d(rhs.rep(), ul) < 0;
}
bool operator < (double ul, const Scl<MPFR> &rhs)
{
    return mpfr_cmp_d(rhs.rep(), ul) > 0;
}
bool operator < (int ul, const Scl<MPFR> &rhs)
{
    return mpfr_cmp_si(rhs.rep(), ul) > 0;
}
bool operator > (int ul, const Scl<MPFR> &rhs)
{
    //mdebug()<<"int > RR"<<ul<<rhs<<(mpfr_cmp_si(rhs.rep(), ul) < 0);
    return mpfr_cmp_si(rhs.rep(), ul) < 0;
}
//----------------------------------------------------------------------
Scl<MPFR> fmax(const Scl<MPFR>& a, const Scl<MPFR>& b) {
    if(b>a)
        return b;
    else
        return a;
}

Scl<MPFR> logb(const Scl<MPFR>& x) {
    mpfr_t tmp;
    Scl<MPFR> res;
    mpfr_set_default_prec(mpfr_get_default_prec()+10);
    mpfr_init(tmp);
    mpfr_set(tmp,x.rep(),GMP_RNDN);
    if(x<0)
        mpfr_neg(tmp,tmp,GMP_RNDN);
    mpfr_log2(tmp,tmp,GMP_RNDN);
    mpfr_set(res.rep(),tmp,GMP_RNDN);
    mpfr_clear(tmp);
    return res;
}

int isinf(const Scl<MPFR>& x) {
    return mpfr_inf_p(x.rep());
}

int isfinite(const Scl<MPFR>& x) {
    return (mpfr_inf_p(x.rep())==0?1:0);
}

int isnan(const Scl<MPFR>& x) {
    return mpfr_nan_p(x.rep());
}

Scl<MPFR> fabs(const Scl<MPFR>& x) {
    Scl<MPFR> result(x);
    result.abs();
    return result;
}

Scl<MPFR> copysign(const Scl<MPFR>& a, const Scl<MPFR>& b) {
    Scl<MPFR> res;
    mpfr_copysign(res.rep(),a.rep(),b.rep(), MPFR_RNDN);
    return res;
}
namespace std {
    Scl<MPFR> abs(const Scl<MPFR>& x) {
        return fabs(x);
    }
}
#endif // //SCL_MPF_H