apps/nick/mpreal.cc

/*
  This file is part of MADNESS.

  Copyright (C) 2007,2010 Oak Ridge National Laboratory

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

  For more information please contact:

  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367

  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680

  $Id$
*/
/*
	Multi-precision real number class. C++ wrapper fo MPFR library.
	Project homepage: http://www.holoborodko.com/pavel/
	Contact e-mail:   pavel@holoborodko.com

	Copyright (c) 2008 Pavel Holoborodko

	This library is free software; you can redistribute it and/or
	modify it under the terms of the GNU Lesser General Public
	License as published by the Free Software Foundation; either
	version 2.1 of the License, or (at your option) any later version.

	This library is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	Lesser General Public License for more details.

	You should have received a copy of the GNU Lesser General Public
	License along with this library; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

	Contributors:
	Brain Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze.
*/

#include <cstring>
#include <sstream>
#include "mpreal.h"

using std::ws;
using std::cerr;
using std::endl;
using std::string;
using std::ostream;
using std::istream;

namespace mpfr{

mp_rnd_t   mpreal::default_rnd  = mpfr_get_default_rounding_mode();
mp_prec_t  mpreal::default_prec = mpfr_get_default_prec();
int		   mpreal::default_base = 10;
int        mpreal::double_bits = -1;

// Default constructor: creates mp number and initializes it to 0.
mpreal::mpreal()
{
	mpfr_init2(mp,default_prec);
	mpfr_set_ui(mp,0,default_rnd);
}

mpreal::mpreal(const mpreal& u)
{
	mpfr_init2(mp,mpfr_get_prec(u.mp));
	mpfr_set(mp,u.mp,default_rnd);
}

mpreal::mpreal(const mpfr_t u)
{
	mpfr_init2(mp,mpfr_get_prec(u));
	mpfr_set(mp,u,default_rnd);
}

mpreal::mpreal(const mpf_t u)
{
	mpfr_init2(mp,mpf_get_prec(u));
	mpfr_set_f(mp,u,default_rnd);
}

mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
{
	mpfr_init2(mp,prec);
	mpfr_set_z(mp,u,default_rnd);
}

mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
{
	mpfr_init2(mp,prec);
	mpfr_set_q(mp,u,default_rnd);
}

mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
{
    if(double_bits == -1 || fits_in_bits(u, double_bits))
    {
    	mpfr_init2(mp,prec);
	    mpfr_set_d(mp,u,mode);
    }
    else
        throw conversion_overflow();
}

mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
{
    mpfr_init2(mp,prec);
	mpfr_set_ld(mp,u,mode);
}

mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
{
	mpfr_init2(mp,prec);
	mpfr_set_ui(mp,u,mode);
}

mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
{
	mpfr_init2(mp,prec);
	mpfr_set_ui(mp,u,mode);
}

mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
{
	mpfr_init2(mp,prec);
	mpfr_set_si(mp,u,mode);
}

mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
{
	mpfr_init2(mp,prec);
	mpfr_set_si(mp,u,mode);
}

mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
{
	mpfr_init2(mp,prec);
	mpfr_set_str(mp, s, base, mode);
}

mpreal::~mpreal()
{
	mpfr_clear(mp);
}

// Operators - Assignment
mpreal& mpreal::operator=(const char* s)
{
	mpfr_t t;
	if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
	{
		mpfr_set(mp,t,mpreal::default_rnd);
		mpfr_clear(t);
	}else{
		mpfr_clear(t);
		// cerr<<"fail to convert string"<<endl;
	}

	return *this;
}

const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
{
	mpreal a;
	mp_prec_t p1, p2, p3;

	p1 = v1.get_prec();
	p2 = v2.get_prec();
	p3 = v3.get_prec();

	a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));

	mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
	return a;
}

const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
{
	mpreal a;
	mp_prec_t p1, p2, p3;

	p1 = v1.get_prec();
	p2 = v2.get_prec();
	p3 = v3.get_prec();

	a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));

	mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
	return a;
}

const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
{
	mpreal a;
	mp_prec_t p1, p2;

	p1 = v1.get_prec();
	p2 = v2.get_prec();

	a.set_prec(p1>p2?p1:p2);

	mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);

	return a;
}

const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{
	mpreal a;
	mp_prec_t yp, xp;

	yp = y.get_prec();
	xp = x.get_prec();

	a.set_prec(yp>xp?yp:xp);

	mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);

	return a;
}

const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
{
	mpreal x;
	mpfr_ptr* t;
	unsigned long int i;

	t = new mpfr_ptr[n];
	for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp;
	mpfr_sum(x.mp,t,n,rnd_mode);
	delete[] t;
	return x;
}

const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{
	mpreal a;
	mp_prec_t yp, xp;

	yp = y.get_prec();
	xp = x.get_prec();

	a.set_prec(yp>xp?yp:xp);

	mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);

	return a;
}

const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{
	mpreal a;
	mp_prec_t yp, xp;

	yp = y.get_prec();
	xp = x.get_prec();

	a.set_prec(yp>xp?yp:xp);

	mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);

	return a;
}

template <class T>
std::string to_string(T t, std::ios_base & (*f)(std::ios_base&))
{
	std::ostringstream oss;
	oss << f << t;
	return oss.str();
}

string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
{
	char *s, *ns = nullptr;
	size_t slen, nslen;
	mp_exp_t exp;
	string out;

	if(mpfr_inf_p(mp))
	{
		if(mpfr_sgn(mp)>0) return "+@Inf@";
		else			   return "-@Inf@";
	}

	if(mpfr_zero_p(mp)) return "0";
	if(mpfr_nan_p(mp))  return "@NaN@";


	s  = mpfr_get_str(nullptr,&exp,b,0,mp,mode);
	ns = mpfr_get_str(nullptr,&exp,b,n,mp,mode);

	if(s!=nullptr && ns!=nullptr)
	{
		slen  = strlen(s);
		nslen = strlen(ns);
		if(nslen<=slen)
		{
			mpfr_free_str(s);
			s = ns;
			slen = nslen;
		}
		else {
			mpfr_free_str(ns);
		}

		// Make human eye-friendly formatting if possible
		if (exp>0 && static_cast<size_t>(exp)<slen)
		{
			if(s[0]=='-')
			{
				// Remove zeros starting from right end
				char* ptr = s+slen-1;
				while (*ptr=='0' && ptr>s+exp) ptr--;

				if(ptr==s+exp) out = string(s,exp+1);
				else		   out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1);

				//out = string(s,exp+1)+'.'+string(s+exp+1);
			}
			else
			{
				// Remove zeros starting from right end
				char* ptr = s+slen-1;
				while (*ptr=='0' && ptr>s+exp-1) ptr--;

				if(ptr==s+exp-1) out = string(s,exp);
				else		     out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1);

				//out = string(s,exp)+'.'+string(s+exp);
			}

		}else{ // exp<0 || exp>slen
			if(s[0]=='-')
			{
				// Remove zeros starting from right end
				char* ptr = s+slen-1;
				while (*ptr=='0' && ptr>s+1) ptr--;

				if(ptr==s+1) out = string(s,2);
				else		 out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1);

				//out = string(s,2)+'.'+string(s+2);
			}
			else
			{
				// Remove zeros starting from right end
				char* ptr = s+slen-1;
				while (*ptr=='0' && ptr>s) ptr--;

				if(ptr==s) out = string(s,1);
				else	   out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1);

				//out = string(s,1)+'.'+string(s+1);
			}

			// Make final string
			if(--exp)
				if(exp>0) out += "e+"+mpfr::to_string<mp_exp_t>(exp,std::dec);
				else 	  out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec);
		}

		mpfr_free_str(s);
		return out;
	}else{
		return "conversion error!";
	}
}

//////////////////////////////////////////////////////////////////////////
// I/O
ostream& operator<<(ostream& os, const mpreal& v)
{
	return os<<v.to_string(os.precision());
}

istream& operator>>(istream &is, mpreal& v)
{
	char c;
	string s = "";
	mpfr_t t;

	if(is.good())
	{
		is>>ws;
		while ((c = is.get())!=EOF)
		{
			if(c ==' ' || c == '\t' || c == '\n' || c == '\r')
			{
				is.putback(c);
				break;
			}
			s += c;
		}

		if(s.size() != 0)
		{
			// Protect current value from alternation in case of input error
			// so some error handling(roll back) procedure can be used
			if(0==mpfr_init_set_str(t,s.c_str(),mpreal::default_base,mpreal::default_rnd))
			{
				mpfr_set(v.mp,t,mpreal::default_rnd);
				mpfr_clear(t);

			}else{
				mpfr_clear(t);
				cerr<<"error reading from istream"<<endl;
				// throw an exception
			}
		}
	}
	return is;
}
}