]> www.ginac.de Git - ginac.git/blobdiff - ginac/numeric.cpp
gcc 2.95 doesn't have the <limits> header
[ginac.git] / ginac / numeric.cpp
index 0274d63af6fcee55c222c507253a7f124a53d8c5..0aeef27a69a6a9249712d6485f05e6cc6143437a 100644 (file)
@@ -7,7 +7,7 @@
  *  of special functions or implement the interface to the bignum package. */
 
 /*
- *  GiNaC Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -313,7 +313,7 @@ DEFAULT_UNARCHIVE(numeric)
  *  want to visibly distinguish from cl_LF.
  *
  *  @see numeric::print() */
-static void print_real_number(const print_context & c, const cln::cl_R &x)
+static void print_real_number(const print_context & c, const cln::cl_R & x)
 {
        cln::cl_print_flags ourflags;
        if (cln::instanceof(x, cln::cl_RA_ring)) {
@@ -322,8 +322,10 @@ static void print_real_number(const print_context & c, const cln::cl_R &x)
                    !is_a<print_latex>(c)) {
                        cln::print_real(c.s, ourflags, x);
                } else {  // rational output in LaTeX context
+                       if (x < 0)
+                               c.s << "-";
                        c.s << "\\frac{";
-                       cln::print_real(c.s, ourflags, cln::numerator(cln::the<cln::cl_RA>(x)));
+                       cln::print_real(c.s, ourflags, cln::abs(cln::numerator(cln::the<cln::cl_RA>(x))));
                        c.s << "}{";
                        cln::print_real(c.s, ourflags, cln::denominator(cln::the<cln::cl_RA>(x)));
                        c.s << '}';
@@ -337,6 +339,82 @@ static void print_real_number(const print_context & c, const cln::cl_R &x)
        }
 }
 
+/** Helper function to print integer number in C++ source format.
+ *
+ *  @see numeric::print() */
+static void print_integer_csrc(const print_context & c, const cln::cl_I & x)
+{
+       // Print small numbers in compact float format, but larger numbers in
+       // scientific format
+       const int max_cln_int = 536870911; // 2^29-1
+       if (x >= cln::cl_I(-max_cln_int) && x <= cln::cl_I(max_cln_int))
+               c.s << cln::cl_I_to_int(x) << ".0";
+       else
+               c.s << cln::double_approx(x);
+}
+
+/** Helper function to print real number in C++ source format.
+ *
+ *  @see numeric::print() */
+static void print_real_csrc(const print_context & c, const cln::cl_R & x)
+{
+       if (cln::instanceof(x, cln::cl_I_ring)) {
+
+               // Integer number
+               print_integer_csrc(c, cln::the<cln::cl_I>(x));
+
+       } else if (cln::instanceof(x, cln::cl_RA_ring)) {
+
+               // Rational number
+               const cln::cl_I numer = cln::numerator(cln::the<cln::cl_RA>(x));
+               const cln::cl_I denom = cln::denominator(cln::the<cln::cl_RA>(x));
+               if (cln::plusp(x) > 0) {
+                       c.s << "(";
+                       print_integer_csrc(c, numer);
+               } else {
+                       c.s << "-(";
+                       print_integer_csrc(c, -numer);
+               }
+               c.s << "/";
+               print_integer_csrc(c, denom);
+               c.s << ")";
+
+       } else {
+
+               // Anything else
+               c.s << cln::double_approx(x);
+       }
+}
+
+/** Helper function to print real number in C++ source format using cl_N types.
+ *
+ *  @see numeric::print() */
+static void print_real_cl_N(const print_context & c, const cln::cl_R & x)
+{
+       if (cln::instanceof(x, cln::cl_I_ring)) {
+
+               // Integer number
+               c.s << "cln::cl_I(\"";
+               print_real_number(c, x);
+               c.s << "\")";
+
+       } else if (cln::instanceof(x, cln::cl_RA_ring)) {
+
+               // Rational number
+               cln::cl_print_flags ourflags;
+               c.s << "cln::cl_RA(\"";
+               cln::print_rational(c.s, ourflags, cln::the<cln::cl_RA>(x));
+               c.s << "\")";
+
+       } else {
+
+               // Anything else
+               c.s << "cln::cl_F(\"";
+               print_real_number(c, cln::cl_float(1.0, cln::default_float_format) * x);
+               c.s << "_" << Digits << "\")";
+       }
+}
+
 /** This method adds to the output so it blends more consistently together
  *  with the other routines and produces something compatible to ginsh input.
  *  
@@ -350,51 +428,69 @@ void numeric::print(const print_context & c, unsigned level) const
                    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
                    << std::endl;
 
+       } else if (is_a<print_csrc_cl_N>(c)) {
+
+               // CLN output
+               if (this->is_real()) {
+
+                       // Real number
+                       print_real_cl_N(c, cln::the<cln::cl_R>(value));
+
+               } else {
+
+                       // Complex number
+                       c.s << "cln::complex(";
+                       print_real_cl_N(c, cln::realpart(cln::the<cln::cl_N>(value)));
+                       c.s << ",";
+                       print_real_cl_N(c, cln::imagpart(cln::the<cln::cl_N>(value)));
+                       c.s << ")";
+               }
+
        } else if (is_a<print_csrc>(c)) {
 
+               // C++ source output
                std::ios::fmtflags oldflags = c.s.flags();
                c.s.setf(std::ios::scientific);
                int oldprec = c.s.precision();
+
+               // Set precision
                if (is_a<print_csrc_double>(c))
                        c.s.precision(16);
                else
                        c.s.precision(7);
-               if (this->is_rational() && !this->is_integer()) {
-                       if (compare(_num0) > 0) {
-                               c.s << "(";
-                               if (is_a<print_csrc_cl_N>(c))
-                                       c.s << "cln::cl_F(\"" << numer().evalf() << "\")";
-                               else
-                                       c.s << numer().to_double();
-                       } else {
-                               c.s << "-(";
-                               if (is_a<print_csrc_cl_N>(c))
-                                       c.s << "cln::cl_F(\"" << -numer().evalf() << "\")";
-                               else
-                                       c.s << -numer().to_double();
-                       }
-                       c.s << "/";
-                       if (is_a<print_csrc_cl_N>(c))
-                               c.s << "cln::cl_F(\"" << denom().evalf() << "\")";
-                       else
-                               c.s << denom().to_double();
-                       c.s << ")";
+
+               if (this->is_real()) {
+
+                       // Real number
+                       print_real_csrc(c, cln::the<cln::cl_R>(value));
+
                } else {
-                       if (is_a<print_csrc_cl_N>(c))
-                               c.s << "cln::cl_F(\"" << evalf() << "\")";
+
+                       // Complex number
+                       c.s << "std::complex<";
+                       if (is_a<print_csrc_double>(c))
+                               c.s << "double>(";
                        else
-                               c.s << to_double();
+                               c.s << "float>(";
+
+                       print_real_csrc(c, cln::realpart(cln::the<cln::cl_N>(value)));
+                       c.s << ",";
+                       print_real_csrc(c, cln::imagpart(cln::the<cln::cl_N>(value)));
+                       c.s << ")";
                }
+
                c.s.flags(oldflags);
                c.s.precision(oldprec);
 
        } else {
+
                const std::string par_open  = is_a<print_latex>(c) ? "{(" : "(";
                const std::string par_close = is_a<print_latex>(c) ? ")}" : ")";
                const std::string imag_sym  = is_a<print_latex>(c) ? "i" : "I";
                const std::string mul_sym   = is_a<print_latex>(c) ? " " : "*";
                const cln::cl_R r = cln::realpart(cln::the<cln::cl_N>(value));
                const cln::cl_R i = cln::imagpart(cln::the<cln::cl_N>(value));
+
                if (is_a<print_python_repr>(c))
                        c.s << class_name() << "('";
                if (cln::zerop(i)) {
@@ -1688,9 +1784,12 @@ const numeric smod(const numeric &a, const numeric &b)
  *  In general, mod(a,b) has the sign of b or is zero, and irem(a,b) has the
  *  sign of a or is zero.
  *
- *  @return remainder of a/b if both are integer, 0 otherwise. */
+ *  @return remainder of a/b if both are integer, 0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric irem(const numeric &a, const numeric &b)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::irem(): division by zero");
        if (a.is_integer() && b.is_integer())
                return cln::rem(cln::the<cln::cl_I>(a.to_cl_N()),
                                cln::the<cln::cl_I>(b.to_cl_N()));
@@ -1705,9 +1804,12 @@ const numeric irem(const numeric &a, const numeric &b)
  *  and irem(a,b) has the sign of a or is zero.  
  *
  *  @return remainder of a/b and quotient stored in q if both are integer,
- *  0 otherwise. */
+ *  0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric irem(const numeric &a, const numeric &b, numeric &q)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::irem(): division by zero");
        if (a.is_integer() && b.is_integer()) {
                const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
                                                               cln::the<cln::cl_I>(b.to_cl_N()));
@@ -1723,9 +1825,12 @@ const numeric irem(const numeric &a, const numeric &b, numeric &q)
 /** Numeric integer quotient.
  *  Equivalent to Maple's iquo as far as sign conventions are concerned.
  *  
- *  @return truncated quotient of a/b if both are integer, 0 otherwise. */
+ *  @return truncated quotient of a/b if both are integer, 0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric iquo(const numeric &a, const numeric &b)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::iquo(): division by zero");
        if (a.is_integer() && b.is_integer())
                return cln::truncate1(cln::the<cln::cl_I>(a.to_cl_N()),
                                  cln::the<cln::cl_I>(b.to_cl_N()));
@@ -1739,9 +1844,12 @@ const numeric iquo(const numeric &a, const numeric &b)
  *  r == a - iquo(a,b,r)*b.
  *
  *  @return truncated quotient of a/b and remainder stored in r if both are
- *  integer, 0 otherwise. */
+ *  integer, 0 otherwise.
+ *  @exception overflow_error (division by zero) if b is zero. */
 const numeric iquo(const numeric &a, const numeric &b, numeric &r)
 {
+       if (b.is_zero())
+               throw std::overflow_error("numeric::iquo(): division by zero");
        if (a.is_integer() && b.is_integer()) {
                const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
                                                               cln::the<cln::cl_I>(b.to_cl_N()));