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 62#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 63#include "dlmalloc.h" 64#endif 65 66using std::ws; 67using std::cerr; 68using std::endl; 69using std::string; 70using std::ostream; 71using std::istream; 72 73namespace mpfr{ 74 75mp_rnd_t mpreal::default_rnd = MPFR_RNDN; //(mpfr_get_default_rounding_mode)(); 76mp_prec_t mpreal::default_prec = 64; //(mpfr_get_default_prec)(); 77int mpreal::default_base = 10; 78int mpreal::double_bits = -1; 79 80#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 81bool mpreal::is_custom_malloc = false; 82#endif 83 84// Default constructor: creates mp number and initializes it to 0. 85mpreal::mpreal() 86{ 87 88#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 89 set_custom_malloc(); 90#endif 91 92 mpfr_init2(mp,default_prec); 93 mpfr_set_ui(mp,0,default_rnd); 94 95 MPREAL_MSVC_DEBUGVIEW_CODE; 96} 97 98mpreal::mpreal(const mpreal& u) 99{ 100 101#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 102 set_custom_malloc(); 103#endif 104 105 mpfr_init2(mp,mpfr_get_prec(u.mp)); 106 mpfr_set(mp,u.mp,default_rnd); 107 108 MPREAL_MSVC_DEBUGVIEW_CODE; 109} 110 111mpreal::mpreal(const mpfr_t u) 112{ 113 114#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 115 set_custom_malloc(); 116#endif 117 118 mpfr_init2(mp,mpfr_get_prec(u)); 119 mpfr_set(mp,u,default_rnd); 120 121 MPREAL_MSVC_DEBUGVIEW_CODE; 122} 123 124mpreal::mpreal(const mpf_t u) 125{ 126 127#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 128 set_custom_malloc(); 129#endif 130 131 mpfr_init2(mp,(mp_prec_t) mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t) 132 mpfr_set_f(mp,u,default_rnd); 133 134 MPREAL_MSVC_DEBUGVIEW_CODE; 135} 136 137mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode) 138{ 139 140#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 141 set_custom_malloc(); 142#endif 143 144 mpfr_init2(mp,prec); 145 mpfr_set_z(mp,u,mode); 146 147 MPREAL_MSVC_DEBUGVIEW_CODE; 148} 149 150mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode) 151{ 152 153#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 154 set_custom_malloc(); 155#endif 156 157 mpfr_init2(mp,prec); 158 mpfr_set_q(mp,u,mode); 159 160 MPREAL_MSVC_DEBUGVIEW_CODE; 161} 162 163mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) 164{ 165 166#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 167 set_custom_malloc(); 168#endif 169 170 if(double_bits == -1 || fits_in_bits(u, double_bits)) 171 { 172 mpfr_init2(mp,prec); 173 mpfr_set_d(mp,u,mode); 174 175 MPREAL_MSVC_DEBUGVIEW_CODE; 176 } 177 else 178 throw conversion_overflow(); 179} 180 181mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) 182{ 183 184#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 185 set_custom_malloc(); 186#endif 187 188 mpfr_init2(mp,prec); 189 mpfr_set_ld(mp,u,mode); 190 191 MPREAL_MSVC_DEBUGVIEW_CODE; 192} 193 194mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) 195{ 196 197#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 198 set_custom_malloc(); 199#endif 200 201 mpfr_init2(mp,prec); 202 mpfr_set_ui(mp,u,mode); 203 204 MPREAL_MSVC_DEBUGVIEW_CODE; 205} 206 207mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) 208{ 209 210#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 211 set_custom_malloc(); 212#endif 213 214 mpfr_init2(mp,prec); 215 mpfr_set_ui(mp,u,mode); 216 217 MPREAL_MSVC_DEBUGVIEW_CODE; 218} 219 220mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) 221{ 222 223#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 224 set_custom_malloc(); 225#endif 226 227 mpfr_init2(mp,prec); 228 mpfr_set_si(mp,u,mode); 229 230 MPREAL_MSVC_DEBUGVIEW_CODE; 231} 232 233mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) 234{ 235 236#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 237 set_custom_malloc(); 238#endif 239 240 mpfr_init2(mp,prec); 241 mpfr_set_si(mp,u,mode); 242 243 MPREAL_MSVC_DEBUGVIEW_CODE; 244} 245 246#if defined (MPREAL_HAVE_INT64_SUPPORT) 247mpreal::mpreal(const uint64_t u, mp_prec_t prec, mp_rnd_t mode) 248{ 249 250#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 251 set_custom_malloc(); 252#endif 253 254 mpfr_init2(mp,prec); 255 mpfr_set_uj(mp, u, mode); 256 257 MPREAL_MSVC_DEBUGVIEW_CODE; 258} 259 260mpreal::mpreal(const int64_t u, mp_prec_t prec, mp_rnd_t mode) 261{ 262 263#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 264 set_custom_malloc(); 265#endif 266 267 mpfr_init2(mp,prec); 268 mpfr_set_sj(mp, u, mode); 269 270 MPREAL_MSVC_DEBUGVIEW_CODE; 271} 272#endif 273 274mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) 275{ 276 277#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 278 set_custom_malloc(); 279#endif 280 281 mpfr_init2(mp,prec); 282 mpfr_set_str(mp, s, base, mode); 283 284 MPREAL_MSVC_DEBUGVIEW_CODE; 285} 286 287mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode) 288{ 289 290#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 291 set_custom_malloc(); 292#endif 293 294 mpfr_init2(mp,prec); 295 mpfr_set_str(mp, s.c_str(), base, mode); 296 297 MPREAL_MSVC_DEBUGVIEW_CODE; 298} 299 300mpreal::~mpreal() 301{ 302 mpfr_clear(mp); 303} 304 305// Operators - Assignment 306mpreal& mpreal::operator=(const char* s) 307{ 308 mpfr_t t; 309 310#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 311 set_custom_malloc(); 312#endif 313 314 if(0==mpfr_init_set_str(t,s,default_base,default_rnd)) 315 { 316 // We will rewrite mp anyway, so flash it and resize 317 mpfr_set_prec(mp,mpfr_get_prec(t)); 318 mpfr_set(mp,t,mpreal::default_rnd); 319 mpfr_clear(t); 320 321 MPREAL_MSVC_DEBUGVIEW_CODE; 322 323 }else{ 324 mpfr_clear(t); 325 } 326 327 return *this; 328} 329 330const mpreal fma (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) 331{ 332 mpreal a; 333 mp_prec_t p1, p2, p3; 334 335 p1 = v1.get_prec(); 336 p2 = v2.get_prec(); 337 p3 = v3.get_prec(); 338 339 a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); 340 341 mpfr_fma(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); 342 return a; 343} 344 345const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode) 346{ 347 mpreal a; 348 mp_prec_t p1, p2, p3; 349 350 p1 = v1.get_prec(); 351 p2 = v2.get_prec(); 352 p3 = v3.get_prec(); 353 354 a.set_prec(p3>p2?(p3>p1?p3:p1):(p2>p1?p2:p1)); 355 356 mpfr_fms(a.mp,v1.mp,v2.mp,v3.mp,rnd_mode); 357 return a; 358} 359 360const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode) 361{ 362 mpreal a; 363 mp_prec_t p1, p2; 364 365 p1 = v1.get_prec(); 366 p2 = v2.get_prec(); 367 368 a.set_prec(p1>p2?p1:p2); 369 370 mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode); 371 372 return a; 373} 374 375const mpreal sum (const mpreal tab[], unsigned long int n, mp_rnd_t rnd_mode) 376{ 377 mpreal x; 378 mpfr_ptr* t; 379 unsigned long int i; 380 381 t = new mpfr_ptr[n]; 382 for (i=0;i<n;i++) t[i] = (mpfr_ptr)tab[i].mp; 383 mpfr_sum(x.mp,t,n,rnd_mode); 384 delete[] t; 385 return x; 386} 387 388const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode) 389{ 390 mpreal a; 391 mp_prec_t yp, xp; 392 393 yp = y.get_prec(); 394 xp = x.get_prec(); 395 396 a.set_prec(yp>xp?yp:xp); 397 398 mpfr_remquo(a.mp,q, x.mp, y.mp, rnd_mode); 399 400 return a; 401} 402 403template <class T> 404std::string toString(T t, std::ios_base & (*f)(std::ios_base&)) 405{ 406 std::ostringstream oss; 407 oss << f << t; 408 return oss.str(); 409} 410 411#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 412 413std::string mpreal::toString(const std::string& format) const 414{ 415 char *s = NULL; 416 string out; 417 418 if( !format.empty() ) 419 { 420 if(!(mpfr_asprintf(&s,format.c_str(),mp) < 0)) 421 { 422 out = std::string(s); 423 424 mpfr_free_str(s); 425 } 426 } 427 428 return out; 429} 430 431#endif 432 433std::string mpreal::toString(int n, int b, mp_rnd_t mode) const 434{ 435 (void)b; 436 (void)mode; 437#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0)) 438 439 // Use MPFR native function for output 440 char format[128]; 441 int digits; 442 443 digits = n > 0 ? n : bits2digits(mpfr_get_prec(mp)); 444 445 sprintf(format,"%%.%dRNg",digits); // Default format 446 447 return toString(std::string(format)); 448 449#else 450 451 char *s, *ns = NULL; 452 size_t slen, nslen; 453 mp_exp_t exp; 454 string out; 455 456#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 457 set_custom_malloc(); 458#endif 459 460 if(mpfr_inf_p(mp)) 461 { 462 if(mpfr_sgn(mp)>0) return "+Inf"; 463 else return "-Inf"; 464 } 465 466 if(mpfr_zero_p(mp)) return "0"; 467 if(mpfr_nan_p(mp)) return "NaN"; 468 469 s = mpfr_get_str(NULL,&exp,b,0,mp,mode); 470 ns = mpfr_get_str(NULL,&exp,b,n,mp,mode); 471 472 if(s!=NULL && ns!=NULL) 473 { 474 slen = strlen(s); 475 nslen = strlen(ns); 476 if(nslen<=slen) 477 { 478 mpfr_free_str(s); 479 s = ns; 480 slen = nslen; 481 } 482 else { 483 mpfr_free_str(ns); 484 } 485 486 // Make human eye-friendly formatting if possible 487 if (exp>0 && static_cast<size_t>(exp)<slen) 488 { 489 if(s[0]=='-') 490 { 491 // Remove zeros starting from right end 492 char* ptr = s+slen-1; 493 while (*ptr=='0' && ptr>s+exp) ptr--; 494 495 if(ptr==s+exp) out = string(s,exp+1); 496 else out = string(s,exp+1)+'.'+string(s+exp+1,ptr-(s+exp+1)+1); 497 498 //out = string(s,exp+1)+'.'+string(s+exp+1); 499 } 500 else 501 { 502 // Remove zeros starting from right end 503 char* ptr = s+slen-1; 504 while (*ptr=='0' && ptr>s+exp-1) ptr--; 505 506 if(ptr==s+exp-1) out = string(s,exp); 507 else out = string(s,exp)+'.'+string(s+exp,ptr-(s+exp)+1); 508 509 //out = string(s,exp)+'.'+string(s+exp); 510 } 511 512 }else{ // exp<0 || exp>slen 513 if(s[0]=='-') 514 { 515 // Remove zeros starting from right end 516 char* ptr = s+slen-1; 517 while (*ptr=='0' && ptr>s+1) ptr--; 518 519 if(ptr==s+1) out = string(s,2); 520 else out = string(s,2)+'.'+string(s+2,ptr-(s+2)+1); 521 522 //out = string(s,2)+'.'+string(s+2); 523 } 524 else 525 { 526 // Remove zeros starting from right end 527 char* ptr = s+slen-1; 528 while (*ptr=='0' && ptr>s) ptr--; 529 530 if(ptr==s) out = string(s,1); 531 else out = string(s,1)+'.'+string(s+1,ptr-(s+1)+1); 532 533 //out = string(s,1)+'.'+string(s+1); 534 } 535 536 // Make final string 537 if(--exp) 538 { 539 if(exp>0) out += "e+"+mpfr::toString<mp_exp_t>(exp,std::dec); 540 else out += "e"+mpfr::toString<mp_exp_t>(exp,std::dec); 541 } 542 } 543 544 mpfr_free_str(s); 545 return out; 546 }else{ 547 return "conversion error!"; 548 } 549#endif 550} 551 552 553////////////////////////////////////////////////////////////////////////// 554// I/O 555ostream& operator<<(ostream& os, const mpreal& v) 556{ 557 return os<<v.toString(static_cast<int>(os.precision())); 558} 559 560istream& operator>>(istream &is, mpreal& v) 561{ 562 string tmp; 563 is >> tmp; 564 mpfr_set_str(v.mp, tmp.c_str(),mpreal::default_base,mpreal::default_rnd); 565 return is; 566} 567 568 569#if defined (MPREAL_HAVE_CUSTOM_MPFR_MALLOC) 570 // Optimized dynamic memory allocation/(re-)deallocation. 571 void * mpreal::mpreal_allocate(size_t alloc_size) 572 { 573 return(dlmalloc(alloc_size)); 574 } 575 576 void * mpreal::mpreal_reallocate(void *ptr, size_t old_size, size_t new_size) 577 { 578 return(dlrealloc(ptr,new_size)); 579 } 580 581 void mpreal::mpreal_free(void *ptr, size_t size) 582 { 583 dlfree(ptr); 584 } 585 586 inline void mpreal::set_custom_malloc(void) 587 { 588 if(!is_custom_malloc) 589 { 590 mp_set_memory_functions(mpreal_allocate,mpreal_reallocate,mpreal_free); 591 is_custom_malloc = true; 592 } 593 } 594#endif 595 596} 597 598