ginac/numeric.cpp

   1 /** @file numeric.cpp
   2  *
   3  *  This file contains the interface to the underlying bignum package.
   4  *  Its most important design principle is to completely hide the inner
   5  *  working of that other package from the user of GiNaC.  It must either
   6  *  provide implementation of arithmetic operators and numerical evaluation
   7  *  of special functions or implement the interface to the bignum package.
   8  *
   9  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
  10  *
  11  *  This program is free software; you can redistribute it and/or modify
  12  *  it under the terms of the GNU General Public License as published by
  13  *  the Free Software Foundation; either version 2 of the License, or
  14  *  (at your option) any later version.
  15  *
  16  *  This program is distributed in the hope that it will be useful,
  17  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19  *  GNU General Public License for more details.
  20  *
  21  *  You should have received a copy of the GNU General Public License
  22  *  along with this program; if not, write to the Free Software
  23  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  24  */
  25
  26 #include <vector>
  27 #include <stdexcept>
  28
  29 #include "ginac.h"
  30 #include "config.h"
  31
  32 // CLN should not pollute the global namespace, hence we include it here
  33 // instead of in some header file where it would propagate to other parts:
  34 #ifdef HAVE_CLN_CLN_H
  35 #include <CLN/cln.h>
  36 #else
  37 #include <cln.h>
  38 #endif
  39
  40 // linker has no problems finding text symbols for numerator or denominator
  41 //#define SANE_LINKER
  42
  43 //////////
  44 // default constructor, destructor, copy constructor assignment
  45 // operator and helpers
  46 //////////
  47
  48 // public
  49
  50 /** default ctor. Numerically it initializes to an integer zero. */
  51 numeric::numeric() : basic(TINFO_NUMERIC)
  52 {
  53     debugmsg("numeric default constructor", LOGLEVEL_CONSTRUCT);
  54     value = new cl_N;
  55     *value=cl_I(0);
  56     calchash();
  57     setflag(status_flags::evaluated|
  58             status_flags::hash_calculated);
  59 }
  60
  61 numeric::~numeric()
  62 {
  63     debugmsg("numeric destructor" ,LOGLEVEL_DESTRUCT);
  64     destroy(0);
  65 }
  66
  67 numeric::numeric(numeric const & other)
  68 {
  69     debugmsg("numeric copy constructor", LOGLEVEL_CONSTRUCT);
  70     copy(other);
  71 }
  72
  73 numeric const & numeric::operator=(numeric const & other)
  74 {
  75     debugmsg("numeric operator=", LOGLEVEL_ASSIGNMENT);
  76     if (this != &other) {
  77         destroy(1);
  78         copy(other);
  79     }
  80     return *this;
  81 }
  82
  83 // protected
  84
  85 void numeric::copy(numeric const & other)
  86 {
  87     basic::copy(other);
  88     value = new cl_N(*other.value);
  89 }
  90
  91 void numeric::destroy(bool call_parent)
  92 {
  93     delete value;
  94     if (call_parent) basic::destroy(call_parent);
  95 }
  96
  97 //////////
  98 // other constructors
  99 //////////
 100
 101 // public
 102
 103 numeric::numeric(int i) : basic(TINFO_NUMERIC)
 104 {
 105     debugmsg("numeric constructor from int",LOGLEVEL_CONSTRUCT);
 106     // Not the whole int-range is available if we don't cast to long
 107     // first. This is due to the behaviour of the cl_I-ctor, which
 108     // emphasizes efficiency:
 109     value = new cl_I((long) i);
 110     calchash();
 111     setflag(status_flags::evaluated|
 112             status_flags::hash_calculated);
 113 }
 114
 115 numeric::numeric(unsigned int i) : basic(TINFO_NUMERIC)
 116 {
 117     debugmsg("numeric constructor from uint",LOGLEVEL_CONSTRUCT);
 118     // Not the whole uint-range is available if we don't cast to ulong
 119     // first. This is due to the behaviour of the cl_I-ctor, which
 120     // emphasizes efficiency:
 121     value = new cl_I((unsigned long)i);
 122     calchash();
 123     setflag(status_flags::evaluated|
 124             status_flags::hash_calculated);
 125 }
 126
 127 numeric::numeric(long i) : basic(TINFO_NUMERIC)
 128 {
 129     debugmsg("numeric constructor from long",LOGLEVEL_CONSTRUCT);
 130     value = new cl_I(i);
 131     calchash();
 132     setflag(status_flags::evaluated|
 133             status_flags::hash_calculated);
 134 }
 135
 136 numeric::numeric(unsigned long i) : basic(TINFO_NUMERIC)
 137 {
 138     debugmsg("numeric constructor from ulong",LOGLEVEL_CONSTRUCT);
 139     value = new cl_I(i);
 140     calchash();
 141     setflag(status_flags::evaluated|
 142             status_flags::hash_calculated);
 143 }
 144
 145 /** Ctor for rational numerics a/b.
 146  *
 147  *  @exception overflow_error (division by zero) */
 148 numeric::numeric(long numer, long denom) : basic(TINFO_NUMERIC)
 149 {
 150     debugmsg("numeric constructor from long/long",LOGLEVEL_CONSTRUCT);
 151     if (!denom)
 152         throw (std::overflow_error("division by zero"));
 153     value = new cl_I(numer);
 154     *value = *value / cl_I(denom);
 155     calchash();
 156     setflag(status_flags::evaluated|
 157             status_flags::hash_calculated);
 158 }
 159
 160 numeric::numeric(double d) : basic(TINFO_NUMERIC)
 161 {
 162     debugmsg("numeric constructor from double",LOGLEVEL_CONSTRUCT);
 163     // We really want to explicitly use the type cl_LF instead of the
 164     // more general cl_F, since that would give us a cl_DF only which
 165     // will not be promoted to cl_LF if overflow occurs:
 166     value = new cl_N;
 167     *value = cl_float(d, cl_default_float_format);
 168     calchash();
 169     setflag(status_flags::evaluated|
 170             status_flags::hash_calculated);
 171 }
 172
 173 numeric::numeric(char const *s) : basic(TINFO_NUMERIC)
 174 {   // MISSING: treatment of complex and ints and rationals.
 175     debugmsg("numeric constructor from string",LOGLEVEL_CONSTRUCT);
 176     if (strchr(s, '.'))
 177         value = new cl_LF(s);
 178     else
 179         value = new cl_I(s);
 180     calchash();
 181     setflag(status_flags::evaluated|
 182             status_flags::hash_calculated);
 183 }
 184
 185 /** Ctor from CLN types.  This is for the initiated user or internal use
 186  *  only. */
 187 numeric::numeric(cl_N const & z) : basic(TINFO_NUMERIC)
 188 {
 189     debugmsg("numeric constructor from cl_N", LOGLEVEL_CONSTRUCT);
 190     value = new cl_N(z);
 191     calchash();
 192     setflag(status_flags::evaluated|
 193             status_flags::hash_calculated);
 194 }
 195
 196 //////////
 197 // functions overriding virtual functions from bases classes
 198 //////////
 199
 200 // public
 201
 202 basic * numeric::duplicate() const
 203 {
 204     debugmsg("numeric duplicate", LOGLEVEL_DUPLICATE);
 205     return new numeric(*this);
 206 }
 207
 208 // The method printraw doesn't do much, it simply uses CLN's operator<<() for
 209 // output, which is ugly but reliable. Examples:
 210 // 2+2i
 211 void numeric::printraw(ostream & os) const
 212 {
 213     debugmsg("numeric printraw", LOGLEVEL_PRINT);
 214     os << "numeric(" << *value << ")";
 215 }
 216
 217 // The method print adds to the output so it blends more consistently together
 218 // with the other routines and produces something compatible to Maple input.
 219 void numeric::print(ostream & os, unsigned upper_precedence) const
 220 {
 221     debugmsg("numeric print", LOGLEVEL_PRINT);
 222     if (is_real()) {
 223         // case 1, real:  x  or  -x
 224         if ((precedence<=upper_precedence) && (!is_pos_integer())) {
 225             os << "(" << *value << ")";
 226         } else {
 227             os << *value;
 228         }
 229     } else {
 230         // case 2, imaginary:  y*I  or  -y*I
 231         if (realpart(*value) == 0) {
 232             if ((precedence<=upper_precedence) && (imagpart(*value) < 0)) {
 233                 if (imagpart(*value) == -1) {
 234                     os << "(-I)";
 235                 } else {
 236                     os << "(" << imagpart(*value) << "*I)";
 237                 }
 238             } else {
 239                 if (imagpart(*value) == 1) {
 240                     os << "I";
 241                 } else {
 242                     if (imagpart (*value) == -1) {
 243                         os << "-I";
 244                     } else {
 245                         os << imagpart(*value) << "*I";
 246                     }
 247                 }
 248             }
 249         } else {
 250             // case 3, complex:  x+y*I  or  x-y*I  or  -x+y*I  or  -x-y*I
 251             if (precedence <= upper_precedence) os << "(";
 252             os << realpart(*value);
 253             if (imagpart(*value) < 0) {
 254                 if (imagpart(*value) == -1) {
 255                     os << "-I";
 256                 } else {
 257                     os << imagpart(*value) << "*I";
 258                 }
 259             } else {
 260                 if (imagpart(*value) == 1) {
 261                     os << "+I";
 262                 } else {
 263                     os << "+" << imagpart(*value) << "*I";
 264                 }
 265             }
 266             if (precedence <= upper_precedence) os << ")";
 267         }
 268     }
 269 }
 270
 271 bool numeric::info(unsigned inf) const
 272 {
 273     switch (inf) {
 274     case info_flags::numeric:
 275     case info_flags::polynomial:
 276     case info_flags::rational_function:
 277         return true;
 278     case info_flags::real:
 279         return is_real();
 280     case info_flags::rational:
 281     case info_flags::rational_polynomial:
 282         return is_rational();
 283     case info_flags::integer:
 284     case info_flags::integer_polynomial:
 285         return is_integer();
 286     case info_flags::positive:
 287         return is_positive();
 288     case info_flags::negative:
 289         return is_negative();
 290     case info_flags::nonnegative:
 291         return compare(numZERO())>=0;
 292     case info_flags::posint:
 293         return is_pos_integer();
 294     case info_flags::negint:
 295         return is_integer() && (compare(numZERO())<0);
 296     case info_flags::nonnegint:
 297         return is_nonneg_integer();
 298     case info_flags::even:
 299         return is_even();
 300     case info_flags::odd:
 301         return is_odd();
 302     case info_flags::prime:
 303         return is_prime();
 304     }
 305     return false;
 306 }
 307
 308 /** Cast numeric into a floating-point object.  For example exact numeric(1) is
 309  *  returned as a 1.0000000000000000000000 and so on according to how Digits is
 310  *  currently set.
 311  *
 312  *  @param level  ignored, but needed for overriding basic::evalf.
 313  *  @return an ex-handle to a numeric. */
 314 ex numeric::evalf(int level) const
 315 {
 316     // level can safely be discarded for numeric objects.
 317     return numeric(cl_float(1.0, cl_default_float_format) * (*value));  // -> CLN
 318 }
 319
 320 // protected
 321
 322 int numeric::compare_same_type(basic const & other) const
 323 {
 324     ASSERT(is_exactly_of_type(other, numeric));
 325     numeric const & o = static_cast<numeric &>(const_cast<basic &>(other));
 326
 327     if (*value == *o.value) {
 328         return 0;
 329     }
 330
 331     return compare(o);
 332 }
 333
 334 bool numeric::is_equal_same_type(basic const & other) const
 335 {
 336     ASSERT(is_exactly_of_type(other,numeric));
 337     numeric const *o = static_cast<numeric const *>(&other);
 338
 339     return is_equal(*o);
 340 }
 341
 342 /*
 343 unsigned numeric::calchash(void) const
 344 {
 345     double d=to_double();
 346     int s=d>0 ? 1 : -1;
 347     d=fabs(d);
 348     if (d>0x07FF0000) {
 349         d=0x07FF0000;
 350     }
 351     return 0x88000000U+s*unsigned(d/0x07FF0000);
 352 }
 353 */
 354
 355
 356 //////////
 357 // new virtual functions which can be overridden by derived classes
 358 //////////
 359
 360 // none
 361
 362 //////////
 363 // non-virtual functions in this class
 364 //////////
 365
 366 // public
 367
 368 /** Numerical addition method.  Adds argument to *this and returns result as
 369  *  a new numeric object. */
 370 numeric numeric::add(numeric const & other) const
 371 {
 372     return numeric((*value)+(*other.value));
 373 }
 374
 375 /** Numerical subtraction method.  Subtracts argument from *this and returns
 376  *  result as a new numeric object. */
 377 numeric numeric::sub(numeric const & other) const
 378 {
 379     return numeric((*value)-(*other.value));
 380 }
 381
 382 /** Numerical multiplication method.  Multiplies *this and argument and returns
 383  *  result as a new numeric object. */
 384 numeric numeric::mul(numeric const & other) const
 385 {
 386     static const numeric * numONEp=&numONE();
 387     if (this==numONEp) {
 388         return other;
 389     } else if (&other==numONEp) {
 390         return *this;
 391     }
 392     return numeric((*value)*(*other.value));
 393 }
 394
 395 /** Numerical division method.  Divides *this by argument and returns result as
 396  *  a new numeric object.
 397  *
 398  *  @exception overflow_error (division by zero) */
 399 numeric numeric::div(numeric const & other) const
 400 {
 401     if (zerop(*other.value))
 402         throw (std::overflow_error("division by zero"));
 403     return numeric((*value)/(*other.value));
 404 }
 405
 406 numeric numeric::power(numeric const & other) const
 407 {
 408     static const numeric * numONEp=&numONE();
 409     if (&other==numONEp) {
 410         return *this;
 411     }
 412     if (zerop(*value) && other.is_real() && minusp(realpart(*other.value)))
 413         throw (std::overflow_error("division by zero"));
 414     return numeric(expt(*value,*other.value));
 415 }
 416
 417 /** Inverse of a number. */
 418 numeric numeric::inverse(void) const
 419 {
 420     return numeric(recip(*value));  // -> CLN
 421 }
 422
 423 numeric const & numeric::add_dyn(numeric const & other) const
 424 {
 425     return static_cast<numeric const &>((new numeric((*value)+(*other.value)))->
 426                                         setflag(status_flags::dynallocated));
 427 }
 428
 429 numeric const & numeric::sub_dyn(numeric const & other) const
 430 {
 431     return static_cast<numeric const &>((new numeric((*value)-(*other.value)))->
 432                                         setflag(status_flags::dynallocated));
 433 }
 434
 435 numeric const & numeric::mul_dyn(numeric const & other) const
 436 {
 437     static const numeric * numONEp=&numONE();
 438     if (this==numONEp) {
 439         return other;
 440     } else if (&other==numONEp) {
 441         return *this;
 442     }
 443     return static_cast<numeric const &>((new numeric((*value)*(*other.value)))->
 444                                         setflag(status_flags::dynallocated));
 445 }
 446
 447 numeric const & numeric::div_dyn(numeric const & other) const
 448 {
 449     if (zerop(*other.value))
 450         throw (std::overflow_error("division by zero"));
 451     return static_cast<numeric const &>((new numeric((*value)/(*other.value)))->
 452                                         setflag(status_flags::dynallocated));
 453 }
 454
 455 numeric const & numeric::power_dyn(numeric const & other) const
 456 {
 457     static const numeric * numONEp=&numONE();
 458     if (&other==numONEp) {
 459         return *this;
 460     }
 461     // The ifs are only a workaround for a bug in CLN. It gets stuck otherwise:
 462     if ( !other.is_integer() &&
 463          other.is_rational() &&
 464          (*this).is_nonneg_integer() ) {
 465         if ( !zerop(*value) ) {
 466             return static_cast<numeric const &>((new numeric(exp(*other.value * log(*value))))->
 467                                                 setflag(status_flags::dynallocated));
 468         } else {
 469             if ( !zerop(*other.value) ) {  // 0^(n/m)
 470                 return static_cast<numeric const &>((new numeric(0))->
 471                                                     setflag(status_flags::dynallocated));
 472             } else {                       // raise FPE (0^0 requested)
 473                 return static_cast<numeric const &>((new numeric(1/(*other.value)))->
 474                                                     setflag(status_flags::dynallocated));
 475             }
 476         }
 477     } else {                               // default -> CLN
 478         return static_cast<numeric const &>((new numeric(expt(*value,*other.value)))->
 479                                             setflag(status_flags::dynallocated));
 480     }
 481 }
 482
 483 numeric const & numeric::operator=(int i)
 484 {
 485     return operator=(numeric(i));
 486 }
 487
 488 numeric const & numeric::operator=(unsigned int i)
 489 {
 490     return operator=(numeric(i));
 491 }
 492
 493 numeric const & numeric::operator=(long i)
 494 {
 495     return operator=(numeric(i));
 496 }
 497
 498 numeric const & numeric::operator=(unsigned long i)
 499 {
 500     return operator=(numeric(i));
 501 }
 502
 503 numeric const & numeric::operator=(double d)
 504 {
 505     return operator=(numeric(d));
 506 }
 507
 508 numeric const & numeric::operator=(char const * s)
 509 {
 510     return operator=(numeric(s));
 511 }
 512
 513 /** This method establishes a canonical order on all numbers.  For complex
 514  *  numbers this is not possible in a mathematically consistent way but we need
 515  *  to establish some order and it ought to be fast.  So we simply define it
 516  *  similar to Maple's csgn. */
 517 int numeric::compare(numeric const & other) const
 518 {
 519     // Comparing two real numbers?
 520     if (is_real() && other.is_real())
 521         // Yes, just compare them
 522         return cl_compare(The(cl_R)(*value), The(cl_R)(*other.value));
 523     else {
 524         // No, first compare real parts
 525         cl_signean real_cmp = cl_compare(realpart(*value), realpart(*other.value));
 526         if (real_cmp)
 527             return real_cmp;
 528
 529         return cl_compare(imagpart(*value), imagpart(*other.value));
 530     }
 531 }
 532
 533 bool numeric::is_equal(numeric const & other) const
 534 {
 535     return (*value == *other.value);
 536 }
 537
 538 /** True if object is zero. */
 539 bool numeric::is_zero(void) const
 540 {
 541     return zerop(*value);  // -> CLN
 542 }
 543
 544 /** True if object is not complex and greater than zero. */
 545 bool numeric::is_positive(void) const
 546 {
 547     if (is_real()) {
 548         return plusp(The(cl_R)(*value));  // -> CLN
 549     }
 550     return false;
 551 }
 552
 553 /** True if object is not complex and less than zero. */
 554 bool numeric::is_negative(void) const
 555 {
 556     if (is_real()) {
 557         return minusp(The(cl_R)(*value));  // -> CLN
 558     }
 559     return false;
 560 }
 561
 562 /** True if object is a non-complex integer. */
 563 bool numeric::is_integer(void) const
 564 {
 565     return (bool)instanceof(*value, cl_I_ring);  // -> CLN
 566 }
 567
 568 /** True if object is an exact integer greater than zero. */
 569 bool numeric::is_pos_integer(void) const
 570 {
 571     return (is_integer() &&
 572             plusp(The(cl_I)(*value)));  // -> CLN
 573 }
 574
 575 /** True if object is an exact integer greater or equal zero. */
 576 bool numeric::is_nonneg_integer(void) const
 577 {
 578     return (is_integer() &&
 579             !minusp(The(cl_I)(*value)));  // -> CLN
 580 }
 581
 582 /** True if object is an exact even integer. */
 583 bool numeric::is_even(void) const
 584 {
 585     return (is_integer() &&
 586             evenp(The(cl_I)(*value)));  // -> CLN
 587 }
 588
 589 /** True if object is an exact odd integer. */
 590 bool numeric::is_odd(void) const
 591 {
 592     return (is_integer() &&
 593             oddp(The(cl_I)(*value)));  // -> CLN
 594 }
 595
 596 /** Probabilistic primality test.
 597  *
 598  *  @return  true if object is exact integer and prime. */
 599 bool numeric::is_prime(void) const
 600 {
 601     return (is_integer() &&
 602             isprobprime(The(cl_I)(*value)));  // -> CLN
 603 }
 604
 605 /** True if object is an exact rational number, may even be complex
 606  *  (denominator may be unity). */
 607 bool numeric::is_rational(void) const
 608 {
 609     if (instanceof(*value, cl_RA_ring)) {
 610         return true;
 611     } else if (!is_real()) {  // complex case, handle Q(i):
 612         if ( instanceof(realpart(*value), cl_RA_ring) &&
 613              instanceof(imagpart(*value), cl_RA_ring) )
 614             return true;
 615     }
 616     return false;
 617 }
 618
 619 /** True if object is a real integer, rational or float (but not complex). */
 620 bool numeric::is_real(void) const
 621 {
 622     return (bool)instanceof(*value, cl_R_ring);  // -> CLN
 623 }
 624
 625 bool numeric::operator==(numeric const & other) const
 626 {
 627     return (*value == *other.value);  // -> CLN
 628 }
 629
 630 bool numeric::operator!=(numeric const & other) const
 631 {
 632     return (*value != *other.value);  // -> CLN
 633 }
 634
 635 /** Numerical comparison: less.
 636  *
 637  *  @exception invalid_argument (complex inequality) */
 638 bool numeric::operator<(numeric const & other) const
 639 {
 640     if ( is_real() && other.is_real() ) {
 641         return (bool)(The(cl_R)(*value) < The(cl_R)(*other.value));  // -> CLN
 642     }
 643     throw (std::invalid_argument("numeric::operator<(): complex inequality"));
 644     return false;  // make compiler shut up
 645 }
 646
 647 /** Numerical comparison: less or equal.
 648  *
 649  *  @exception invalid_argument (complex inequality) */
 650 bool numeric::operator<=(numeric const & other) const
 651 {
 652     if ( is_real() && other.is_real() ) {
 653         return (bool)(The(cl_R)(*value) <= The(cl_R)(*other.value));  // -> CLN
 654     }
 655     throw (std::invalid_argument("numeric::operator<=(): complex inequality"));
 656     return false;  // make compiler shut up
 657 }
 658
 659 /** Numerical comparison: greater.
 660  *
 661  *  @exception invalid_argument (complex inequality) */
 662 bool numeric::operator>(numeric const & other) const
 663 {
 664     if ( is_real() && other.is_real() ) {
 665         return (bool)(The(cl_R)(*value) > The(cl_R)(*other.value));  // -> CLN
 666     }
 667     throw (std::invalid_argument("numeric::operator>(): complex inequality"));
 668     return false;  // make compiler shut up
 669 }
 670
 671 /** Numerical comparison: greater or equal.
 672  *
 673  *  @exception invalid_argument (complex inequality) */
 674 bool numeric::operator>=(numeric const & other) const
 675 {
 676     if ( is_real() && other.is_real() ) {
 677         return (bool)(The(cl_R)(*value) >= The(cl_R)(*other.value));  // -> CLN
 678     }
 679     throw (std::invalid_argument("numeric::operator>=(): complex inequality"));
 680     return false;  // make compiler shut up
 681 }
 682
 683 /** Converts numeric types to machine's int. You should check with is_integer()
 684  *  if the number is really an integer before calling this method. */
 685 int numeric::to_int(void) const
 686 {
 687     ASSERT(is_integer());
 688     return cl_I_to_int(The(cl_I)(*value));
 689 }
 690
 691 /** Converts numeric types to machine's double. You should check with is_real()
 692  *  if the number is really not complex before calling this method. */
 693 double numeric::to_double(void) const
 694 {
 695     ASSERT(is_real());
 696     return cl_double_approx(realpart(*value));
 697 }
 698
 699 /** Real part of a number. */
 700 numeric numeric::real(void) const
 701 {
 702     return numeric(realpart(*value));  // -> CLN
 703 }
 704
 705 /** Imaginary part of a number. */
 706 numeric numeric::imag(void) const
 707 {
 708     return numeric(imagpart(*value));  // -> CLN
 709 }
 710
 711 #ifndef SANE_LINKER
 712 // Unfortunately, CLN did not provide an official way to access the numerator
 713 // or denominator of a rational number (cl_RA). Doing some excavations in CLN
 714 // one finds how it works internally in src/rational/cl_RA.h:
 715 struct cl_heap_ratio : cl_heap {
 716     cl_I numerator;
 717     cl_I denominator;
 718 };
 719
 720 inline cl_heap_ratio* TheRatio (const cl_N& obj)
 721 { return (cl_heap_ratio*)(obj.pointer); }
 722 #endif // ndef SANE_LINKER
 723
 724 /** Numerator.  Computes the numerator of rational numbers, rationalized
 725  *  numerator of complex if real and imaginary part are both rational numbers
 726  *  (i.e numer(4/3+5/6*I) == 8+5*I), the number itself in all other cases. */
 727 numeric numeric::numer(void) const
 728 {
 729     if (is_integer()) {
 730         return numeric(*this);
 731     }
 732 #ifdef SANE_LINKER
 733     else if (instanceof(*value, cl_RA_ring)) {
 734         return numeric(numerator(The(cl_RA)(*value)));
 735     }
 736     else if (!is_real()) {  // complex case, handle Q(i):
 737         cl_R r = realpart(*value);
 738         cl_R i = imagpart(*value);
 739         if (instanceof(r, cl_I_ring) && instanceof(i, cl_I_ring))
 740             return numeric(*this);
 741         if (instanceof(r, cl_I_ring) && instanceof(i, cl_RA_ring))
 742             return numeric(complex(r*denominator(The(cl_RA)(i)), numerator(The(cl_RA)(i))));
 743         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_I_ring))
 744             return numeric(complex(numerator(The(cl_RA)(r)), i*denominator(The(cl_RA)(r))));
 745         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_RA_ring)) {
 746             cl_I s = lcm(denominator(The(cl_RA)(r)), denominator(The(cl_RA)(i)));
 747             return numeric(complex(numerator(The(cl_RA)(r))*(exquo(s,denominator(The(cl_RA)(r)))),
 748                                    numerator(The(cl_RA)(i))*(exquo(s,denominator(The(cl_RA)(i))))));
 749         }
 750     }
 751 #else
 752     else if (instanceof(*value, cl_RA_ring)) {
 753         return numeric(TheRatio(*value)->numerator);
 754     }
 755     else if (!is_real()) {  // complex case, handle Q(i):
 756         cl_R r = realpart(*value);
 757         cl_R i = imagpart(*value);
 758         if (instanceof(r, cl_I_ring) && instanceof(i, cl_I_ring))
 759             return numeric(*this);
 760         if (instanceof(r, cl_I_ring) && instanceof(i, cl_RA_ring))
 761             return numeric(complex(r*TheRatio(i)->denominator, TheRatio(i)->numerator));
 762         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_I_ring))
 763             return numeric(complex(TheRatio(r)->numerator, i*TheRatio(r)->denominator));
 764         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_RA_ring)) {
 765             cl_I s = lcm(TheRatio(r)->denominator, TheRatio(i)->denominator);
 766             return numeric(complex(TheRatio(r)->numerator*(exquo(s,TheRatio(r)->denominator)),
 767                                    TheRatio(i)->numerator*(exquo(s,TheRatio(i)->denominator))));
 768         }
 769     }
 770 #endif // def SANE_LINKER
 771     // at least one float encountered
 772     return numeric(*this);
 773 }
 774
 775 /** Denominator.  Computes the denominator of rational numbers, common integer
 776  *  denominator of complex if real and imaginary part are both rational numbers
 777  *  (i.e denom(4/3+5/6*I) == 6), one in all other cases. */
 778 numeric numeric::denom(void) const
 779 {
 780     if (is_integer()) {
 781         return numONE();
 782     }
 783 #ifdef SANE_LINKER
 784     if (instanceof(*value, cl_RA_ring)) {
 785         return numeric(denominator(The(cl_RA)(*value)));
 786     }
 787     if (!is_real()) {  // complex case, handle Q(i):
 788         cl_R r = realpart(*value);
 789         cl_R i = imagpart(*value);
 790         if (instanceof(r, cl_I_ring) && instanceof(i, cl_I_ring))
 791             return numONE();
 792         if (instanceof(r, cl_I_ring) && instanceof(i, cl_RA_ring))
 793             return numeric(denominator(The(cl_RA)(i)));
 794         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_I_ring))
 795             return numeric(denominator(The(cl_RA)(r)));
 796         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_RA_ring))
 797             return numeric(lcm(denominator(The(cl_RA)(r)), denominator(The(cl_RA)(i))));
 798     }
 799 #else
 800     if (instanceof(*value, cl_RA_ring)) {
 801         return numeric(TheRatio(*value)->denominator);
 802     }
 803     if (!is_real()) {  // complex case, handle Q(i):
 804         cl_R r = realpart(*value);
 805         cl_R i = imagpart(*value);
 806         if (instanceof(r, cl_I_ring) && instanceof(i, cl_I_ring))
 807             return numONE();
 808         if (instanceof(r, cl_I_ring) && instanceof(i, cl_RA_ring))
 809             return numeric(TheRatio(i)->denominator);
 810         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_I_ring))
 811             return numeric(TheRatio(r)->denominator);
 812         if (instanceof(r, cl_RA_ring) && instanceof(i, cl_RA_ring))
 813             return numeric(lcm(TheRatio(r)->denominator, TheRatio(i)->denominator));
 814     }
 815 #endif // def SANE_LINKER
 816     // at least one float encountered
 817     return numONE();
 818 }
 819
 820 /** Size in binary notation.  For integers, this is the smallest n >= 0 such
 821  *  that -2^n <= x < 2^n. If x > 0, this is the unique n > 0 such that
 822  *  2^(n-1) <= x < 2^n.
 823  *
 824  *  @return  number of bits (excluding sign) needed to represent that number
 825  *  in two's complement if it is an integer, 0 otherwise. */
 826 int numeric::int_length(void) const
 827 {
 828     if (is_integer()) {
 829         return integer_length(The(cl_I)(*value));  // -> CLN
 830     } else {
 831         return 0;
 832     }
 833 }
 834
 835
 836 //////////
 837 // static member variables
 838 //////////
 839
 840 // protected
 841
 842 unsigned numeric::precedence = 30;
 843
 844 //////////
 845 // global constants
 846 //////////
 847
 848 const numeric some_numeric;
 849 type_info const & typeid_numeric=typeid(some_numeric);
 850 /** Imaginary unit.  This is not a constant but a numeric since we are
 851  *  natively handing complex numbers anyways. */
 852 const numeric I = (complex(cl_I(0),cl_I(1)));
 853
 854 //////////
 855 // global functions
 856 //////////
 857
 858 numeric const & numZERO(void)
 859 {
 860     const static ex eZERO = ex((new numeric(0))->setflag(status_flags::dynallocated));
 861     const static numeric * nZERO = static_cast<const numeric *>(eZERO.bp);
 862     return *nZERO;
 863 }
 864
 865 numeric const & numONE(void)
 866 {
 867     const static ex eONE = ex((new numeric(1))->setflag(status_flags::dynallocated));
 868     const static numeric * nONE = static_cast<const numeric *>(eONE.bp);
 869     return *nONE;
 870 }
 871
 872 numeric const & numTWO(void)
 873 {
 874     const static ex eTWO = ex((new numeric(2))->setflag(status_flags::dynallocated));
 875     const static numeric * nTWO = static_cast<const numeric *>(eTWO.bp);
 876     return *nTWO;
 877 }
 878
 879 numeric const & numTHREE(void)
 880 {
 881     const static ex eTHREE = ex((new numeric(3))->setflag(status_flags::dynallocated));
 882     const static numeric * nTHREE = static_cast<const numeric *>(eTHREE.bp);
 883     return *nTHREE;
 884 }
 885
 886 numeric const & numMINUSONE(void)
 887 {
 888     const static ex eMINUSONE = ex((new numeric(-1))->setflag(status_flags::dynallocated));
 889     const static numeric * nMINUSONE = static_cast<const numeric *>(eMINUSONE.bp);
 890     return *nMINUSONE;
 891 }
 892
 893 numeric const & numHALF(void)
 894 {
 895     const static ex eHALF = ex((new numeric(1, 2))->setflag(status_flags::dynallocated));
 896     const static numeric * nHALF = static_cast<const numeric *>(eHALF.bp);
 897     return *nHALF;
 898 }
 899
 900 /** Exponential function.
 901  *
 902  *  @return  arbitrary precision numerical exp(x). */
 903 numeric exp(numeric const & x)
 904 {
 905     return exp(*x.value);  // -> CLN
 906 }
 907
 908 /** Natural logarithm.
 909  *
 910  *  @param z complex number
 911  *  @return  arbitrary precision numerical log(x).
 912  *  @exception overflow_error (logarithmic singularity) */
 913 numeric log(numeric const & z)
 914 {
 915     if (z.is_zero())
 916         throw (std::overflow_error("log(): logarithmic singularity"));
 917     return log(*z.value);  // -> CLN
 918 }
 919
 920 /** Numeric sine (trigonometric function).
 921  *
 922  *  @return  arbitrary precision numerical sin(x). */
 923 numeric sin(numeric const & x)
 924 {
 925     return sin(*x.value);  // -> CLN
 926 }
 927
 928 /** Numeric cosine (trigonometric function).
 929  *
 930  *  @return  arbitrary precision numerical cos(x). */
 931 numeric cos(numeric const & x)
 932 {
 933     return cos(*x.value);  // -> CLN
 934 }
 935
 936 /** Numeric tangent (trigonometric function).
 937  *
 938  *  @return  arbitrary precision numerical tan(x). */
 939 numeric tan(numeric const & x)
 940 {
 941     return tan(*x.value);  // -> CLN
 942 }
 943
 944 /** Numeric inverse sine (trigonometric function).
 945  *
 946  *  @return  arbitrary precision numerical asin(x). */
 947 numeric asin(numeric const & x)
 948 {
 949     return asin(*x.value);  // -> CLN
 950 }
 951
 952 /** Numeric inverse cosine (trigonometric function).
 953  *
 954  *  @return  arbitrary precision numerical acos(x). */
 955 numeric acos(numeric const & x)
 956 {
 957     return acos(*x.value);  // -> CLN
 958 }
 959
 960 /** Arcustangents.
 961  *
 962  *  @param z complex number
 963  *  @return atan(z)
 964  *  @exception overflow_error (logarithmic singularity) */
 965 numeric atan(numeric const & x)
 966 {
 967     if (!x.is_real() &&
 968         x.real().is_zero() &&
 969         !abs(x.imag()).is_equal(numONE()))
 970         throw (std::overflow_error("atan(): logarithmic singularity"));
 971     return atan(*x.value);  // -> CLN
 972 }
 973
 974 /** Arcustangents.
 975  *
 976  *  @param x real number
 977  *  @param y real number
 978  *  @return atan(y/x) */
 979 numeric atan(numeric const & y, numeric const & x)
 980 {
 981     if (x.is_real() && y.is_real())
 982         return atan(realpart(*x.value), realpart(*y.value));  // -> CLN
 983     else
 984         throw (std::invalid_argument("numeric::atan(): complex argument"));
 985 }
 986
 987 /** Numeric hyperbolic sine (trigonometric function).
 988  *
 989  *  @return  arbitrary precision numerical sinh(x). */
 990 numeric sinh(numeric const & x)
 991 {
 992     return sinh(*x.value);  // -> CLN
 993 }
 994
 995 /** Numeric hyperbolic cosine (trigonometric function).
 996  *
 997  *  @return  arbitrary precision numerical cosh(x). */
 998 numeric cosh(numeric const & x)
 999 {
1000     return cosh(*x.value);  // -> CLN
1001 }
1002
1003 /** Numeric hyperbolic tangent (trigonometric function).
1004  *
1005  *  @return  arbitrary precision numerical tanh(x). */
1006 numeric tanh(numeric const & x)
1007 {
1008     return tanh(*x.value);  // -> CLN
1009 }
1010
1011 /** Numeric inverse hyperbolic sine (trigonometric function).
1012  *
1013  *  @return  arbitrary precision numerical asinh(x). */
1014 numeric asinh(numeric const & x)
1015 {
1016     return asinh(*x.value);  // -> CLN
1017 }
1018
1019 /** Numeric inverse hyperbolic cosine (trigonometric function).
1020  *
1021  *  @return  arbitrary precision numerical acosh(x). */
1022 numeric acosh(numeric const & x)
1023 {
1024     return acosh(*x.value);  // -> CLN
1025 }
1026
1027 /** Numeric inverse hyperbolic tangent (trigonometric function).
1028  *
1029  *  @return  arbitrary precision numerical atanh(x). */
1030 numeric atanh(numeric const & x)
1031 {
1032     return atanh(*x.value);  // -> CLN
1033 }
1034
1035 /** The gamma function.
1036  *  stub stub stub stub stub stub! */
1037 numeric gamma(numeric const & x)
1038 {
1039     clog << "gamma(): Nobody expects the Spanish inquisition" << endl;
1040     return numeric(0);
1041 }
1042
1043 /** Factorial combinatorial function.
1044  *
1045  *  @exception range_error (argument must be integer >= 0) */
1046 numeric factorial(numeric const & nn)
1047 {
1048     if ( !nn.is_nonneg_integer() ) {
1049         throw (std::range_error("numeric::factorial(): argument must be integer >= 0"));
1050     }
1051
1052     return numeric(factorial(nn.to_int()));  // -> CLN
1053 }
1054
1055 /** The double factorial combinatorial function.  (Scarcely used, but still
1056  *  useful in cases, like for exact results of Gamma(n+1/2) for instance.)
1057  *
1058  *  @param n  integer argument >= -1
1059  *  @return n!! == n * (n-2) * (n-4) * ... * ({1|2}) with 0!! == 1 == (-1)!!
1060  *  @exception range_error (argument must be integer >= -1) */
1061 numeric doublefactorial(numeric const & nn)
1062 {
1063     // We store the results separately for even and odd arguments.  This has
1064     // the advantage that we don't have to compute any even result at all if
1065     // the function is always called with odd arguments and vice versa.  There
1066     // is no tradeoff involved in this, it is guaranteed to save time as well
1067     // as memory.  (If this is not enough justification consider the Gamma
1068     // function of half integer arguments: it only needs odd doublefactorials.)
1069     static vector<numeric> evenresults;
1070     static int highest_evenresult = -1;
1071     static vector<numeric> oddresults;
1072     static int highest_oddresult = -1;
1073
1074     if ( nn == numeric(-1) ) {
1075         return numONE();
1076     }
1077     if ( !nn.is_nonneg_integer() ) {
1078         throw (std::range_error("numeric::doublefactorial(): argument must be integer >= -1"));
1079     }
1080     if ( nn.is_even() ) {
1081         int n = nn.div(numTWO()).to_int();
1082         if ( n <= highest_evenresult ) {
1083             return evenresults[n];
1084         }
1085         if ( evenresults.capacity() < (unsigned)(n+1) ) {
1086             evenresults.reserve(n+1);
1087         }
1088         if ( highest_evenresult < 0 ) {
1089             evenresults.push_back(numONE());
1090             highest_evenresult=0;
1091         }
1092         for (int i=highest_evenresult+1; i<=n; i++) {
1093             evenresults.push_back(numeric(evenresults[i-1].mul(numeric(i*2))));
1094         }
1095         highest_evenresult=n;
1096         return evenresults[n];
1097     } else {
1098         int n = nn.sub(numONE()).div(numTWO()).to_int();
1099         if ( n <= highest_oddresult ) {
1100             return oddresults[n];
1101         }
1102         if ( oddresults.capacity() < (unsigned)n ) {
1103             oddresults.reserve(n+1);
1104         }
1105         if ( highest_oddresult < 0 ) {
1106             oddresults.push_back(numONE());
1107             highest_oddresult=0;
1108         }
1109         for (int i=highest_oddresult+1; i<=n; i++) {
1110             oddresults.push_back(numeric(oddresults[i-1].mul(numeric(i*2+1))));
1111         }
1112         highest_oddresult=n;
1113         return oddresults[n];
1114     }
1115 }
1116
1117 /** The Binomial function. It computes the binomial coefficients. If the
1118  *  arguments are both nonnegative integers and 0 <= k <= n, then
1119  *  binomial(n, k) = n!/k!/(n-k)! which is the number of ways of choosing k
1120  *  objects from n distinct objects. If k > n, then binomial(n,k) returns 0. */
1121 numeric binomial(numeric const & n, numeric const & k)
1122 {
1123     if (n.is_nonneg_integer() && k.is_nonneg_integer()) {
1124         return numeric(binomial(n.to_int(),k.to_int()));  // -> CLN
1125     } else {
1126         // should really be gamma(n+1)/(gamma(r+1)/gamma(n-r+1)
1127         return numeric(0);
1128     }
1129     // return factorial(n).div(factorial(k).mul(factorial(n.sub(k))));
1130 }
1131
1132 /** Absolute value. */
1133 numeric abs(numeric const & x)
1134 {
1135     return abs(*x.value);  // -> CLN
1136 }
1137
1138 /** Modulus (in positive representation).
1139  *  In general, mod(a,b) has the sign of b or is zero, and rem(a,b) has the
1140  *  sign of a or is zero. This is different from Maple's modp, where the sign
1141  *  of b is ignored. It is in agreement with Mathematica's Mod.
1142  *
1143  *  @return a mod b in the range [0,abs(b)-1] with sign of b if both are
1144  *  integer, 0 otherwise. */
1145 numeric mod(numeric const & a, numeric const & b)
1146 {
1147     if (a.is_integer() && b.is_integer()) {
1148         return mod(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
1149     }
1150     else {
1151         return numZERO();  // Throw?
1152     }
1153 }
1154
1155 /** Modulus (in symmetric representation).
1156  *  Equivalent to Maple's mods.
1157  *
1158  *  @return a mod b in the range [-iquo(abs(m)-1,2), iquo(abs(m),2)]. */
1159 numeric smod(numeric const & a, numeric const & b)
1160 {
1161     if (a.is_integer() && b.is_integer()) {
1162         cl_I b2 = The(cl_I)(ceiling1(The(cl_I)(*b.value) / 2)) - 1;
1163         return mod(The(cl_I)(*a.value) + b2, The(cl_I)(*b.value)) - b2;
1164     } else {
1165         return numZERO();  // Throw?
1166     }
1167 }
1168
1169 /** Numeric integer remainder.
1170  *  Equivalent to Maple's irem(a,b) as far as sign conventions are concerned.
1171  *  In general, mod(a,b) has the sign of b or is zero, and irem(a,b) has the
1172  *  sign of a or is zero.
1173  *
1174  *  @return remainder of a/b if both are integer, 0 otherwise. */
1175 numeric irem(numeric const & a, numeric const & b)
1176 {
1177     if (a.is_integer() && b.is_integer()) {
1178         return rem(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
1179     }
1180     else {
1181         return numZERO();  // Throw?
1182     }
1183 }
1184
1185 /** Numeric integer remainder.
1186  *  Equivalent to Maple's irem(a,b,'q') it obeyes the relation
1187  *  irem(a,b,q) == a - q*b.  In general, mod(a,b) has the sign of b or is zero,
1188  *  and irem(a,b) has the sign of a or is zero.
1189  *
1190  *  @return remainder of a/b and quotient stored in q if both are integer,
1191  *  0 otherwise. */
1192 numeric irem(numeric const & a, numeric const & b, numeric & q)
1193 {
1194     if (a.is_integer() && b.is_integer()) {  // -> CLN
1195         cl_I_div_t rem_quo = truncate2(The(cl_I)(*a.value), The(cl_I)(*b.value));
1196         q = rem_quo.quotient;
1197         return rem_quo.remainder;
1198     }
1199     else {
1200         q = numZERO();
1201         return numZERO();  // Throw?
1202     }
1203 }
1204
1205 /** Numeric integer quotient.
1206  *  Equivalent to Maple's iquo as far as sign conventions are concerned.
1207  *
1208  *  @return truncated quotient of a/b if both are integer, 0 otherwise. */
1209 numeric iquo(numeric const & a, numeric const & b)
1210 {
1211     if (a.is_integer() && b.is_integer()) {
1212         return truncate1(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
1213     } else {
1214         return numZERO();  // Throw?
1215     }
1216 }
1217
1218 /** Numeric integer quotient.
1219  *  Equivalent to Maple's iquo(a,b,'r') it obeyes the relation
1220  *  r == a - iquo(a,b,r)*b.
1221  *
1222  *  @return truncated quotient of a/b and remainder stored in r if both are
1223  *  integer, 0 otherwise. */
1224 numeric iquo(numeric const & a, numeric const & b, numeric & r)
1225 {
1226     if (a.is_integer() && b.is_integer()) {  // -> CLN
1227         cl_I_div_t rem_quo = truncate2(The(cl_I)(*a.value), The(cl_I)(*b.value));
1228         r = rem_quo.remainder;
1229         return rem_quo.quotient;
1230     } else {
1231         r = numZERO();
1232         return numZERO();  // Throw?
1233     }
1234 }
1235
1236 /** Numeric square root.
1237  *  If possible, sqrt(z) should respect squares of exact numbers, i.e. sqrt(4)
1238  *  should return integer 2.
1239  *
1240  *  @param z numeric argument
1241  *  @return square root of z. Branch cut along negative real axis, the negative
1242  *  real axis itself where imag(z)==0 and real(z)<0 belongs to the upper part
1243  *  where imag(z)>0. */
1244 numeric sqrt(numeric const & z)
1245 {
1246     return sqrt(*z.value);  // -> CLN
1247 }
1248
1249 /** Integer numeric square root. */
1250 numeric isqrt(numeric const & x)
1251 {
1252         if (x.is_integer()) {
1253                 cl_I root;
1254                 isqrt(The(cl_I)(*x.value), &root);      // -> CLN
1255                 return root;
1256         } else
1257                 return numZERO();  // Throw?
1258 }
1259
1260 /** Greatest Common Divisor.
1261  *
1262  *  @return  The GCD of two numbers if both are integer, a numerical 1
1263  *  if they are not. */
1264 numeric gcd(numeric const & a, numeric const & b)
1265 {
1266     if (a.is_integer() && b.is_integer())
1267         return gcd(The(cl_I)(*a.value), The(cl_I)(*b.value));   // -> CLN
1268     else
1269         return numONE();
1270 }
1271
1272 /** Least Common Multiple.
1273  *
1274  *  @return  The LCM of two numbers if both are integer, the product of those
1275  *  two numbers if they are not. */
1276 numeric lcm(numeric const & a, numeric const & b)
1277 {
1278     if (a.is_integer() && b.is_integer())
1279         return lcm(The(cl_I)(*a.value), The(cl_I)(*b.value));   // -> CLN
1280     else
1281         return *a.value * *b.value;
1282 }
1283
1284 ex PiEvalf(void)
1285 {
1286     return numeric(cl_pi(cl_default_float_format));  // -> CLN
1287 }
1288
1289 ex EulerGammaEvalf(void)
1290 {
1291     return numeric(cl_eulerconst(cl_default_float_format));  // -> CLN
1292 }
1293
1294 ex CatalanEvalf(void)
1295 {
1296     return numeric(cl_catalanconst(cl_default_float_format));  // -> CLN
1297 }
1298
1299 // It initializes to 17 digits, because in CLN cl_float_format(17) turns out to
1300 // be 61 (<64) while cl_float_format(18)=65.  We want to have a cl_LF instead
1301 // of cl_SF, cl_FF or cl_DF but everything else is basically arbitrary.
1302 _numeric_digits::_numeric_digits()
1303     : digits(17)
1304 {
1305     assert(!too_late);
1306     too_late = true;
1307     cl_default_float_format = cl_float_format(17);
1308 }
1309
1310 _numeric_digits& _numeric_digits::operator=(long prec)
1311 {
1312     digits=prec;
1313     cl_default_float_format = cl_float_format(prec);
1314     return *this;
1315 }
1316
1317 _numeric_digits::operator long()
1318 {
1319     return (long)digits;
1320 }
1321
1322 void _numeric_digits::print(ostream & os) const
1323 {
1324     debugmsg("_numeric_digits print", LOGLEVEL_PRINT);
1325     os << digits;
1326 }
1327
1328 ostream& operator<<(ostream& os, _numeric_digits const & e)
1329 {
1330     e.print(os);
1331     return os;
1332 }
1333
1334 //////////
1335 // static member variables
1336 //////////
1337
1338 // private
1339
1340 bool _numeric_digits::too_late = false;
1341
1342 /** Accuracy in decimal digits.  Only object of this type!  Can be set using
1343  *  assignment from C++ unsigned ints and evaluated like any built-in type. */
1344 _numeric_digits Digits;