]> www.ginac.de Git - ginac.git/blobdiff - ginac/numeric.cpp
* irem(), iquo(): throw an exception, when second argument vanishes.
[ginac.git] / ginac / numeric.cpp
index e621607a730cabb8878cd5b5cf7d2c0b0dd71b54..858cbbd578ee0ba2fd332876249ddef43e8c026a 100644 (file)
@@ -7,7 +7,7 @@
  *  of special functions or implement the interface to the bignum package. */
 
 /*
- *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2002 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
@@ -90,10 +90,10 @@ numeric::numeric(int i) : basic(TINFO_numeric)
 {
        // Not the whole int-range is available if we don't cast to long
        // first.  This is due to the behaviour of the cl_I-ctor, which
-       // emphasizes efficiency.  However, if the integer is small enough
-       // i.e. satisfies cl_immediate_p(), we save space and dereferences by
-       // using an immediate type:
-       if (cln::cl_immediate_p(i))
+       // emphasizes efficiency.  However, if the integer is small enough
+       // we save space and dereferences by using an immediate type.
+       // (C.f. <cln/object.h>)
+       if (i < (1U<<cl_value_len-1))
                value = cln::cl_I(i);
        else
                value = cln::cl_I((long) i);
@@ -105,10 +105,10 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 {
        // Not the whole uint-range is available if we don't cast to ulong
        // first.  This is due to the behaviour of the cl_I-ctor, which
-       // emphasizes efficiency.  However, if the integer is small enough
-       // i.e. satisfies cl_immediate_p(), we save space and dereferences by
-       // using an immediate type:
-       if (cln::cl_immediate_p(i))
+       // emphasizes efficiency.  However, if the integer is small enough
+       // we save space and dereferences by using an immediate type.
+       // (C.f. <cln/object.h>)
+       if (i < (1U<<cl_value_len-1))
                value = cln::cl_I(i);
        else
                value = cln::cl_I((unsigned long) i);
@@ -354,6 +354,11 @@ void numeric::print(const print_context & c, unsigned level) const
 
                std::ios::fmtflags oldflags = c.s.flags();
                c.s.setf(std::ios::scientific);
+               int oldprec = c.s.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 << "(";
@@ -381,6 +386,7 @@ void numeric::print(const print_context & c, unsigned level) const
                                c.s << to_double();
                }
                c.s.flags(oldflags);
+               c.s.precision(oldprec);
 
        } else {
                const std::string par_open  = is_a<print_latex>(c) ? "{(" : "(";
@@ -403,25 +409,19 @@ void numeric::print(const print_context & c, unsigned level) const
                } else {
                        if (cln::zerop(r)) {
                                // case 2, imaginary:  y*I  or  -y*I
-                               if ((precedence() <= level) && (i < 0)) {
-                                       if (i == -1) {
-                                               c.s << par_open+imag_sym+par_close;
-                                       } else {
+                               if (i==1)
+                                       c.s << imag_sym;
+                               else {
+                                       if (precedence()<=level)
                                                c.s << par_open;
+                                       if (i == -1)
+                                               c.s << "-" << imag_sym;
+                                       else {
                                                print_real_number(c, i);
-                                               c.s << mul_sym+imag_sym+par_close;
-                                       }
-                               } else {
-                                       if (i == 1) {
-                                               c.s << imag_sym;
-                                       } else {
-                                               if (i == -1) {
-                                                       c.s << "-" << imag_sym;
-                                               } else {
-                                                       print_real_number(c, i);
-                                                       c.s << mul_sym+imag_sym;
-                                               }
+                                               c.s << mul_sym+imag_sym;
                                        }
+                                       if (precedence()<=level)
+                                               c.s << par_close;
                                }
                        } else {
                                // case 3, complex:  x+y*I  or  x-y*I  or  -x+y*I  or  -x-y*I
@@ -498,6 +498,21 @@ bool numeric::info(unsigned inf) const
        return false;
 }
 
+int numeric::degree(const ex & s) const
+{
+       return 0;
+}
+
+int numeric::ldegree(const ex & s) const
+{
+       return 0;
+}
+
+ex numeric::coeff(const ex & s, int n) const
+{
+       return n==0 ? *this : _ex0;
+}
+
 /** Disassemble real part and imaginary part to scan for the occurrence of a
  *  single number.  Also handles the imaginary unit.  It ignores the sign on
  *  both this and the argument, which may lead to what might appear as funny
@@ -1540,21 +1555,33 @@ const numeric bernoulli(const numeric &nn)
                results.push_back(cln::recip(cln::cl_RA(6)));
                next_r = 4;
        }
+       if (n<next_r)
+               return results[n/2-1];
+
+       results.reserve(n/2);
        for (unsigned p=next_r; p<=n;  p+=2) {
-               cln::cl_I  c = 1;
+               cln::cl_I  c = 1;  // seed for binonmial coefficients
                cln::cl_RA b = cln::cl_RA(1-p)/2;
                const unsigned p3 = p+3;
-               const unsigned p2 = p+2;
                const unsigned pm = p-2;
-               unsigned i, k;
-               for (i=2, k=0; i <= pm; i += 2, k++) {
-                       c = cln::exquo(c * ((p3 - i)*(p2 - i)), (i - 1)*i);
-                       b = b + c * results[k];
-               }
-               results.push_back(-b / (p+1));
-               next_r += 2;
+               unsigned i, k, p_2;
+               // test if intermediate unsigned int can be represented by immediate
+               // objects by CLN (i.e. < 2^29 for 32 Bit machines, see <cln/object.h>)
+               if (p < (1UL<<cl_value_len/2)) {
+                       for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+                               c = cln::exquo(c * ((p3-i) * p_2), (i-1)*k);
+                               b = b + c*results[k-1];
+                       }
+               } else {
+                       for (i=2, k=1, p_2=p/2; i<=pm; i+=2, ++k, --p_2) {
+                               c = cln::exquo((c * (p3-i)) * p_2, cln::cl_I(i-1)*k);
+                               b = b + c*results[k-1];
+                       }
+               }
+               results.push_back(-b/(p+1));
        }
-       return results[n/2 - 1];
+       next_r = n+2;
+       return results[n/2-1];
 }
 
 
@@ -1661,9 +1688,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()));
@@ -1678,9 +1708,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()));
@@ -1696,9 +1729,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()));
@@ -1712,9 +1748,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()));