Updated multiprecision module to support the most recent version of MPFR C++

unsupported/Eigen/MPRealSupport

@@ -1,21 +1,17 @@
 // This file is part of a joint effort between Eigen, a lightweight C++ template library
 // for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/)
 //
-// Copyright (C) 2010 Pavel Holoborodko <pavel@holoborodko.com>
+// Copyright (C) 2010-2012 Pavel Holoborodko <pavel@holoborodko.com>
 // Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com>
 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-// 
-// Contributors:
-// Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, Pere Constans
 
 #ifndef EIGEN_MPREALSUPPORT_MODULE_H
 #define EIGEN_MPREALSUPPORT_MODULE_H
 
-#include <ctime>
 #include <mpreal.h>
 #include <Eigen/Core>
 
@@ -78,71 +74,65 @@ int main()
     typedef mpfr::mpreal Real;
     typedef mpfr::mpreal NonInteger;
     
-    inline static mpfr::mpreal highest() { return  mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); }
+    inline static Real highest   (long Precision = mpfr::mpreal::get_default_prec())  { return  mpfr::maxval(Precision); }
-    inline static mpfr::mpreal lowest()  { return -mpfr::mpreal_max(mpfr::mpreal::get_default_prec()); }
+    inline static Real lowest    (long Precision = mpfr::mpreal::get_default_prec())  { return -mpfr::maxval(Precision); }
+
+    // Constants
+    inline static Real Pi       (long Precision = mpfr::mpreal::get_default_prec())     {    return mpfr::const_pi(Precision);        }
+    inline static Real Euler    (long Precision = mpfr::mpreal::get_default_prec())     {    return mpfr::const_euler(Precision);     }
+    inline static Real Log2     (long Precision = mpfr::mpreal::get_default_prec())     {    return mpfr::const_log2(Precision);      }
+    inline static Real Catalan  (long Precision = mpfr::mpreal::get_default_prec())     {    return mpfr::const_catalan(Precision);   }
+
+    inline static Real epsilon  (long Precision = mpfr::mpreal::get_default_prec())     {    return mpfr::machine_epsilon(Precision); }
+    inline static Real epsilon  (const Real& x)                                         {    return mpfr::machine_epsilon(x); }
 
-    inline static Real epsilon()
-    {
-      return mpfr::machine_epsilon(mpfr::mpreal::get_default_prec());
-    }
     inline static Real dummy_precision()   
     { 
-      unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1)*90)/100;
+        unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100;
-      return mpfr::machine_epsilon(weak_prec);
+        return mpfr::machine_epsilon(weak_prec);
     }
   };
 
-namespace internal {
+  namespace internal {
 
-  template<> mpfr::mpreal random<mpfr::mpreal>()
+  template<> inline mpfr::mpreal random<mpfr::mpreal>()
   {
-#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+    return mpfr::random();
-    static gmp_randstate_t state;
-    static bool isFirstTime = true;
-
-    if(isFirstTime)
-    {
-      gmp_randinit_default(state);
-      gmp_randseed_ui(state,(unsigned)time(NULL));
-      isFirstTime = false;
-    }
-
-    return mpfr::urandom(state)*2-1;
-#else
-    return mpfr::mpreal(random<double>());
-#endif
   }
 
-  template<> mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
+  template<> inline mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b)
   {
     return a + (b-a) * random<mpfr::mpreal>();
   }
 
-  bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec)
+  inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
   {
-    return mpfr::abs(a) <= mpfr::abs(b) * prec;
+    return mpfr::abs(a) <= mpfr::abs(b) * eps;
   }
 
-  inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec)
+  inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
   {
-    return mpfr::abs(a - b) <= (mpfr::min)(mpfr::abs(a), mpfr::abs(b)) * prec;
+    return mpfr::isEqualFuzzy(a,b,eps);
   }
 
-  inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& prec)
+  inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps)
   {
-    return a <= b || isApprox(a, b, prec);
+    return a <= b || mpfr::isEqualFuzzy(a,b,eps);
   }
 
   template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x)
   { return x.toLDouble(); }
+
   template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x)
   { return x.toDouble(); }
+
   template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x)
   { return x.toLong(); }
+
   template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x)
   { return int(x.toLong()); }
 
-} // end namespace internal
+  } // end namespace internal
 }
 
 #endif // EIGEN_MPREALSUPPORT_MODULE_H

unsupported/test/mpreal/dlmalloc.h

@@ -1,562 +0,0 @@
-/*
-  Default header file for malloc-2.8.x, written by Doug Lea
-  and released to the public domain, as explained at
-  http://creativecommons.org/licenses/publicdomain. 
- 
-  last update: Wed May 27 14:25:17 2009  Doug Lea  (dl at gee)
-
-  This header is for ANSI C/C++ only.  You can set any of
-  the following #defines before including:
-
-  * If USE_DL_PREFIX is defined, it is assumed that malloc.c 
-    was also compiled with this option, so all routines
-    have names starting with "dl".
-
-  * If HAVE_USR_INCLUDE_MALLOC_H is defined, it is assumed that this
-    file will be #included AFTER <malloc.h>. This is needed only if
-    your system defines a struct mallinfo that is incompatible with the
-    standard one declared here.  Otherwise, you can include this file
-    INSTEAD of your system system <malloc.h>.  At least on ANSI, all
-    declarations should be compatible with system versions
-
-  * If MSPACES is defined, declarations for mspace versions are included.
-*/
-
-#ifndef MALLOC_280_H
-#define MALLOC_280_H
-
-#define USE_DL_PREFIX
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#include <stddef.h>   /* for size_t */
-
-#ifndef ONLY_MSPACES
-#define ONLY_MSPACES 0     /* define to a value */
-#endif  /* ONLY_MSPACES */
-#ifndef NO_MALLINFO
-#define NO_MALLINFO 0
-#endif  /* NO_MALLINFO */
-
-
-#if !ONLY_MSPACES
-
-#ifndef USE_DL_PREFIX
-#define dlcalloc               calloc
-#define dlfree                 free
-#define dlmalloc               malloc
-#define dlmemalign             memalign
-#define dlrealloc              realloc
-#define dlvalloc               valloc
-#define dlpvalloc              pvalloc
-#define dlmallinfo             mallinfo
-#define dlmallopt              mallopt
-#define dlmalloc_trim          malloc_trim
-#define dlmalloc_stats         malloc_stats
-#define dlmalloc_usable_size   malloc_usable_size
-#define dlmalloc_footprint     malloc_footprint
-#define dlindependent_calloc   independent_calloc
-#define dlindependent_comalloc independent_comalloc
-#endif /* USE_DL_PREFIX */
-#if !NO_MALLINFO 
-#ifndef HAVE_USR_INCLUDE_MALLOC_H
-#ifndef _MALLOC_H
-#ifndef MALLINFO_FIELD_TYPE
-#define MALLINFO_FIELD_TYPE size_t
-#endif /* MALLINFO_FIELD_TYPE */
-#ifndef STRUCT_MALLINFO_DECLARED
-#define STRUCT_MALLINFO_DECLARED 1
-struct mallinfo {
-  MALLINFO_FIELD_TYPE arena;    /* non-mmapped space allocated from system */
-  MALLINFO_FIELD_TYPE ordblks;  /* number of free chunks */
-  MALLINFO_FIELD_TYPE smblks;   /* always 0 */
-  MALLINFO_FIELD_TYPE hblks;    /* always 0 */
-  MALLINFO_FIELD_TYPE hblkhd;   /* space in mmapped regions */
-  MALLINFO_FIELD_TYPE usmblks;  /* maximum total allocated space */
-  MALLINFO_FIELD_TYPE fsmblks;  /* always 0 */
-  MALLINFO_FIELD_TYPE uordblks; /* total allocated space */
-  MALLINFO_FIELD_TYPE fordblks; /* total free space */
-  MALLINFO_FIELD_TYPE keepcost; /* releasable (via malloc_trim) space */
-};
-#endif /* STRUCT_MALLINFO_DECLARED */
-#endif  /* _MALLOC_H */
-#endif  /* HAVE_USR_INCLUDE_MALLOC_H */
-#endif  /* !NO_MALLINFO */
-
-/*
-  malloc(size_t n)
-  Returns a pointer to a newly allocated chunk of at least n bytes, or
-  null if no space is available, in which case errno is set to ENOMEM
-  on ANSI C systems.
-
-  If n is zero, malloc returns a minimum-sized chunk. (The minimum
-  size is 16 bytes on most 32bit systems, and 32 bytes on 64bit
-  systems.)  Note that size_t is an unsigned type, so calls with
-  arguments that would be negative if signed are interpreted as
-  requests for huge amounts of space, which will often fail. The
-  maximum supported value of n differs across systems, but is in all
-  cases less than the maximum representable value of a size_t.
-*/
-void* dlmalloc(size_t);
-
-/*
-  free(void* p)
-  Releases the chunk of memory pointed to by p, that had been previously
-  allocated using malloc or a related routine such as realloc.
-  It has no effect if p is null. If p was not malloced or already
-  freed, free(p) will by default cuase the current program to abort.
-*/
-void  dlfree(void*);
-
-/*
-  calloc(size_t n_elements, size_t element_size);
-  Returns a pointer to n_elements * element_size bytes, with all locations
-  set to zero.
-*/
-void* dlcalloc(size_t, size_t);
-
-/*
-  realloc(void* p, size_t n)
-  Returns a pointer to a chunk of size n that contains the same data
-  as does chunk p up to the minimum of (n, p's size) bytes, or null
-  if no space is available.
-
-  The returned pointer may or may not be the same as p. The algorithm
-  prefers extending p in most cases when possible, otherwise it
-  employs the equivalent of a malloc-copy-free sequence.
-
-  If p is null, realloc is equivalent to malloc.
-
-  If space is not available, realloc returns null, errno is set (if on
-  ANSI) and p is NOT freed.
-
-  if n is for fewer bytes than already held by p, the newly unused
-  space is lopped off and freed if possible.  realloc with a size
-  argument of zero (re)allocates a minimum-sized chunk.
-
-  The old unix realloc convention of allowing the last-free'd chunk
-  to be used as an argument to realloc is not supported.
-*/
-
-void* dlrealloc(void*, size_t);
-
-/*
-  memalign(size_t alignment, size_t n);
-  Returns a pointer to a newly allocated chunk of n bytes, aligned
-  in accord with the alignment argument.
-
-  The alignment argument should be a power of two. If the argument is
-  not a power of two, the nearest greater power is used.
-  8-byte alignment is guaranteed by normal malloc calls, so don't
-  bother calling memalign with an argument of 8 or less.
-
-  Overreliance on memalign is a sure way to fragment space.
-*/
-void* dlmemalign(size_t, size_t);
-
-/*
-  valloc(size_t n);
-  Equivalent to memalign(pagesize, n), where pagesize is the page
-  size of the system. If the pagesize is unknown, 4096 is used.
-*/
-void* dlvalloc(size_t);
-
-/*
-  mallopt(int parameter_number, int parameter_value)
-  Sets tunable parameters The format is to provide a
-  (parameter-number, parameter-value) pair.  mallopt then sets the
-  corresponding parameter to the argument value if it can (i.e., so
-  long as the value is meaningful), and returns 1 if successful else
-  0.  SVID/XPG/ANSI defines four standard param numbers for mallopt,
-  normally defined in malloc.h.  None of these are use in this malloc,
-  so setting them has no effect. But this malloc also supports other
-  options in mallopt:
-
-  Symbol            param #  default    allowed param values
-  M_TRIM_THRESHOLD     -1   2*1024*1024   any   (-1U disables trimming)
-  M_GRANULARITY        -2     page size   any power of 2 >= page size
-  M_MMAP_THRESHOLD     -3      256*1024   any   (or 0 if no MMAP support)
-*/
-int dlmallopt(int, int);
-
-#define M_TRIM_THRESHOLD     (-1)
-#define M_GRANULARITY        (-2)
-#define M_MMAP_THRESHOLD     (-3)
-
-
-/*
-  malloc_footprint();
-  Returns the number of bytes obtained from the system.  The total
-  number of bytes allocated by malloc, realloc etc., is less than this
-  value. Unlike mallinfo, this function returns only a precomputed
-  result, so can be called frequently to monitor memory consumption.
-  Even if locks are otherwise defined, this function does not use them,
-  so results might not be up to date.
-*/
-size_t dlmalloc_footprint();
-
-#if !NO_MALLINFO
-/*
-  mallinfo()
-  Returns (by copy) a struct containing various summary statistics:
-
-  arena:     current total non-mmapped bytes allocated from system
-  ordblks:   the number of free chunks
-  smblks:    always zero.
-  hblks:     current number of mmapped regions
-  hblkhd:    total bytes held in mmapped regions
-  usmblks:   the maximum total allocated space. This will be greater
-                than current total if trimming has occurred.
-  fsmblks:   always zero
-  uordblks:  current total allocated space (normal or mmapped)
-  fordblks:  total free space
-  keepcost:  the maximum number of bytes that could ideally be released
-               back to system via malloc_trim. ("ideally" means that
-               it ignores page restrictions etc.)
-
-  Because these fields are ints, but internal bookkeeping may
-  be kept as longs, the reported values may wrap around zero and
-  thus be inaccurate.
-*/
-
-struct mallinfo dlmallinfo(void);
-#endif  /* NO_MALLINFO */
-
-/*
-  independent_calloc(size_t n_elements, size_t element_size, void* chunks[]);
-
-  independent_calloc is similar to calloc, but instead of returning a
-  single cleared space, it returns an array of pointers to n_elements
-  independent elements that can hold contents of size elem_size, each
-  of which starts out cleared, and can be independently freed,
-  realloc'ed etc. The elements are guaranteed to be adjacently
-  allocated (this is not guaranteed to occur with multiple callocs or
-  mallocs), which may also improve cache locality in some
-  applications.
-
-  The "chunks" argument is optional (i.e., may be null, which is
-  probably the most typical usage). If it is null, the returned array
-  is itself dynamically allocated and should also be freed when it is
-  no longer needed. Otherwise, the chunks array must be of at least
-  n_elements in length. It is filled in with the pointers to the
-  chunks.
-
-  In either case, independent_calloc returns this pointer array, or
-  null if the allocation failed.  If n_elements is zero and "chunks"
-  is null, it returns a chunk representing an array with zero elements
-  (which should be freed if not wanted).
-
-  Each element must be individually freed when it is no longer
-  needed. If you'd like to instead be able to free all at once, you
-  should instead use regular calloc and assign pointers into this
-  space to represent elements.  (In this case though, you cannot
-  independently free elements.)
-
-  independent_calloc simplifies and speeds up implementations of many
-  kinds of pools.  It may also be useful when constructing large data
-  structures that initially have a fixed number of fixed-sized nodes,
-  but the number is not known at compile time, and some of the nodes
-  may later need to be freed. For example:
-
-  struct Node { int item; struct Node* next; };
-
-  struct Node* build_list() {
-    struct Node** pool;
-    int n = read_number_of_nodes_needed();
-    if (n <= 0) return 0;
-    pool = (struct Node**)(independent_calloc(n, sizeof(struct Node), 0);
-    if (pool == 0) die();
-    // organize into a linked list...
-    struct Node* first = pool[0];
-    for (i = 0; i < n-1; ++i)
-      pool[i]->next = pool[i+1];
-    free(pool);     // Can now free the array (or not, if it is needed later)
-    return first;
-  }
-*/
-void** dlindependent_calloc(size_t, size_t, void**);
-
-/*
-  independent_comalloc(size_t n_elements, size_t sizes[], void* chunks[]);
-
-  independent_comalloc allocates, all at once, a set of n_elements
-  chunks with sizes indicated in the "sizes" array.    It returns
-  an array of pointers to these elements, each of which can be
-  independently freed, realloc'ed etc. The elements are guaranteed to
-  be adjacently allocated (this is not guaranteed to occur with
-  multiple callocs or mallocs), which may also improve cache locality
-  in some applications.
-
-  The "chunks" argument is optional (i.e., may be null). If it is null
-  the returned array is itself dynamically allocated and should also
-  be freed when it is no longer needed. Otherwise, the chunks array
-  must be of at least n_elements in length. It is filled in with the
-  pointers to the chunks.
-
-  In either case, independent_comalloc returns this pointer array, or
-  null if the allocation failed.  If n_elements is zero and chunks is
-  null, it returns a chunk representing an array with zero elements
-  (which should be freed if not wanted).
-
-  Each element must be individually freed when it is no longer
-  needed. If you'd like to instead be able to free all at once, you
-  should instead use a single regular malloc, and assign pointers at
-  particular offsets in the aggregate space. (In this case though, you
-  cannot independently free elements.)
-
-  independent_comallac differs from independent_calloc in that each
-  element may have a different size, and also that it does not
-  automatically clear elements.
-
-  independent_comalloc can be used to speed up allocation in cases
-  where several structs or objects must always be allocated at the
-  same time.  For example:
-
-  struct Head { ... }
-  struct Foot { ... }
-
-  void send_message(char* msg) {
-    int msglen = strlen(msg);
-    size_t sizes[3] = { sizeof(struct Head), msglen, sizeof(struct Foot) };
-    void* chunks[3];
-    if (independent_comalloc(3, sizes, chunks) == 0)
-      die();
-    struct Head* head = (struct Head*)(chunks[0]);
-    char*        body = (char*)(chunks[1]);
-    struct Foot* foot = (struct Foot*)(chunks[2]);
-    // ...
-  }
-
-  In general though, independent_comalloc is worth using only for
-  larger values of n_elements. For small values, you probably won't
-  detect enough difference from series of malloc calls to bother.
-
-  Overuse of independent_comalloc can increase overall memory usage,
-  since it cannot reuse existing noncontiguous small chunks that
-  might be available for some of the elements.
-*/
-void** dlindependent_comalloc(size_t, size_t*, void**);
-
-
-/*
-  pvalloc(size_t n);
-  Equivalent to valloc(minimum-page-that-holds(n)), that is,
-  round up n to nearest pagesize.
- */
-void*  dlpvalloc(size_t);
-
-/*
-  malloc_trim(size_t pad);
-
-  If possible, gives memory back to the system (via negative arguments
-  to sbrk) if there is unused memory at the `high' end of the malloc
-  pool or in unused MMAP segments. You can call this after freeing
-  large blocks of memory to potentially reduce the system-level memory
-  requirements of a program. However, it cannot guarantee to reduce
-  memory. Under some allocation patterns, some large free blocks of
-  memory will be locked between two used chunks, so they cannot be
-  given back to the system.
-
-  The `pad' argument to malloc_trim represents the amount of free
-  trailing space to leave untrimmed. If this argument is zero, only
-  the minimum amount of memory to maintain internal data structures
-  will be left. Non-zero arguments can be supplied to maintain enough
-  trailing space to service future expected allocations without having
-  to re-obtain memory from the system.
-
-  Malloc_trim returns 1 if it actually released any memory, else 0.
-*/
-int  dlmalloc_trim(size_t);
-
-/*
-  malloc_stats();
-  Prints on stderr the amount of space obtained from the system (both
-  via sbrk and mmap), the maximum amount (which may be more than
-  current if malloc_trim and/or munmap got called), and the current
-  number of bytes allocated via malloc (or realloc, etc) but not yet
-  freed. Note that this is the number of bytes allocated, not the
-  number requested. It will be larger than the number requested
-  because of alignment and bookkeeping overhead. Because it includes
-  alignment wastage as being in use, this figure may be greater than
-  zero even when no user-level chunks are allocated.
-
-  The reported current and maximum system memory can be inaccurate if
-  a program makes other calls to system memory allocation functions
-  (normally sbrk) outside of malloc.
-
-  malloc_stats prints only the most commonly interesting statistics.
-  More information can be obtained by calling mallinfo.
-*/
-void  dlmalloc_stats();
-
-#endif /* !ONLY_MSPACES */
-
-/*
-  malloc_usable_size(void* p);
-
-  Returns the number of bytes you can actually use in
-  an allocated chunk, which may be more than you requested (although
-  often not) due to alignment and minimum size constraints.
-  You can use this many bytes without worrying about
-  overwriting other allocated objects. This is not a particularly great
-  programming practice. malloc_usable_size can be more useful in
-  debugging and assertions, for example:
-
-  p = malloc(n);
-  assert(malloc_usable_size(p) >= 256);
-*/
-size_t dlmalloc_usable_size(void*);
-
-
-#if MSPACES
-
-/*
-  mspace is an opaque type representing an independent
-  region of space that supports mspace_malloc, etc.
-*/
-typedef void* mspace;
-
-/*
-  create_mspace creates and returns a new independent space with the
-  given initial capacity, or, if 0, the default granularity size.  It
-  returns null if there is no system memory available to create the
-  space.  If argument locked is non-zero, the space uses a separate
-  lock to control access. The capacity of the space will grow
-  dynamically as needed to service mspace_malloc requests.  You can
-  control the sizes of incremental increases of this space by
-  compiling with a different DEFAULT_GRANULARITY or dynamically
-  setting with mallopt(M_GRANULARITY, value).
-*/
-mspace create_mspace(size_t capacity, int locked);
-
-/*
-  destroy_mspace destroys the given space, and attempts to return all
-  of its memory back to the system, returning the total number of
-  bytes freed. After destruction, the results of access to all memory
-  used by the space become undefined.
-*/
-size_t destroy_mspace(mspace msp);
-
-/*
-  create_mspace_with_base uses the memory supplied as the initial base
-  of a new mspace. Part (less than 128*sizeof(size_t) bytes) of this
-  space is used for bookkeeping, so the capacity must be at least this
-  large. (Otherwise 0 is returned.) When this initial space is
-  exhausted, additional memory will be obtained from the system.
-  Destroying this space will deallocate all additionally allocated
-  space (if possible) but not the initial base.
-*/
-mspace create_mspace_with_base(void* base, size_t capacity, int locked);
-
-/*
-  mspace_track_large_chunks controls whether requests for large chunks
-  are allocated in their own untracked mmapped regions, separate from
-  others in this mspace. By default large chunks are not tracked,
-  which reduces fragmentation. However, such chunks are not
-  necessarily released to the system upon destroy_mspace.  Enabling
-  tracking by setting to true may increase fragmentation, but avoids
-  leakage when relying on destroy_mspace to release all memory
-  allocated using this space.  The function returns the previous
-  setting.
-*/
-int mspace_track_large_chunks(mspace msp, int enable);
-
-/*
-  mspace_malloc behaves as malloc, but operates within
-  the given space.
-*/
-void* mspace_malloc(mspace msp, size_t bytes);
-
-/*
-  mspace_free behaves as free, but operates within
-  the given space.
-
-  If compiled with FOOTERS==1, mspace_free is not actually needed.
-  free may be called instead of mspace_free because freed chunks from
-  any space are handled by their originating spaces.
-*/
-void mspace_free(mspace msp, void* mem);
-
-/*
-  mspace_realloc behaves as realloc, but operates within
-  the given space.
-
-  If compiled with FOOTERS==1, mspace_realloc is not actually
-  needed.  realloc may be called instead of mspace_realloc because
-  realloced chunks from any space are handled by their originating
-  spaces.
-*/
-void* mspace_realloc(mspace msp, void* mem, size_t newsize);
-
-/*
-  mspace_calloc behaves as calloc, but operates within
-  the given space.
-*/
-void* mspace_calloc(mspace msp, size_t n_elements, size_t elem_size);
-
-/*
-  mspace_memalign behaves as memalign, but operates within
-  the given space.
-*/
-void* mspace_memalign(mspace msp, size_t alignment, size_t bytes);
-
-/*
-  mspace_independent_calloc behaves as independent_calloc, but
-  operates within the given space.
-*/
-void** mspace_independent_calloc(mspace msp, size_t n_elements,
-                                 size_t elem_size, void* chunks[]);
-
-/*
-  mspace_independent_comalloc behaves as independent_comalloc, but
-  operates within the given space.
-*/
-void** mspace_independent_comalloc(mspace msp, size_t n_elements,
-                                   size_t sizes[], void* chunks[]);
-
-/*
-  mspace_footprint() returns the number of bytes obtained from the
-  system for this space.
-*/
-size_t mspace_footprint(mspace msp);
-
-
-#if !NO_MALLINFO
-/*
-  mspace_mallinfo behaves as mallinfo, but reports properties of
-  the given space.
-*/
-struct mallinfo mspace_mallinfo(mspace msp);
-#endif /* NO_MALLINFO */
-
-/*
-  malloc_usable_size(void* p) behaves the same as malloc_usable_size;
-*/
- size_t mspace_usable_size(void* mem);
-
-/*
-  mspace_malloc_stats behaves as malloc_stats, but reports
-  properties of the given space.
-*/
-void mspace_malloc_stats(mspace msp);
-
-/*
-  mspace_trim behaves as malloc_trim, but
-  operates within the given space.
-*/
-int mspace_trim(mspace msp, size_t pad);
-
-/*
-  An alias for mallopt.
-*/
-int mspace_mallopt(int, int);
-
-#endif  /* MSPACES */
-
-#ifdef __cplusplus
-};  /* end of extern "C" */
-#endif
-
-#endif /* MALLOC_280_H */

unsupported/test/mpreal/mpreal.cpp

@@ -1,597 +0,0 @@
-/*
-	Multi-precision real number class. C++ interface fo MPFR library.
-	Project homepage: http://www.holoborodko.com/pavel/
-	Contact e-mail:   pavel@holoborodko.com
-
-	Copyright (c) 2008-2011 Pavel Holoborodko
-
-	Core Developers: 
-	Pavel Holoborodko, Dmitriy Gubanov, Konstantin Holoborodko. 
-
-	Contributors:
-	Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, 
-	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
-	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
-
-	****************************************************************************
-	****************************************************************************
-	Redistribution and use in source and binary forms, with or without
-	modification, are permitted provided that the following conditions
-	are met:
-
-	1. Redistributions of source code must retain the above copyright
-	notice, this list of conditions and the following disclaimer.
-
-	2. Redistributions in binary form must reproduce the above copyright
-	notice, this list of conditions and the following disclaimer in the
-	documentation and/or other materials provided with the distribution.
-
-	3. The name of the author may be used to endorse or promote products
-	derived from this software without specific prior written permission.
-
-	THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
-	ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-	ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
-	FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-	DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
-	OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
-	HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-	LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
-	OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
-	SUCH DAMAGE.
-*/
-#include <cstring>
-#include "mpreal.h"
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-#include "dlmalloc.h"
-#endif
-
-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_RNDN;	//(mpfr_get_default_rounding_mode)();	
-mp_prec_t  mpreal::default_prec = 64;			//(mpfr_get_default_prec)();	
-int		   mpreal::default_base = 10;
-int        mpreal::double_bits = -1;
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-bool       mpreal::is_custom_malloc = false;
-#endif
-
-// Default constructor: creates mp number and initializes it to 0.
-mpreal::mpreal() 
-{ 
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-	mpfr_init2(mp,default_prec); 
-	mpfr_set_ui(mp,0,default_rnd);
-
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const mpreal& u) 
-{
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-	mpfr_init2(mp,mpfr_get_prec(u.mp));
-	mpfr_set(mp,u.mp,default_rnd);
-
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const mpfr_t u)
-{
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-	mpfr_init2(mp,mpfr_get_prec(u));
-	mpfr_set(mp,u,default_rnd);
-	
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const mpf_t u)
-{
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#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);
-
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const mpz_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_z(mp,u,mode);
-	
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const mpq_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_q(mp,u,mode);
-	
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
-{
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-    if(double_bits == -1 || fits_in_bits(u, double_bits))
-    {
-    	mpfr_init2(mp,prec);
-	    mpfr_set_d(mp,u,mode);
-		
-		MPREAL_MSVC_DEBUGVIEW_CODE;
-    }
-    else
-        throw conversion_overflow();
-}
-
-mpreal::mpreal(const long double 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_ld(mp,u,mode);
-	
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-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();
-#endif
-
-	mpfr_init2(mp,prec);
-	mpfr_set_ui(mp,u,mode);
-	
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const unsigned int 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_ui(mp,u,mode);
-
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const long int 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_si(mp,u,mode);
-
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::mpreal(const int 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_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)
-{
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-	mpfr_init2(mp,prec);
-	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)
-{
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-	mpfr_init2(mp,prec);
-	mpfr_set_str(mp, s.c_str(), base, mode); 
-
-	MPREAL_MSVC_DEBUGVIEW_CODE;
-}
-
-mpreal::~mpreal() 
-{ 
-	mpfr_clear(mp);
-}                           
-
-// Operators - Assignment
-mpreal& mpreal::operator=(const char* s)
-{
-	mpfr_t t;
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-	if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
-	{
-		// We will rewrite mp anyway, so flash it and resize
-		mpfr_set_prec(mp,mpfr_get_prec(t)); 
-		mpfr_set(mp,t,mpreal::default_rnd);
-		mpfr_clear(t);
-
-		MPREAL_MSVC_DEBUGVIEW_CODE;
-
-	}else{
-		mpfr_clear(t);
-	}
-
-	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 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 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 toString(T t, std::ios_base & (*f)(std::ios_base&))
-{
-	std::ostringstream oss;
-	oss << f << t;
-	return oss.str();
-}
-
-#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
-
-std::string mpreal::toString(const std::string& format) const
-{
-	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;
-}
-
-#endif
-
-std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
-{
-  (void)b;
-  (void)mode;
-#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;
-	mp_exp_t exp;
-	string out;
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	set_custom_malloc();
-#endif
-
-	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(NULL,&exp,b,0,mp,mode);
-	ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
-
-	if(s!=NULL && ns!=NULL)
-	{
-		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::toString<mp_exp_t>(exp,std::dec);
-				else 	  out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec);
-			}
-		}
-
-		mpfr_free_str(s);
-		return out;
-	}else{
-		return "conversion error!";
-	}
-#endif
-}
-
-
-//////////////////////////////////////////////////////////////////////////
-// I/O
-ostream& operator<<(ostream& os, const mpreal& v)
-{
-	return os<<v.toString(static_cast<int>(os.precision()));
-}
-
-istream& operator>>(istream &is, mpreal& v)
-{
-	string tmp;
-	is >> tmp;
-	mpfr_set_str(v.mp, tmp.c_str(),mpreal::default_base,mpreal::default_rnd);
-	return is;
-}
-
-
-#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC)
-	// 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)
-	{
-		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
-
-}
-

unsupported/test/mpreal_support.cpp

@@ -57,8 +57,3 @@ void test_mpreal_support()
     stream << A;
   }
 }
-
-extern "C" {
-#include "mpreal/dlmalloc.c"
-}
-#include "mpreal/mpreal.cpp"