Ticket #6296: mpreal.cpp

File mpreal.cpp, 12.9 KB (added by anonymous, 11 years ago)
Line 
1/*
2 Multi-precision real number class. C++ interface fo MPFR library.
3 Project homepage: http://www.holoborodko.com/pavel/
4 Contact e-mail: pavel@holoborodko.com
5
6 Copyright (c) 2008-2011 Pavel Holoborodko
7
8 Core Developers:
9 Pavel Holoborodko, Dmitriy Gubanov, Konstantin Holoborodko.
10
11 Contributors:
12 Brian Gladman, Helmut Jarausch, Fokko Beekhof, Ulrich Mutze,
13 Heinz van Saanen, Pere Constans, Peter van Hoof, Gael Guennebaud,
14 Tsai Chia Cheng, Alexei Zubanov.
15
16 ****************************************************************************
17 This library is free software; you can redistribute it and/or
18 modify it under the terms of the GNU Lesser General Public
19 License as published by the Free Software Foundation; either
20 version 2.1 of the License, or (at your option) any later version.
21
22 This library is distributed in the hope that it will be useful,
23 but WITHOUT ANY WARRANTY; without even the implied warranty of
24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25 Lesser General Public License for more details.
26
27 You should have received a copy of the GNU Lesser General Public
28 License along with this library; if not, write to the Free Software
29 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
30
31 ****************************************************************************
32 ****************************************************************************
33 Redistribution and use in source and binary forms, with or without
34 modification, are permitted provided that the following conditions
35 are met:
36
37 1. Redistributions of source code must retain the above copyright
38 notice, this list of conditions and the following disclaimer.
39
40 2. Redistributions in binary form must reproduce the above copyright
41 notice, this list of conditions and the following disclaimer in the
42 documentation and/or other materials provided with the distribution.
43
44 3. The name of the author may be used to endorse or promote products
45 derived from this software without specific prior written permission.
46
47 THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
48 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
49 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
50 ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
51 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
52 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
53 OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
54 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
55 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
56 OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
57 SUCH DAMAGE.
58*/
59#include <cstring>
60#include "mpreal.h"
61//#include "dlmalloc.h"
62#include <stdlib.h>
63
64using std::ws;
65using std::cerr;
66using std::endl;
67using std::string;
68using std::ostream;
69using std::istream;
70
71namespace mpfr{
72
73mp_rnd_t mpreal::default_rnd = MPFR_RNDN; //(mpfr_get_default_rounding_mode)();
74mp_prec_t mpreal::default_prec = 64; //(mpfr_get_default_prec)();
75int mpreal::default_base = 10;
76int mpreal::double_bits = -1;
77bool mpreal::is_custom_malloc = false;
78
79// Default constructor: creates mp number and initializes it to 0.
80mpreal::mpreal()
81{
82 set_custom_malloc();
83 mpfr_init2(mp,default_prec);
84 mpfr_set_ui(mp,0,default_rnd);
85
86#if defined(_MSC_VER) && defined(_DEBUG)
87 setDebugView();
88#endif
89
90}
91
92mpreal::mpreal(const mpreal& u)
93{
94 set_custom_malloc();
95 mpfr_init2(mp,mpfr_get_prec(u.mp));
96 mpfr_set(mp,u.mp,default_rnd);
97
98#if defined(_MSC_VER) && defined(_DEBUG)
99 setDebugView();
100#endif
101
102}
103
104mpreal::mpreal(const mpfr_t u)
105{
106 set_custom_malloc();
107 mpfr_init2(mp,mpfr_get_prec(u));
108 mpfr_set(mp,u,default_rnd);
109
110#if defined(_MSC_VER) && defined(_DEBUG)
111 setDebugView();
112#endif
113
114}
115
116mpreal::mpreal(const mpf_t u)
117{
118 set_custom_malloc();
119 mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t)
120 mpfr_set_f(mp,u,default_rnd);
121
122#if defined(_MSC_VER) && defined(_DEBUG)
123 setDebugView();
124#endif
125
126}
127
128mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode)
129{
130 set_custom_malloc();
131 mpfr_init2(mp,prec);
132 mpfr_set_z(mp,u,mode);
133
134#if defined(_MSC_VER) && defined(_DEBUG)
135 setDebugView();
136#endif
137
138}
139
140mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode)
141{
142 set_custom_malloc();
143 mpfr_init2(mp,prec);
144 mpfr_set_q(mp,u,mode);
145
146#if defined(_MSC_VER) && defined(_DEBUG)
147 setDebugView();
148#endif
149
150}
151
152mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
153{
154 set_custom_malloc();
155 if(double_bits == -1 || fits_in_bits(u, double_bits))
156 {
157 mpfr_init2(mp,prec);
158 mpfr_set_d(mp,u,mode);
159
160#if defined(_MSC_VER) && defined(_DEBUG)
161 setDebugView();
162#endif
163
164 }
165 else
166 throw conversion_overflow();
167}
168
169mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode)
170{
171 set_custom_malloc();
172 mpfr_init2(mp,prec);
173 mpfr_set_ld(mp,u,mode);
174
175#if defined(_MSC_VER) && defined(_DEBUG)
176 setDebugView();
177#endif
178
179}
180
181mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode)
182{
183 set_custom_malloc();
184 mpfr_init2(mp,prec);
185 mpfr_set_ui(mp,u,mode);
186
187#if defined(_MSC_VER) && defined(_DEBUG)
188 setDebugView();
189#endif
190
191}
192
193mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode)
194{
195 set_custom_malloc();
196 mpfr_init2(mp,prec);
197 mpfr_set_ui(mp,u,mode);
198
199#if defined(_MSC_VER) && defined(_DEBUG)
200 setDebugView();
201#endif
202
203}
204
205mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode)
206{
207 set_custom_malloc();
208 mpfr_init2(mp,prec);
209 mpfr_set_si(mp,u,mode);
210
211#if defined(_MSC_VER) && defined(_DEBUG)
212 setDebugView();
213#endif
214
215}
216
217mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode)
218{
219 set_custom_malloc();
220 mpfr_init2(mp,prec);
221 mpfr_set_si(mp,u,mode);
222
223#if defined(_MSC_VER) && defined(_DEBUG)
224 setDebugView();
225#endif
226
227}
228
229#if defined (MPREAL_HAVE_INT64_SUPPORT)
230mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode)
231{
232 set_custom_malloc();
233 mpfr_init2(mp,prec);
234 mpfr_set_uj(mp, u, mode);
235
236#if defined(_MSC_VER) && defined(_DEBUG)
237 setDebugView();
238#endif
239
240}
241
242mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode)
243{
244 set_custom_malloc();
245 mpfr_init2(mp,prec);
246 mpfr_set_sj(mp, u, mode);
247
248#if defined(_MSC_VER) && defined(_DEBUG)
249 setDebugView();
250#endif
251
252}
253#endif
254
255mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode)
256{
257 set_custom_malloc();
258 mpfr_init2(mp,prec);
259 mpfr_set_str(mp, s, base, mode);
260
261#if defined(_MSC_VER) && defined(_DEBUG)
262 setDebugView();
263#endif
264
265}
266
267mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode)
268{
269 set_custom_malloc();
270 mpfr_init2(mp,prec);
271 mpfr_set_str(mp, s.c_str(), base, mode);
272
273#if defined(_MSC_VER) && defined(_DEBUG)
274 setDebugView();
275#endif
276
277}
278
279mpreal::~mpreal()
280{
281 mpfr_clear(mp);
282}
283
284// Operators - Assignment
285mpreal& mpreal::operator=(const char* s)
286{
287 mpfr_t t;
288
289 set_custom_malloc();
290
291 if(0==mpfr_init_set_str(t,s,default_base,default_rnd))
292 {
293 // We will rewrite mp anyway, so flash it and resize
294 mpfr_set_prec(mp,mpfr_get_prec(t));
295 mpfr_set(mp,t,mpreal::default_rnd);
296 mpfr_clear(t);
297
298#if defined(_MSC_VER) && defined(_DEBUG)
299 setDebugView();
300#endif
301
302 }else{
303 mpfr_clear(t);
304 }
305
306 return *this;
307}
308
309const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
310{
311 mpreal a;
312 mp_prec_t p1, p2, p3;
313
314 p1 = v1.get_prec();
315 p2 = v2.get_prec();
316 p3 = v3.get_prec();
317
318 a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
319
320 mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
321 return a;
322}
323
324const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode)
325{
326 mpreal a;
327 mp_prec_t p1, p2, p3;
328
329 p1 = v1.get_prec();
330 p2 = v2.get_prec();
331 p3 = v3.get_prec();
332
333 a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1));
334
335 mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode);
336 return a;
337}
338
339const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode)
340{
341 mpreal a;
342 mp_prec_t p1, p2;
343
344 p1 = v1.get_prec();
345 p2 = v2.get_prec();
346
347 a.set_prec(p1>p2?p1:p2);
348
349 mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode);
350
351 return a;
352}
353
354const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode)
355{
356 mpreal x;
357 mpfr_ptr* t;
358 unsigned long int i;
359
360 t = new mpfr_ptr[n];
361 for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp;
362 mpfr_sum(x.mp,t,n,rnd_mode);
363 delete[] t;
364 return x;
365}
366
367const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode)
368{
369 mpreal a;
370 mp_prec_t yp, xp;
371
372 yp = y.get_prec();
373 xp = x.get_prec();
374
375 a.set_prec(yp>xp?yp:xp);
376
377 mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode);
378
379 return a;
380}
381
382template <class T>
383std::string toString(T t, std::ios_base & (*f)(std::ios_base&))
384{
385 std::ostringstream oss;
386 oss << f << t;
387 return oss.str();
388}
389
390mpreal::operator std::string() const
391{
392 return toString();
393}
394
395#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
396
397std::string mpreal::toString(const std::string& format) const
398{
399 char *s;
400 string out;
401
402 if( !format.empty() )
403 {
404 if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0))
405 {
406 out = std::string(s);
407
408 mpfr_free_str(s);
409 }
410 }
411
412 return out;
413}
414
415#endif
416
417std::string mpreal::toString(size_t n, int b, mp_rnd_t mode) const
418{
419
420#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
421
422 // Use MPFR native function for output
423 char format[128];
424 int digits;
425
426 digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp));
427
428 sprintf(format,"%%.%dRNg",digits); // Default format
429
430 return toString(std::string(format));
431
432#else
433
434 char *s, *ns = NULL;
435 size_t slen, nslen;
436 mp_exp_t exp;
437 string out;
438
439 set_custom_malloc();
440
441 if(mpfr_inf_p(mp))
442 {
443 if(mpfr_sgn(mp)>0) return "+Inf";
444 else return "-Inf";
445 }
446
447 if(mpfr_zero_p(mp)) return "0";
448 if(mpfr_nan_p(mp)) return "NaN";
449
450 s = mpfr_get_str(NULL,&exp,b,0,mp,mode);
451 ns = mpfr_get_str(NULL,&exp,b,n,mp,mode);
452
453 if(s!=NULL && ns!=NULL)
454 {
455 slen = strlen(s);
456 nslen = strlen(ns);
457 if(nslen<=slen)
458 {
459 mpfr_free_str(s);
460 s = ns;
461 slen = nslen;
462 }
463 else {
464 mpfr_free_str(ns);
465 }
466
467 // Make human eye-friendly formatting if possible
468 if (exp>0 && static_cast<size_t>(exp)<slen)
469 {
470 if(s[0]=='-')
471 {
472 // Remove zeros starting from right end
473 char* ptr = s+slen-1;
474 while (*ptr=='0' && ptr>s+exp) ptr--;
475
476 if(ptr==s+exp) out = string(s,exp+1);
477 else out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1);
478
479 //out = string(s,exp+1)+'.'+string(s+exp+1);
480 }
481 else
482 {
483 // Remove zeros starting from right end
484 char* ptr = s+slen-1;
485 while (*ptr=='0' && ptr>s+exp-1) ptr--;
486
487 if(ptr==s+exp-1) out = string(s,exp);
488 else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1);
489
490 //out = string(s,exp)+'.'+string(s+exp);
491 }
492
493 }else{ // exp<0 || exp>slen
494 if(s[0]=='-')
495 {
496 // Remove zeros starting from right end
497 char* ptr = s+slen-1;
498 while (*ptr=='0' && ptr>s+1) ptr--;
499
500 if(ptr==s+1) out = string(s,2);
501 else out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1);
502
503 //out = string(s,2)+'.'+string(s+2);
504 }
505 else
506 {
507 // Remove zeros starting from right end
508 char* ptr = s+slen-1;
509 while (*ptr=='0' && ptr>s) ptr--;
510
511 if(ptr==s) out = string(s,1);
512 else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1);
513
514 //out = string(s,1)+'.'+string(s+1);
515 }
516
517 // Make final string
518 if(--exp)
519 {
520 if(exp>0) out += "e+"+mpfr::toString<mp_exp_t>(exp,std::dec);
521 else out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec);
522 }
523 }
524
525 mpfr_free_str(s);
526 return out;
527 }else{
528 return "conversion error!";
529 }
530#endif
531}
532
533
534//////////////////////////////////////////////////////////////////////////
535// I/O
536ostream& operator<<(ostream& os, const mpreal& v)
537{
538 return os<<v.toString(static_cast<size_t>(os.precision()));
539}
540
541istream& operator>>(istream &is, mpreal& v)
542{
543 string tmp;
544 is >> tmp;
545 mpfr_set_str(v.mp, tmp.c_str(),mpreal::default_base,mpreal::default_rnd);
546 return is;
547}
548
549// Optimized dynamic memory allocation/(re-)deallocation.
550void * mpreal::mpreal_allocate(size_t alloc_size)
551{
552 //return(dlmalloc(alloc_size));
553 return(malloc(alloc_size));
554}
555
556void * mpreal::mpreal_reallocate(void *ptr, size_t old_size, size_t new_size)
557{
558 //return(dlrealloc(ptr,new_size));
559 return(realloc(ptr,new_size));
560}
561
562void mpreal::mpreal_free(void *ptr, size_t size)
563{
564 //dlfree(ptr);
565 free(ptr); ptr = NULL;
566}
567
568inline void mpreal::set_custom_malloc(void)
569{
570 if(!is_custom_malloc)
571 {
572 mp_set_memory_functions(mpreal_allocate,mpreal_reallocate,mpreal_free);
573 is_custom_malloc = true;
574 }
575}
576}
577