update internal mpfr C++ copy

This commit is contained in:
Gael Guennebaud 2012-06-21 09:56:54 +02:00
parent 8c71d7314b
commit b5093e2585
4 changed files with 825 additions and 1216 deletions

View File

@ -1267,7 +1267,7 @@ int mspace_mallopt(int, int);
#endif /* MSPACES */ #endif /* MSPACES */
#ifdef __cplusplus #ifdef __cplusplus
} /* end of extern "C" */ }; /* end of extern "C" */
#endif /* __cplusplus */ #endif /* __cplusplus */
/* /*

View File

@ -556,7 +556,7 @@ int mspace_mallopt(int, int);
#endif /* MSPACES */ #endif /* MSPACES */
#ifdef __cplusplus #ifdef __cplusplus
} /* end of extern "C" */ }; /* end of extern "C" */
#endif #endif
#endif /* MALLOC_280_H */ #endif /* MALLOC_280_H */

View File

@ -3,14 +3,15 @@
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-2011 Pavel Holoborodko
Core Developers: Core Developers:
Pavel Holoborodko, Dmitriy Gubanov, Konstantin Holoborodko. Pavel Holoborodko, Dmitriy Gubanov, Konstantin Holoborodko.
Contributors: Contributors:
Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze,
Heinz van Saanen, Pere Constans, Peter van Hoof. Heinz van Saanen, Pere Constans, Peter van Hoof, Gael Guennebaud,
Tsai Chia Cheng, Alexei Zubanov.
**************************************************************************** ****************************************************************************
This library is free software; you can redistribute it and/or This library is free software; you can redistribute it and/or
@ -27,31 +28,21 @@
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
****************************************************************************
**************************************************************************** ****************************************************************************
Redistribution and use in source and binary forms, with or without Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions modification, are permitted provided that the following conditions
are met: are met:
1. Redistributions of source code must retain the above copyright 1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer. notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright 2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution. documentation and/or other materials provided with the distribution.
3. Redistributions of any form whatsoever must retain the following
acknowledgment:
"
This product includes software developed by Pavel Holoborodko
Web: http://www.holoborodko.com/pavel/
e-mail: pavel@holoborodko.com
"
4. This software cannot be, by any means, used for any commercial 3. The name of the author may be used to endorse or promote products
purpose without the prior permission of the copyright holder. derived from this software without specific prior written permission.
Any of the above conditions can be waived if you get permission from
the copyright holder.
THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
@ -66,9 +57,11 @@
SUCH DAMAGE. SUCH DAMAGE.
*/ */
#include <cstring> #include <cstring>
#include <cstdlib>
#include "mpreal.h" #include "mpreal.h"
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
#include "dlmalloc.h" #include "dlmalloc.h"
#endif
using std::ws; using std::ws;
using std::cerr; using std::cerr;
@ -79,62 +72,107 @@ 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_RNDN; //(mpfr_get_default_rounding_mode)();
mp_prec_t mpreal::default_prec = mpfr_get_default_prec(); mp_prec_t mpreal::default_prec = 64; //(mpfr_get_default_prec)();
int mpreal::default_base = 10; int mpreal::default_base = 10;
int mpreal::double_bits = -1; int mpreal::double_bits = -1;
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
bool mpreal::is_custom_malloc = false; bool mpreal::is_custom_malloc = false;
#endif
// Default constructor: creates mp number and initializes it to 0. // Default constructor: creates mp number and initializes it to 0.
mpreal::mpreal() mpreal::mpreal()
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
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_MSVC_DEBUGVIEW_CODE;
} }
mpreal::mpreal(const mpreal& u) mpreal::mpreal(const mpreal& u)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
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_MSVC_DEBUGVIEW_CODE;
} }
mpreal::mpreal(const mpfr_t u) mpreal::mpreal(const mpfr_t u)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
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_MSVC_DEBUGVIEW_CODE;
} }
mpreal::mpreal(const mpf_t u) mpreal::mpreal(const mpf_t u)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
mpfr_init2(mp,mpf_get_prec(u)); #endif
mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
mpfr_set_f(mp,u,default_rnd); mpfr_set_f(mp,u,default_rnd);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_z(mp,u,mode); mpfr_set_z(mp,u,mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_q(mp,u,mode); mpfr_set_q(mp,u,mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
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 defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
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);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
else else
throw conversion_overflow(); throw conversion_overflow();
@ -142,51 +180,121 @@ mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_ld(mp,u,mode); mpfr_set_ld(mp,u,mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_ui(mp,u,mode); mpfr_set_ui(mp,u,mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_ui(mp,u,mode); mpfr_set_ui(mp,u,mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_si(mp,u,mode); mpfr_set_si(mp,u,mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_si(mp,u,mode); mpfr_set_si(mp,u,mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
#if defined (MPREAL_HAVE_INT64_SUPPORT)
mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode)
{
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc();
#endif
mpfr_init2(mp,prec);
mpfr_set_uj(mp, u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)
{
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc();
#endif
mpfr_init2(mp,prec);
mpfr_set_sj(mp, u, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
}
#endif
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)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_str(mp, s, base, mode); mpfr_set_str(mp, s, base, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode) mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)
{ {
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
mpfr_init2(mp,prec); mpfr_init2(mp,prec);
mpfr_set_str(mp, s.c_str(), base, mode); mpfr_set_str(mp, s.c_str(), base, mode);
MPREAL_MSVC_DEBUGVIEW_CODE;
} }
mpreal::~mpreal() mpreal::~mpreal()
@ -198,18 +306,22 @@ mpreal::~mpreal()
mpreal& mpreal::operator=(const char* s) mpreal& mpreal::operator=(const char* s)
{ {
mpfr_t t; mpfr_t t;
#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
if(0==mpfr_init_set_str(t,s,default_base,default_rnd)) if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
{ {
// We will rewrite mp anyway, so use flash it and resize // We will rewrite mp anyway, so flash it and resize
mpfr_set_prec(mp,mpfr_get_prec(t)); //<- added 01.04.2011 mpfr_set_prec(mp,mpfr_get_prec(t));
mpfr_set(mp,t,mpreal::default_rnd); mpfr_set(mp,t,mpreal::default_rnd);
mpfr_clear(t); mpfr_clear(t);
MPREAL_MSVC_DEBUGVIEW_CODE;
}else{ }else{
mpfr_clear(t); mpfr_clear(t);
// cerr<<"fail to convert string"<<endl;
} }
return *this; return *this;
@ -260,21 +372,6 @@ const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
return a; 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) const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
{ {
mpreal x; mpreal x;
@ -288,21 +385,6 @@ const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
return x; 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) const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ {
mpreal a; mpreal a;
@ -319,36 +401,70 @@ const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mod
} }
template <class T> template <class T>
std::string to_string(T t, std::ios_base & (*f)(std::ios_base&)) std::string toString(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 #if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
std::string mpreal::toString(const std::string& format) const
{ {
return to_string(); char *s = NULL;
string out;
if( !format.empty() )
{
if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0))
{
out = std::string(s);
mpfr_free_str(s);
}
}
return out;
} }
std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const #endif
std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
{ {
char *s, *ns = NULL;
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
// Use MPFR native function for output
char format[128];
int digits;
digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp));
sprintf(format,"%%.%dRNg",digits); // Default format
return toString(std::string(format));
#else
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 defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
set_custom_malloc(); set_custom_malloc();
#endif
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);
@ -419,8 +535,8 @@ std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
// 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::toString<mp_exp_t>(exp,std::dec);
else out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec); else out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec);
} }
} }
@ -429,79 +545,52 @@ std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
}else{ }else{
return "conversion error!"; return "conversion error!";
} }
#endif
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
// I/O // I/O
ostream& operator<<(ostream& os, const mpreal& v) ostream& operator<<(ostream& os, const mpreal& v)
{ {
return os<<v.to_string(static_cast<size_t>(os.precision())); return os<<v.toString(static_cast<int>(os.precision()));
} }
istream& operator>>(istream &is, mpreal& v) istream& operator>>(istream &is, mpreal& v)
{ {
char c; string tmp;
string s = ""; is >> tmp;
mpfr_t t; mpfr_set_str(v.mp, tmp.c_str(),mpreal::default_base,mpreal::default_rnd);
mpreal::set_custom_malloc();
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; return is;
} }
// Optimized dynamic memory allocation/(re-)deallocation.
void * mpreal::mpreal_allocate(size_t alloc_size)
{
return(dlmalloc(alloc_size));
}
void * mpreal::mpreal_reallocate(void *ptr, size_t /*old_size*/, size_t new_size) #if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
{ // Optimized dynamic memory allocation/(re-)deallocation.
return(dlrealloc(ptr,new_size)); void * mpreal::mpreal_allocate(size_t alloc_size)
}
void mpreal::mpreal_free(void *ptr, size_t /*size*/)
{
dlfree(ptr);
}
inline void mpreal::set_custom_malloc(void)
{
if(!is_custom_malloc)
{ {
mp_set_memory_functions(mpreal_allocate,mpreal_reallocate,mpreal_free); return(dlmalloc(alloc_size));
is_custom_malloc = true;
} }
}
void * mpreal::mpreal_reallocate(void *ptr, size_t old_size, size_t new_size)
{
return(dlrealloc(ptr,new_size));
}
void mpreal::mpreal_free(void *ptr, size_t size)
{
dlfree(ptr);
}
inline void mpreal::set_custom_malloc(void)
{
if(!is_custom_malloc)
{
mp_set_memory_functions(mpreal_allocate,mpreal_reallocate,mpreal_free);
is_custom_malloc = true;
}
}
#endif
} }

File diff suppressed because it is too large Load Diff