GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-17 06:05:53 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
(transplanted from 79ad55a901434066e31523fcb47ce80cbc6897d5
)
This commit is contained in:
Gael Guennebaud 2011-08-19 15:03:45 +02:00
parent 5d5cf478ab
commit 4ab20b4cae
5 changed files with 9990 additions and 3530 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5703
unsupported/test/dlmalloc.c Executable file
View File

File diff suppressed because it is too large Load Diff

562
unsupported/test/dlmalloc.h Executable file
View File

@ -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 */

915
unsupported/test/mpreal.cpp
View File

@ -1,408 +1,507 @@
/* /*
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/ Project homepage: http://www.holoborodko.com/pavel/
Contact e-mail: pavel@holoborodko.com Contact e-mail: pavel@holoborodko.com
Copyright (c) 2008-2010 Pavel Holoborodko Copyright (c) 2008-2010 Pavel Holoborodko
This library is free software; you can redistribute it and/or Core Developers:
modify it under the terms of the GNU Lesser General Public Pavel Holoborodko, Dmitriy Gubanov, Konstantin Holoborodko.
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version. Contributors:
Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze,
This library is distributed in the hope that it will be useful, Heinz van Saanen, Pere Constans, Peter van Hoof.
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. This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
You should have received a copy of the GNU Lesser General Public License as published by the Free Software Foundation; either
License along with this library; if not, write to the Free Software version 2.1 of the License, or (at your option) any later version.
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
This library is distributed in the hope that it will be useful,
Contributors: but WITHOUT ANY WARRANTY; without even the implied warranty of
Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Heinz van Saanen, Pere Constans, Dmitriy Gubanov Lesser General Public License for more details.
*/
You should have received a copy of the GNU Lesser General Public
#include <cstring> License along with this library; if not, write to the Free Software
#include "mpreal.h" Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
using std::ws; ****************************************************************************
using std::cerr; Redistribution and use in source and binary forms, with or without
using std::endl; modification, are permitted provided that the following conditions
using std::string; are met:
using std::ostream;
using std::istream; 1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
namespace mpfr{
2. Redistributions in binary form must reproduce the above copyright
mp_rnd_t mpreal::default_rnd = mpfr_get_default_rounding_mode(); notice, this list of conditions and the following disclaimer in the
mp_prec_t mpreal::default_prec = mpfr_get_default_prec(); documentation and/or other materials provided with the distribution.
int mpreal::default_base = 10;
int mpreal::double_bits = -1; 3. Redistributions of any form whatsoever must retain the following
acknowledgment:
// Default constructor: creates mp number and initializes it to 0. "
mpreal::mpreal() This product includes software developed by Pavel Holoborodko
{ Web: http://www.holoborodko.com/pavel/
mpfr_init2(mp,default_prec); e-mail: pavel@holoborodko.com
mpfr_set_ui(mp,0,default_rnd); "
}
4. This software cannot be, by any means, used for any commercial
mpreal::mpreal(const mpreal& u) purpose without the prior permission of the copyright holder.
{
mpfr_init2(mp,mpfr_get_prec(u.mp)); Any of the above conditions can be waived if you get permission from
mpfr_set(mp,u.mp,default_rnd); the copyright holder.
}
THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
mpreal::mpreal(const mpfr_t u) ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
{ IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
mpfr_init2(mp,mpfr_get_prec(u)); ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
mpfr_set(mp,u,default_rnd); 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)
mpreal::mpreal(const mpf_t u) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
{ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
mpfr_init2(mp,mpf_get_prec(u)); OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
mpfr_set_f(mp,u,default_rnd); SUCH DAMAGE.
} */
#include <cstring>
mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode) #include <cstdlib>
{ #include "mpreal.h"
mpfr_init2(mp,prec); #include "dlmalloc.h"
mpfr_set_z(mp,u,mode);
} using std::ws;
using std::cerr;
mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode) using std::endl;
{ using std::string;
mpfr_init2(mp,prec); using std::ostream;
mpfr_set_q(mp,u,mode); using std::istream;
}
namespace mpfr{
mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
{ mp_rnd_t mpreal::default_rnd = mpfr_get_default_rounding_mode();
if(double_bits == -1 || fits_in_bits(u, double_bits)) mp_prec_t mpreal::default_prec = mpfr_get_default_prec();
{ int mpreal::default_base = 10;
mpfr_init2(mp,prec); int mpreal::double_bits = -1;
mpfr_set_d(mp,u,mode); bool mpreal::is_custom_malloc = false;
}
else // Default constructor: creates mp number and initializes it to 0.
throw conversion_overflow(); mpreal::mpreal()
} {
set_custom_malloc();
mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) mpfr_init2(mp,default_prec);
{ mpfr_set_ui(mp,0,default_rnd);
mpfr_init2(mp,prec); }
mpfr_set_ld(mp,u,mode);
} mpreal::mpreal(const mpreal& u)
{
mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) set_custom_malloc();
{ mpfr_init2(mp,mpfr_get_prec(u.mp));
mpfr_init2(mp,prec); mpfr_set(mp,u.mp,default_rnd);
mpfr_set_ui(mp,u,mode); }
}
mpreal::mpreal(const mpfr_t u)
mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) {
{ set_custom_malloc();
mpfr_init2(mp,prec); mpfr_init2(mp,mpfr_get_prec(u));
mpfr_set_ui(mp,u,mode); mpfr_set(mp,u,default_rnd);
} }
mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) mpreal::mpreal(const mpf_t u)
{ {
mpfr_init2(mp,prec); set_custom_malloc();
mpfr_set_si(mp,u,mode); mpfr_init2(mp,mpf_get_prec(u));
} mpfr_set_f(mp,u,default_rnd);
}
mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
{ mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
mpfr_init2(mp,prec); {
mpfr_set_si(mp,u,mode); set_custom_malloc();
} mpfr_init2(mp,prec);
mpfr_set_z(mp,u,mode);
mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) }
{
mpfr_init2(mp,prec); mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
mpfr_set_str(mp, s, base, mode); {
} set_custom_malloc();
mpfr_init2(mp,prec);
mpreal::~mpreal() mpfr_set_q(mp,u,mode);
{ }
mpfr_clear(mp);
} mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
{
// Operators - Assignment set_custom_malloc();
mpreal& mpreal::operator=(const char* s) if(double_bits == -1 || fits_in_bits(u, double_bits))
{ {
mpfr_t t; mpfr_init2(mp,prec);
if(0==mpfr_init_set_str(t,s,default_base,default_rnd)) mpfr_set_d(mp,u,mode);
{ }
mpfr_set(mp,t,mpreal::default_rnd); else
mpfr_clear(t); throw conversion_overflow();
}else{ }
mpfr_clear(t);
// cerr<<"fail to convert string"<<endl; mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
} {
set_custom_malloc();
return *this; mpfr_init2(mp,prec);
} mpfr_set_ld(mp,u,mode);
}
const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
{ mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
mpreal a; {
mp_prec_t p1, p2, p3; set_custom_malloc();
mpfr_init2(mp,prec);
p1 = v1.get_prec(); mpfr_set_ui(mp,u,mode);
p2 = v2.get_prec(); }
p3 = v3.get_prec();
mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); {
set_custom_malloc();
mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); mpfr_init2(mp,prec);
return a; mpfr_set_ui(mp,u,mode);
} }
const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
{ {
mpreal a; set_custom_malloc();
mp_prec_t p1, p2, p3; mpfr_init2(mp,prec);
mpfr_set_si(mp,u,mode);
p1 = v1.get_prec(); }
p2 = v2.get_prec();
p3 = v3.get_prec(); mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
{
a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); set_custom_malloc();
mpfr_init2(mp,prec);
mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); mpfr_set_si(mp,u,mode);
return a; }
}
mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode) {
{ set_custom_malloc();
mpreal a; mpfr_init2(mp,prec);
mp_prec_t p1, p2; mpfr_set_str(mp, s, base, mode);
}
p1 = v1.get_prec();
p2 = v2.get_prec(); mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)
{
a.set_prec(p1>p2?p1:p2); set_custom_malloc();
mpfr_init2(mp,prec);
mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode); mpfr_set_str(mp, s.c_str(), base, mode);
}
return a;
} mpreal::~mpreal()
{
const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) mpfr_clear(mp);
{ }
mpreal a;
mp_prec_t yp, xp; // Operators - Assignment
mpreal& mpreal::operator=(const char* s)
yp = y.get_prec(); {
xp = x.get_prec(); mpfr_t t;
a.set_prec(yp>xp?yp:xp); set_custom_malloc();
mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode); if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
{
return a; // 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);
const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode) mpfr_clear(t);
{ }else{
mpreal x; mpfr_clear(t);
mpfr_ptr* t; // cerr<<"fail to convert string"<<endl;
unsigned long int i; }
t = new mpfr_ptr[n]; return *this;
for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp; }
mpfr_sum(x.mp,t,n,rnd_mode);
delete[] t; const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
return x; {
} mpreal a;
mp_prec_t p1, p2, p3;
const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ p1 = v1.get_prec();
mpreal a; p2 = v2.get_prec();
mp_prec_t yp, xp; p3 = v3.get_prec();
yp = y.get_prec(); a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
xp = x.get_prec();
mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
a.set_prec(yp>xp?yp:xp); return a;
}
mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);
const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
return a; {
} mpreal a;
mp_prec_t p1, p2, p3;
const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{ p1 = v1.get_prec();
mpreal a; p2 = v2.get_prec();
mp_prec_t yp, xp; p3 = v3.get_prec();
yp = y.get_prec(); a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
xp = x.get_prec();
mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
a.set_prec(yp>xp?yp:xp); return a;
}
mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);
const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
return a; {
} mpreal a;
mp_prec_t p1, p2;
template <class T>
std::string to_string(T t, std::ios_base & (*f)(std::ios_base&)) p1 = v1.get_prec();
{ p2 = v2.get_prec();
std::ostringstream oss;
oss << f << t; a.set_prec(p1>p2?p1:p2);
return oss.str();
} mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);
mpreal::operator std::string() const return a;
{ }
return to_string();
} const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{
string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const mpreal a;
{ mp_prec_t yp, xp;
char *s, *ns = NULL;
size_t slen, nslen; yp = y.get_prec();
mp_exp_t exp; xp = x.get_prec();
string out;
a.set_prec(yp>xp?yp:xp);
if(mpfr_inf_p(mp))
{ mpfr_hypot(a.mp, x.mp, y.mp, rnd_mode);
if(mpfr_sgn(mp)>0) return "+@Inf@";
else return "-@Inf@"; return a;
} }
if(mpfr_zero_p(mp)) return "0"; const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
if(mpfr_nan_p(mp)) return "@NaN@"; {
mpreal x;
mpfr_ptr* t;
s = mpfr_get_str(NULL,&exp,b,0,mp,mode); unsigned long int i;
ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
t = new mpfr_ptr[n];
if(s!=NULL && ns!=NULL) for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp;
{ mpfr_sum(x.mp,t,n,rnd_mode);
slen = strlen(s); delete[] t;
nslen = strlen(ns); return x;
if(nslen<=slen) }
{
mpfr_free_str(s); const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
s = ns; {
slen = nslen; mpreal a;
} mp_prec_t yp, xp;
else {
mpfr_free_str(ns); yp = y.get_prec();
} xp = x.get_prec();
// Make human eye-friendly formatting if possible a.set_prec(yp>xp?yp:xp);
if (exp>0 && static_cast<size_t>(exp)<slen)
{ mpfr_remainder(a.mp, x.mp, y.mp, rnd_mode);
if(s[0]=='-')
{ return a;
// Remove zeros starting from right end }
char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s+exp) ptr--; const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
{
if(ptr==s+exp) out = string(s,exp+1); mpreal a;
else out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1); mp_prec_t yp, xp;
//out = string(s,exp+1)+'.'+string(s+exp+1); yp = y.get_prec();
} xp = x.get_prec();
else
{ a.set_prec(yp>xp?yp:xp);
// Remove zeros starting from right end
char* ptr = s+slen-1; mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);
while (*ptr=='0' && ptr>s+exp-1) ptr--;
return a;
if(ptr==s+exp-1) out = string(s,exp); }
else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1);
template <class T>
//out = string(s,exp)+'.'+string(s+exp); std::string to_string(T t, std::ios_base & (*f)(std::ios_base&))
} {
std::ostringstream oss;
}else{ // exp<0 || exp>slen oss << f << t;
if(s[0]=='-') return oss.str();
{ }
// Remove zeros starting from right end
char* ptr = s+slen-1; mpreal::operator std::string() const
while (*ptr=='0' && ptr>s+1) ptr--; {
return to_string();
if(ptr==s+1) out = string(s,2); }
else out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1);
std::string mpreal::to_string(size_t n, int b, mp_rnd_t mode) const
//out = string(s,2)+'.'+string(s+2); {
} char *s, *ns = NULL;
else size_t slen, nslen;
{ mp_exp_t exp;
// Remove zeros starting from right end string out;
char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s) ptr--; set_custom_malloc();
if(ptr==s) out = string(s,1); if(mpfr_inf_p(mp))
else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1); {
if(mpfr_sgn(mp)>0) return "+@Inf@";
//out = string(s,1)+'.'+string(s+1); else return "-@Inf@";
} }
// Make final string if(mpfr_zero_p(mp)) return "0";
if(--exp) if(mpfr_nan_p(mp)) return "@NaN@";
{
if(exp>0) out += "e+"+mpfr::to_string<mp_exp_t>(exp,std::dec); s = mpfr_get_str(NULL,&exp,b,0,mp,mode);
else out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec); ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
}
} if(s!=NULL && ns!=NULL)
{
mpfr_free_str(s); slen = strlen(s);
return out; nslen = strlen(ns);
}else{ if(nslen<=slen)
return "conversion error!"; {
} mpfr_free_str(s);
} s = ns;
slen = nslen;
////////////////////////////////////////////////////////////////////////// }
// I/O else {
ostream& operator<<(ostream& os, const mpreal& v) mpfr_free_str(ns);
{ }
return os<<v.to_string(os.precision());
} // Make human eye-friendly formatting if possible
if (exp>0 && static_cast<size_t>(exp)<slen)
istream& operator>>(istream &is, mpreal& v) {
{ if(s[0]=='-')
char c; {
string s = ""; // Remove zeros starting from right end
mpfr_t t; char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s+exp) ptr--;
if(is.good())
{ if(ptr==s+exp) out = string(s,exp+1);
is>>ws; else out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1);
while ((c = is.get())!=EOF)
{ //out = string(s,exp+1)+'.'+string(s+exp+1);
if(c ==' ' || c == '\t' || c == '\n' || c == '\r') }
{ else
is.putback(c); {
break; // Remove zeros starting from right end
} char* ptr = s+slen-1;
s += c; while (*ptr=='0' && ptr>s+exp-1) ptr--;
}
if(ptr==s+exp-1) out = string(s,exp);
if(s.size() != 0) else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1);
{
// Protect current value from alternation in case of input error //out = string(s,exp)+'.'+string(s+exp);
// 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))
{ }else{ // exp<0 || exp>slen
mpfr_set(v.mp,t,mpreal::default_rnd); if(s[0]=='-')
mpfr_clear(t); {
// Remove zeros starting from right end
}else{ char* ptr = s+slen-1;
mpfr_clear(t); while (*ptr=='0' && ptr>s+1) ptr--;
cerr<<"error reading from istream"<<endl;
// throw an exception 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);
return is; }
} else
} {
// Remove zeros starting from right end
char* ptr = s+slen-1;
while (*ptr=='0' && ptr>s) ptr--;
if(ptr==s) out = string(s,1);
else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1);
//out = string(s,1)+'.'+string(s+1);
}
// Make final string
if(--exp)
{
if(exp>0) out += "e+"+mpfr::to_string<mp_exp_t>(exp,std::dec);
else out += "e"+mpfr::to_string<mp_exp_t>(exp,std::dec);
}
}
mpfr_free_str(s);
return out;
}else{
return "conversion error!";
}
}
//////////////////////////////////////////////////////////////////////////
// I/O
ostream& operator<<(ostream& os, const mpreal& v)
{
return os<<v.to_string(static_cast<size_t>(os.precision()));
}
istream& operator>>(istream &is, mpreal& v)
{
char c;
string s = "";
mpfr_t t;
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;
}
// 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;
}
}
}

6337
unsupported/test/mpreal.h
View File

File diff suppressed because it is too large Load Diff

3
unsupported/test/mpreal_support.cpp
View File

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