update to latest mpreal and fix a min/max issue in mprel.h

unsupported/test/dlmalloc.h

@@ -0,0 +1,562 @@
+/*
+  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.cpp

@@ -1,10 +1,18 @@
 /*
-	Multi-precision real number class. C++ wrapper fo MPFR library.
+	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-2010 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.
+
+	****************************************************************************
 	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
@@ -19,13 +27,48 @@
 	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:
+	****************************************************************************
-	Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, 
+	Redistribution and use in source and binary forms, with or without
-	Heinz van Saanen, Pere Constans, Dmitriy Gubanov
+	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. 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 
+	purpose without the prior permission of the copyright holder.
+	
+	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
+	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 <cstdlib>
 #include "mpreal.h"
+#include "dlmalloc.h"
 
 using std::ws;
 using std::cerr;
@@ -40,46 +83,54 @@ 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;
+bool       mpreal::is_custom_malloc = false;
 
 // Default constructor: creates mp number and initializes it to 0.
 mpreal::mpreal() 
 { 
+	set_custom_malloc();
 	mpfr_init2(mp,default_prec); 
 	mpfr_set_ui(mp,0,default_rnd);
 }
 
 mpreal::mpreal(const mpreal& u) 
 {
+	set_custom_malloc();
 	mpfr_init2(mp,mpfr_get_prec(u.mp));
 	mpfr_set(mp,u.mp,default_rnd);
 }
 
 mpreal::mpreal(const mpfr_t u)
 {
+	set_custom_malloc();
 	mpfr_init2(mp,mpfr_get_prec(u));
 	mpfr_set(mp,u,default_rnd);
 }
 
 mpreal::mpreal(const mpf_t u)
 {
+	set_custom_malloc();
 	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)
 {
+	set_custom_malloc();
 	mpfr_init2(mp,prec);
 	mpfr_set_z(mp,u,mode);
 }
 
 mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
 {
+	set_custom_malloc();
 	mpfr_init2(mp,prec);
 	mpfr_set_q(mp,u,mode);
 }
 
 mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
 {
+	set_custom_malloc();
     if(double_bits == -1 || fits_in_bits(u, double_bits))
     {
     	mpfr_init2(mp,prec);
@@ -91,40 +142,53 @@ 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)
 { 
+	set_custom_malloc();
     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)
 { 
+	set_custom_malloc();
 	mpfr_init2(mp,prec);
 	mpfr_set_ui(mp,u,mode);
 }
 
 mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
 { 
+	set_custom_malloc();
 	mpfr_init2(mp,prec);
 	mpfr_set_ui(mp,u,mode);
 }
 
 mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
 { 
+	set_custom_malloc();
 	mpfr_init2(mp,prec);
 	mpfr_set_si(mp,u,mode);
 }
 
 mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
 { 
+	set_custom_malloc();
 	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)
 {
+	set_custom_malloc();
 	mpfr_init2(mp,prec);
 	mpfr_set_str(mp, s, base, mode); 
 }
 
+mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)
+{
+	set_custom_malloc();
+	mpfr_init2(mp,prec);
+	mpfr_set_str(mp, s.c_str(), base, mode); 
+}
+
 mpreal::~mpreal() 
 { 
 	mpfr_clear(mp);
@@ -134,8 +198,13 @@ mpreal::~mpreal()
 mpreal& mpreal::operator=(const char* s)
 {
 	mpfr_t t;
+	
+	set_custom_malloc();
+
 	if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
 	{
+		// We will rewrite mp anyway, so use flash it and resize
+		mpfr_set_prec(mp,mpfr_get_prec(t)); //<- added 01.04.2011
 		mpfr_set(mp,t,mpreal::default_rnd);
 		mpfr_clear(t);
 	}else{
@@ -262,13 +331,15 @@ mpreal::operator std::string() const
 	return to_string();
 }
 
-string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
+std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
 {
 	char *s, *ns = NULL;	
 	size_t slen, nslen;
 	mp_exp_t exp;
 	string out;
 	
+	set_custom_malloc();
+	
 	if(mpfr_inf_p(mp))
 	{ 
 		if(mpfr_sgn(mp)>0) return "+@Inf@";
@@ -278,7 +349,6 @@ string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
 	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);
 
@@ -365,7 +435,7 @@ string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
 // I/O
 ostream& operator<<(ostream& os, const mpreal& v)
 {
-	return os<<v.to_string(os.precision());
+	return os<<v.to_string(static_cast<size_t>(os.precision()));
 }
 
 istream& operator>>(istream &is, mpreal& v)
@@ -374,6 +444,8 @@ istream& operator>>(istream &is, mpreal& v)
 	string s = "";
 	mpfr_t t;
 	
+	mpreal::set_custom_malloc();
+	
 	if(is.good())
 	{
 		is>>ws;
@@ -391,6 +463,7 @@ istream& operator>>(istream &is, mpreal& v)
 		{
 			// 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);				
@@ -405,4 +478,30 @@ istream& operator>>(istream &is, mpreal& v)
 	}
 	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)
+{
+	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;
+	}
+}
+}
+

unsupported/test/mpreal.h

@@ -1,10 +1,18 @@
 /*
-	Multi-precision real number class. C++ wrapper fo MPFR library.
+	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-2010 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.
+
+	****************************************************************************
 	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
@@ -19,9 +27,43 @@
 	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:
+	****************************************************************************
-	Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, 
+	Redistribution and use in source and binary forms, with or without
-	Heinz van Saanen, Pere Constans, Dmitriy Gubanov
+	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. 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 
+	purpose without the prior permission of the copyright holder.
+	
+	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
+	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.
 */
 
 #ifndef __MP_REAL_H__
@@ -38,7 +80,7 @@
 
 // Detect compiler using signatures from http://predef.sourceforge.net/
 #if defined(__GNUC__) && defined(__INTEL_COMPILER)
-	#define IsInf(x) isinf(x)                                // GNU C/C++ + Intel ICC compiler
+	#define IsInf(x) isinf(x)								// Intel ICC compiler on Linux 
 
 #elif defined(__GNUC__)
 	#define IsInf(x) std::isinf(x)							// GNU C/C++ 
@@ -47,7 +89,7 @@
 	#define IsInf(x) (!_finite(x))							// Microsoft Visual C++
 
 #else
-	#define IsInf(x) std::isinf(x)				// C99 conformance
+	#define IsInf(x) std::isinf(x)							// Unknown compiler, just hope for C99 conformance
 #endif
 
 namespace mpfr {
@@ -79,6 +121,7 @@ public:
 	mpreal(const long int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
 	mpreal(const int u, mp_prec_t prec = default_prec, mp_rnd_t mode = default_rnd);
 	mpreal(const char* s, mp_prec_t prec = default_prec, int base = default_base, mp_rnd_t mode = default_rnd);
+	mpreal(const std::string& s, mp_prec_t prec = default_prec, int base = default_base, mp_rnd_t mode = default_rnd);
 
 	~mpreal();                           
 
@@ -189,6 +232,7 @@ public:
 	inline operator unsigned long() const;
 	inline operator unsigned int() const;
 	inline operator long() const;
+	inline operator int() const;
 	operator std::string() const;
 	inline operator mpfr_ptr();
 
@@ -205,6 +249,7 @@ public:
 	friend const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
+
 	friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode = mpreal::default_rnd);
@@ -219,6 +264,8 @@ public:
 	friend const mpreal exp  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd); 
 	friend const mpreal exp2 (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
+	friend const mpreal log1p  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
+	friend const mpreal expm1  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 
 	friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
@@ -238,18 +285,17 @@ public:
 	friend const mpreal sech  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal csch  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal coth  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
-
 	friend const mpreal acosh  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal asinh  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal atanh  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
-	friend const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode = mpreal::default_rnd);
+	friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd);
-	friend const mpreal log1p  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
+
-	friend const mpreal expm1  (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
+	friend const mpreal fac_ui (unsigned long int v,  mp_prec_t prec = mpreal::default_prec, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal eint   (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 
 	friend const mpreal gamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
-	friend const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode = mpreal::default_rnd);
+	friend const mpreal lgamma (const mpreal& v, int *signp = 0, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal zeta (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal erf (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal erfc (const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
@@ -262,7 +308,6 @@ public:
 	friend const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::default_rnd);
-	friend const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode = mpreal::default_rnd);
 	friend int sgn(const mpreal& v); // -1 or +1
 
@@ -277,6 +322,7 @@ public:
 // MPFR 3.0.0 Specifics
 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
 	friend const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
+	friend const mpreal ai(const mpreal& v, mp_rnd_t rnd_mode = mpreal::default_rnd);
     friend const mpreal urandom (gmp_randstate_t& state,mp_rnd_t rnd_mode = mpreal::default_rnd); 	// use gmp_randinit_default() to init state, gmp_randclear() to clear
 	friend bool _isregular(const mpreal& v);
 #endif
@@ -391,6 +437,17 @@ public:
 	#ifndef min
 		friend const mpreal min(const mpreal& x, const mpreal& y);
 	#endif
+	
+	friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = default_rnd);
+	friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = default_rnd);
+
+private:
+	// Optimized dynamic memory allocation/(re-)deallocation.
+	static bool is_custom_malloc;
+	static void *mpreal_allocate			(size_t alloc_size);
+	static void *mpreal_reallocate			(void *ptr, size_t old_size, size_t new_size);
+	static void mpreal_free					(void *ptr, size_t size);
+	inline static void set_custom_malloc	(void);
 };
 
 //////////////////////////////////////////////////////////////////////////
@@ -647,9 +704,9 @@ const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::defaul
 
 //////////////////////////////////////////////////////////////////////////
 // Estimate machine epsilon for the given precision
-inline const mpreal machine_epsilon(mp_prec_t prec);
+inline const mpreal machine_epsilon(mp_prec_t prec = mpreal::default_prec);
-inline const mpreal mpreal_min(mp_prec_t prec);
+inline const mpreal mpreal_min(mp_prec_t prec = mpreal::default_prec);
-inline const mpreal mpreal_max(mp_prec_t prec);
+inline const mpreal mpreal_max(mp_prec_t prec = mpreal::default_prec);
 
 //////////////////////////////////////////////////////////////////////////
 // Implementation of inline functions
@@ -659,7 +716,11 @@ inline const mpreal mpreal_max(mp_prec_t prec);
 // Operators - Assignment
 inline mpreal& mpreal::operator=(const mpreal& v)
 {
-	if (this!= &v)	mpfr_set(mp,v.mp,default_rnd);
+	if (this!= &v)
+	{
+		mpfr_set_prec(mp,mpfr_get_prec(v.mp)); 
+		mpfr_set(mp,v.mp,default_rnd);
+	}
 	return *this;
 }
 
@@ -1972,17 +2033,22 @@ inline mpreal::operator long double() const
 
 inline mpreal::operator unsigned long() const
 {
-	return mpfr_get_ui(mp,default_rnd);	
+	return mpfr_get_ui(mp,GMP_RNDZ);	
 }
 
 inline mpreal::operator unsigned int() const
 {
-	return static_cast<unsigned int>(mpfr_get_ui(mp,default_rnd));	
+	return static_cast<unsigned int>(mpfr_get_ui(mp,GMP_RNDZ));	
 }
 
 inline mpreal::operator long() const
 {
-	return mpfr_get_si(mp,default_rnd);	
+	return mpfr_get_si(mp,GMP_RNDZ);	
+}
+
+inline mpreal::operator int() const
+{
+	return static_cast<int>(mpfr_get_si(mp,GMP_RNDZ));	
 }
 
 inline mpreal::operator mpfr_ptr()
@@ -2052,7 +2118,7 @@ inline const mpreal ldexp(const mpreal& v, mp_exp_t exp)
 
 inline const mpreal machine_epsilon(mp_prec_t prec)
 {
-	// smallest eps such that 1.0+eps != 1.0
+	// the smallest eps such that 1.0+eps != 1.0
 	// depends (of cause) on the precision
 	mpreal x(1,prec); 
 	return nextabove(x)-x;
@@ -2415,9 +2481,9 @@ inline const mpreal atanh  (const mpreal& v, mp_rnd_t rnd_mode)
 	return x;
 }
 
-inline const mpreal fac_ui (unsigned long int v, mp_rnd_t rnd_mode)
+inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode)
 {
-	mpreal x;
+	mpreal x(0,prec);
 	mpfr_fac_ui(x.mp,v,rnd_mode);
 	return x;
 }
@@ -2460,7 +2526,13 @@ inline const mpreal lngamma (const mpreal& v, mp_rnd_t rnd_mode)
 inline const mpreal lgamma (const mpreal& v, int *signp, mp_rnd_t rnd_mode)
 {
 	mpreal x(v);
+	int tsignp;
+
+	if(signp)
 		mpfr_lgamma(x.mp,signp,v.mp,rnd_mode);
+	else
+		mpfr_lgamma(x.mp,&tsignp,v.mp,rnd_mode);
+
 	return x;
 }
 
@@ -2569,12 +2641,21 @@ inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode)
 //////////////////////////////////////////////////////////////////////////
 // MPFR 3.0.0 Specifics
 #if (MPFR_VERSION >= MPFR_VERSION_NUM(3,0,0))
+
 inline const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode)
 {
 	mpreal x(v);
 	mpfr_digamma(x.mp,v.mp,rnd_mode);
 	return x;
 }
+
+inline const mpreal ai(const mpreal& v, mp_rnd_t rnd_mode)
+{
+	mpreal x(v);
+	mpfr_ai(x.mp,v.mp,rnd_mode);
+	return x;
+}
+
 #endif // MPFR 3.0.0 Specifics
 
 //////////////////////////////////////////////////////////////////////////
@@ -2699,19 +2780,30 @@ inline void swap(mpreal& a, mpreal& b)
 	mpfr_swap(a.mp,b.mp);
 }
 
-#ifndef max
+
-inline const mpreal max(const mpreal& x, const mpreal& y)
+inline const mpreal (max)(const mpreal& x, const mpreal& y)
 {
 	return (x>y?x:y);
 }
-#endif
 
-#ifndef min
+inline const mpreal (min)(const mpreal& x, const mpreal& y)
-inline const mpreal min(const mpreal& x, const mpreal& y)
 {
 	return (x<y?x:y);
 }
-#endif
+
+inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
+{
+	mpreal a;
+	mpfr_max(a.mp,x.mp,y.mp,rnd_mode);
+	return a;
+}
+
+inline const mpreal fmin(const mpreal& x, const mpreal& y,  mp_rnd_t rnd_mode)
+{
+	mpreal a;
+	mpfr_min(a.mp,x.mp,y.mp,rnd_mode);
+	return a;
+}
 
 inline const mpreal nexttoward (const mpreal& x, const mpreal& y)
 {
@@ -3105,6 +3197,7 @@ inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode)
 {
 	return pow(mpreal(a),static_cast<long int>(b),rnd_mode); // mpfr_pow_si
 }
+
 }
 
 // Explicit specialization of std::swap for mpreal numbers

unsupported/test/mpreal_support.cpp

@@ -41,4 +41,7 @@ void test_mpreal_support()
   }
 }
 
+extern "C" {
+#include "dlmalloc.c"
+}
 #include "mpreal.cpp"