This commit is contained in:
Benoit Jacob 2010-11-26 10:10:26 -05:00
parent e868b6736a
commit 139392488d

View File

@ -1,408 +1,408 @@
/* /*
Multi-precision real number class. C++ wrapper fo MPFR library. Multi-precision real number class. C++ wrapper fo MPFR library.
Project homepage: http://www.holoborodko.com/pavel/ Project homepage: http://www.holoborodko.com/pavel/
Contact e-mail: pavel@holoborodko.com Contact e-mail: pavel@holoborodko.com
Copyright (c) 2008-2010 Pavel Holoborodko Copyright (c) 2008-2010 Pavel Holoborodko
This library is free software; you can redistribute it and/or This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version. 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, This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details. Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Contributors: Contributors:
Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze,
Heinz van Saanen, Pere Constans, Dmitriy Gubanov Heinz van Saanen, Pere Constans, Dmitriy Gubanov
*/ */
#include <cstring> #include <cstring>
#include "mpreal.h" #include "mpreal.h"
using std::ws; using std::ws;
using std::cerr; using std::cerr;
using std::endl; using std::endl;
using std::string; using std::string;
using std::ostream; using std::ostream;
using std::istream; using std::istream;
namespace mpfr{ namespace mpfr{
mp_rnd_t mpreal::default_rnd = mpfr_get_default_rounding_mode(); mp_rnd_t mpreal::default_rnd = mpfr_get_default_rounding_mode();
mp_prec_t mpreal::default_prec = mpfr_get_default_prec(); mp_prec_t mpreal::default_prec = mpfr_get_default_prec();
int mpreal::default_base = 10; int mpreal::default_base = 10;
int mpreal::double_bits = -1; int mpreal::double_bits = -1;
// Default constructor: creates mp number and initializes it to 0. // Default constructor: creates mp number and initializes it to 0.
mpreal::mpreal() mpreal::mpreal()
{ {
mpfr_init2(mp,default_prec); mpfr_init2(mp,default_prec);
mpfr_set_ui(mp,0,default_rnd); mpfr_set_ui(mp,0,default_rnd);
} }
mpreal::mpreal(const mpreal& u) mpreal::mpreal(const mpreal& u)
{ {
mpfr_init2(mp,mpfr_get_prec(u.mp)); mpfr_init2(mp,mpfr_get_prec(u.mp));
mpfr_set(mp,u.mp,default_rnd); mpfr_set(mp,u.mp,default_rnd);
} }
mpreal::mpreal(const mpfr_t u) mpreal::mpreal(const mpfr_t u)
{ {
mpfr_init2(mp,mpfr_get_prec(u)); mpfr_init2(mp,mpfr_get_prec(u));
mpfr_set(mp,u,default_rnd); mpfr_set(mp,u,default_rnd);
} }
mpreal::mpreal(const mpf_t u) mpreal::mpreal(const mpf_t u)
{ {
mpfr_init2(mp,mpf_get_prec(u)); mpfr_init2(mp,mpf_get_prec(u));
mpfr_set_f(mp,u,default_rnd); mpfr_set_f(mp,u,default_rnd);
} }
mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_z(mp,u,mode); mpfr_set_z(mp,u,mode);
} }
mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_q(mp,u,mode); mpfr_set_q(mp,u,mode);
} }
mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
{ {
if(double_bits == -1 || fits_in_bits(u, double_bits)) if(double_bits == -1 || fits_in_bits(u, double_bits))
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_d(mp,u,mode); mpfr_set_d(mp,u,mode);
} }
else else
throw conversion_overflow(); throw conversion_overflow();
} }
mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_ld(mp,u,mode); mpfr_set_ld(mp,u,mode);
} }
mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_ui(mp,u,mode); mpfr_set_ui(mp,u,mode);
} }
mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_ui(mp,u,mode); mpfr_set_ui(mp,u,mode);
} }
mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_si(mp,u,mode); mpfr_set_si(mp,u,mode);
} }
mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_si(mp,u,mode); mpfr_set_si(mp,u,mode);
} }
mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
{ {
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_str(mp, s, base, mode); mpfr_set_str(mp, s, base, mode);
} }
mpreal::~mpreal() mpreal::~mpreal()
{ {
mpfr_clear(mp); mpfr_clear(mp);
} }
// Operators - Assignment // Operators - Assignment
mpreal& mpreal::operator=(const char* s) mpreal& mpreal::operator=(const char* s)
{ {
mpfr_t t; mpfr_t t;
if(0==mpfr_init_set_str(t,s,default_base,default_rnd)) if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
{ {
mpfr_set(mp,t,mpreal::default_rnd); mpfr_set(mp,t,mpreal::default_rnd);
mpfr_clear(t); mpfr_clear(t);
}else{ }else{
mpfr_clear(t); mpfr_clear(t);
// cerr<<"fail to convert string"<<endl; // cerr<<"fail to convert string"<<endl;
} }
return *this; return *this;
} }
const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a;
mp_prec_t p1, p2, p3; mp_prec_t p1, p2, p3;
p1 = v1.get_prec(); p1 = v1.get_prec();
p2 = v2.get_prec(); p2 = v2.get_prec();
p3 = v3.get_prec(); p3 = v3.get_prec();
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); 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); mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
return a; return a;
} }
const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a;
mp_prec_t p1, p2, p3; mp_prec_t p1, p2, p3;
p1 = v1.get_prec(); p1 = v1.get_prec();
p2 = v2.get_prec(); p2 = v2.get_prec();
p3 = v3.get_prec(); p3 = v3.get_prec();
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); 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); mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
return a; return a;
} }
const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode) const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a;
mp_prec_t p1, p2; mp_prec_t p1, p2;
p1 = v1.get_prec(); p1 = v1.get_prec();
p2 = v2.get_prec(); p2 = v2.get_prec();
a.set_prec(p1>p2?p1:p2); a.set_prec(p1>p2?p1:p2);
mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode); mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);
return a; return a;
} }
const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a;
mp_prec_t yp, xp; mp_prec_t yp, xp;
yp = y.get_prec(); yp = y.get_prec();
xp = x.get_prec(); xp = x.get_prec();
a.set_prec(yp>xp?yp:xp); a.set_prec(yp>xp?yp:xp);
mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode); mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);
return a; return a;
} }
const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode) const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
{ {
mpreal x; mpreal x;
mpfr_ptr* t; mpfr_ptr* t;
unsigned long int i; unsigned long int i;
t = new mpfr_ptr[n]; t = new mpfr_ptr[n];
for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp; for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp;
mpfr_sum(x.mp,t,n,rnd_mode); mpfr_sum(x.mp,t,n,rnd_mode);
delete[] t; delete[] t;
return x; return x;
} }
const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a;
mp_prec_t yp, xp; mp_prec_t yp, xp;
yp = y.get_prec(); yp = y.get_prec();
xp = x.get_prec(); xp = x.get_prec();
a.set_prec(yp>xp?yp:xp); a.set_prec(yp>xp?yp:xp);
mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode); mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);
return a; return a;
} }
const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a;
mp_prec_t yp, xp; mp_prec_t yp, xp;
yp = y.get_prec(); yp = y.get_prec();
xp = x.get_prec(); xp = x.get_prec();
a.set_prec(yp>xp?yp:xp); a.set_prec(yp>xp?yp:xp);
mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode); mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);
return a; return a;
} }
template <class T> template <class T>
std::string to_string(T t, std::ios_base & (*f)(std::ios_base&)) std::string to_string(T t, std::ios_base & (*f)(std::ios_base&))
{ {
std::ostringstream oss; std::ostringstream oss;
oss << f << t; oss << f << t;
return oss.str(); return oss.str();
} }
mpreal::operator std::string() const mpreal::operator std::string() const
{ {
return to_string(); return to_string();
} }
string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
{ {
char *s, *ns = NULL; char *s, *ns = NULL;
size_t slen, nslen; size_t slen, nslen;
mp_exp_t exp; mp_exp_t exp;
string out; string out;
if(mpfr_inf_p(mp)) if(mpfr_inf_p(mp))
{ {
if(mpfr_sgn(mp)>0) return "+@Inf@"; if(mpfr_sgn(mp)>0) return "+@Inf@";
else return "-@Inf@"; else return "-@Inf@";
} }
if(mpfr_zero_p(mp)) return "0"; if(mpfr_zero_p(mp)) return "0";
if(mpfr_nan_p(mp)) return "@NaN@"; if(mpfr_nan_p(mp)) return "@NaN@";
s = mpfr_get_str(NULL,&exp,b,0,mp,mode); s = mpfr_get_str(NULL,&exp,b,0,mp,mode);
ns = mpfr_get_str(NULL,&exp,b,n,mp,mode); ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
if(s!=NULL && ns!=NULL) if(s!=NULL && ns!=NULL)
{ {
slen = strlen(s); slen = strlen(s);
nslen = strlen(ns); nslen = strlen(ns);
if(nslen<=slen) if(nslen<=slen)
{ {
mpfr_free_str(s); mpfr_free_str(s);
s = ns; s = ns;
slen = nslen; slen = nslen;
} }
else { else {
mpfr_free_str(ns); mpfr_free_str(ns);
} }
// Make human eye-friendly formatting if possible // Make human eye-friendly formatting if possible
if (exp>0 && static_cast<size_t>(exp)<slen) if (exp>0 && static_cast<size_t>(exp)<slen)
{ {
if(s[0]=='-') if(s[0]=='-')
{ {
// Remove zeros starting from right end // Remove zeros starting from right end
char* ptr = s+slen-1; char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s+exp) ptr--; while (*ptr=='0' && ptr>s+exp) ptr--;
if(ptr==s+exp) out = string(s,exp+1); if(ptr==s+exp) out = string(s,exp+1);
else out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+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); //out = string(s,exp+1)+'.'+string(s+exp+1);
} }
else else
{ {
// Remove zeros starting from right end // Remove zeros starting from right end
char* ptr = s+slen-1; char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s+exp-1) ptr--; while (*ptr=='0' && ptr>s+exp-1) ptr--;
if(ptr==s+exp-1) out = string(s,exp); if(ptr==s+exp-1) out = string(s,exp);
else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1); else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1);
//out = string(s,exp)+'.'+string(s+exp); //out = string(s,exp)+'.'+string(s+exp);
} }
}else{ // exp<0 || exp>slen }else{ // exp<0 || exp>slen
if(s[0]=='-') if(s[0]=='-')
{ {
// Remove zeros starting from right end // Remove zeros starting from right end
char* ptr = s+slen-1; char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s+1) ptr--; while (*ptr=='0' && ptr>s+1) ptr--;
if(ptr==s+1) out = string(s,2); if(ptr==s+1) out = string(s,2);
else out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1); else out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1);
//out = string(s,2)+'.'+string(s+2); //out = string(s,2)+'.'+string(s+2);
} }
else else
{ {
// Remove zeros starting from right end // Remove zeros starting from right end
char* ptr = s+slen-1; char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s) ptr--; while (*ptr=='0' && ptr>s) ptr--;
if(ptr==s) out = string(s,1); if(ptr==s) out = string(s,1);
else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1); else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1);
//out = string(s,1)+'.'+string(s+1); //out = string(s,1)+'.'+string(s+1);
} }
// Make final string // Make final string
if(--exp) if(--exp)
{ {
if(exp>0) out += "e+"+mpfr::to_string<mp_exp_t>(exp,std::dec); 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); else out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec);
} }
} }
mpfr_free_str(s); mpfr_free_str(s);
return out; return out;
}else{ }else{
return "conversion error!"; return "conversion error!";
} }
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
// I/O // I/O
ostream& operator<<(ostream& os, const mpreal& v) ostream& operator<<(ostream& os, const mpreal& v)
{ {
return os<<v.to_string(os.precision()); return os<<v.to_string(os.precision());
} }
istream& operator>>(istream &is, mpreal& v) istream& operator>>(istream &is, mpreal& v)
{ {
char c; char c;
string s = ""; string s = "";
mpfr_t t; mpfr_t t;
if(is.good()) if(is.good())
{ {
is>>ws; is>>ws;
while ((c = is.get())!=EOF) while ((c = is.get())!=EOF)
{ {
if(c ==' ' || c == '\t' || c == '\n' || c == '\r') if(c ==' ' || c == '\t' || c == '\n' || c == '\r')
{ {
is.putback(c); is.putback(c);
break; break;
} }
s += c; s += c;
} }
if(s.size() != 0) if(s.size() != 0)
{ {
// Protect current value from alternation in case of input error // Protect current value from alternation in case of input error
// so some error handling(roll back) procedure can be used // 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)) 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_set(v.mp,t,mpreal::default_rnd);
mpfr_clear(t); mpfr_clear(t);
}else{ }else{
mpfr_clear(t); mpfr_clear(t);
cerr<<"error reading from istream"<<endl; cerr<<"error reading from istream"<<endl;
// throw an exception // throw an exception
} }
} }
} }
return is; return is;
} }
} }