]> www.ginac.de Git - ginac.git/blobdiff - ginac/numeric.cpp
[BUGFIX] numeric::info(nonnegative): correctly handle complex numbers.
[ginac.git] / ginac / numeric.cpp
index 3ca4b58c5ef31487c1e7a654a47885e67df7c4aa..18725b216cf5e40eeccebe3cc7fe41abec698372 100644 (file)
@@ -7,7 +7,7 @@
  *  of special functions or implement the interface to the bignum package. */
 
 /*
- *  GiNaC Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2011 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
 
 #include "config.h"
 
-#include <vector>
-#include <stdexcept>
-#include <string>
-#include <sstream>
-#include <limits>
-
 #include "numeric.h"
 #include "ex.h"
 #include "operators.h"
 #include "tostring.h"
 #include "utils.h"
 
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
 // CLN should pollute the global namespace as little as possible.  Hence, we
 // include most of it here and include only the part needed for properly
 // declaring cln::cl_number in numeric.h.  This can only be safely done in
@@ -73,7 +73,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(numeric, basic,
 //////////
 
 /** default ctor. Numerically it initializes to an integer zero. */
-numeric::numeric() : basic(&numeric::tinfo_static)
+numeric::numeric()
 {
        value = cln::cl_I(0);
        setflag(status_flags::evaluated | status_flags::expanded);
@@ -85,7 +85,7 @@ numeric::numeric() : basic(&numeric::tinfo_static)
 
 // public
 
-numeric::numeric(int i) : basic(&numeric::tinfo_static)
+numeric::numeric(int i)
 {
        // 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
@@ -106,7 +106,7 @@ numeric::numeric(int i) : basic(&numeric::tinfo_static)
 }
 
 
-numeric::numeric(unsigned int i) : basic(&numeric::tinfo_static)
+numeric::numeric(unsigned int i)
 {
        // 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
@@ -127,14 +127,14 @@ numeric::numeric(unsigned int i) : basic(&numeric::tinfo_static)
 }
 
 
-numeric::numeric(long i) : basic(&numeric::tinfo_static)
+numeric::numeric(long i)
 {
        value = cln::cl_I(i);
        setflag(status_flags::evaluated | status_flags::expanded);
 }
 
 
-numeric::numeric(unsigned long i) : basic(&numeric::tinfo_static)
+numeric::numeric(unsigned long i)
 {
        value = cln::cl_I(i);
        setflag(status_flags::evaluated | status_flags::expanded);
@@ -144,7 +144,7 @@ numeric::numeric(unsigned long i) : basic(&numeric::tinfo_static)
 /** Constructor for rational numerics a/b.
  *
  *  @exception overflow_error (division by zero) */
-numeric::numeric(long numer, long denom) : basic(&numeric::tinfo_static)
+numeric::numeric(long numer, long denom)
 {
        if (!denom)
                throw std::overflow_error("division by zero");
@@ -153,7 +153,7 @@ numeric::numeric(long numer, long denom) : basic(&numeric::tinfo_static)
 }
 
 
-numeric::numeric(double d) : basic(&numeric::tinfo_static)
+numeric::numeric(double d)
 {
        // We really want to explicitly use the type cl_LF instead of the
        // more general cl_F, since that would give us a cl_DF only which
@@ -165,7 +165,7 @@ numeric::numeric(double d) : basic(&numeric::tinfo_static)
 
 /** ctor from C-style string.  It also accepts complex numbers in GiNaC
  *  notation like "2+5*I". */
-numeric::numeric(const char *s) : basic(&numeric::tinfo_static)
+numeric::numeric(const char *s)
 {
        cln::cl_N ctorval = 0;
        // parse complex numbers (functional but not completely safe, unfortunately
@@ -244,7 +244,7 @@ numeric::numeric(const char *s) : basic(&numeric::tinfo_static)
 
 /** Ctor from CLN types.  This is for the initiated user or internal use
  *  only. */
-numeric::numeric(const cln::cl_N &z) : basic(&numeric::tinfo_static)
+numeric::numeric(const cln::cl_N &z)
 {
        value = z;
        setflag(status_flags::evaluated | status_flags::expanded);
@@ -255,66 +255,126 @@ numeric::numeric(const cln::cl_N &z) : basic(&numeric::tinfo_static)
 // archiving
 //////////
 
-numeric::numeric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+/** 
+ * Construct a floating point number from sign, mantissa, and exponent 
+ */
+static const cln::cl_F make_real_float(const cln::cl_idecoded_float& dec)
 {
-       cln::cl_N ctorval = 0;
+       cln::cl_F x = cln::cl_float(dec.mantissa, cln::default_float_format);
+       x = cln::scale_float(x, dec.exponent);
+       cln::cl_F sign = cln::cl_float(dec.sign, cln::default_float_format);
+       x = cln::float_sign(sign, x);
+       return x;
+}
 
+/** 
+ * Read serialized floating point number 
+ */
+static const cln::cl_F read_real_float(std::istream& s)
+{
+       cln::cl_idecoded_float dec;
+       s >> dec.sign >> dec.mantissa >> dec.exponent;
+       const cln::cl_F x = make_real_float(dec);
+       return x;
+}
+
+void numeric::read_archive(const archive_node &n, lst &sym_lst)
+{
+       inherited::read_archive(n, sym_lst);
+       value = 0;
+       
        // Read number as string
        std::string str;
        if (n.find_string("number", str)) {
                std::istringstream s(str);
-               cln::cl_idecoded_float re, im;
+               cln::cl_R re, im;
                char c;
                s.get(c);
                switch (c) {
-                       case 'R':    // Integer-decoded real number
-                               s >> re.sign >> re.mantissa >> re.exponent;
-                               ctorval = re.sign * re.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), re.exponent);
+                       case 'R':
+                               // real FP (floating point) number
+                               re = read_real_float(s);
+                               value = re;
+                               break;
+                       case 'C':
+                               // both real and imaginary part are FP numbers
+                               re = read_real_float(s);
+                               im = read_real_float(s); 
+                               value = cln::complex(re, im);
                                break;
-                       case 'C':    // Integer-decoded complex number
-                               s >> re.sign >> re.mantissa >> re.exponent;
-                               s >> im.sign >> im.mantissa >> im.exponent;
-                               ctorval = cln::complex(re.sign * re.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), re.exponent),
-                                                      im.sign * im.mantissa * cln::expt(cln::cl_float(2.0, cln::default_float_format), im.exponent));
+                       case 'H':
+                               // real part is a rational number,
+                               // imaginary part is a FP number
+                               s >> re;
+                               im = read_real_float(s);
+                               value = cln::complex(re, im);
                                break;
-                       default:    // Ordinary number
+                       case 'J':
+                               // real part is a FP number,
+                               // imaginary part is a rational number
+                               re = read_real_float(s);
+                               s >> im;
+                               value = cln::complex(re, im);
+                               break;
+                       default:
+                               // both real and imaginary parts are rational
                                s.putback(c);
-                               s >> ctorval;
+                               s >> value;
                                break;
                }
        }
-       value = ctorval;
        setflag(status_flags::evaluated | status_flags::expanded);
 }
+GINAC_BIND_UNARCHIVER(numeric);
+
+static void write_real_float(std::ostream& s, const cln::cl_R& n)
+{
+       const cln::cl_idecoded_float dec = cln::integer_decode_float(cln::the<cln::cl_F>(n));
+       s << dec.sign << ' ' << dec.mantissa << ' ' << dec.exponent;
+}
 
 void numeric::archive(archive_node &n) const
 {
        inherited::archive(n);
 
        // Write number as string
+       
+       const cln::cl_R re = cln::realpart(value);
+       const cln::cl_R im = cln::imagpart(value);
+       const bool re_rationalp = cln::instanceof(re, cln::cl_RA_ring);
+       const bool im_rationalp = cln::instanceof(im, cln::cl_RA_ring);
+
+       // Non-rational numbers are written in an integer-decoded format
+       // to preserve the precision
        std::ostringstream s;
-       if (this->is_crational())
+       if (re_rationalp && im_rationalp)
                s << value;
-       else {
-               // Non-rational numbers are written in an integer-decoded format
-               // to preserve the precision
-               if (this->is_real()) {
-                       cln::cl_idecoded_float re = cln::integer_decode_float(cln::the<cln::cl_F>(value));
-                       s << "R";
-                       s << re.sign << " " << re.mantissa << " " << re.exponent;
-               } else {
-                       cln::cl_idecoded_float re = cln::integer_decode_float(cln::the<cln::cl_F>(cln::realpart(cln::the<cln::cl_N>(value))));
-                       cln::cl_idecoded_float im = cln::integer_decode_float(cln::the<cln::cl_F>(cln::imagpart(cln::the<cln::cl_N>(value))));
-                       s << "C";
-                       s << re.sign << " " << re.mantissa << " " << re.exponent << " ";
-                       s << im.sign << " " << im.mantissa << " " << im.exponent;
-               }
+       else if (zerop(im)) {
+               // real FP (floating point) number
+               s << 'R';
+               write_real_float(s, re);
+       } else if (re_rationalp) {
+               s << 'H'; // just any unique character
+               // real part is a rational number,
+               // imaginary part is a FP number
+               s << re << ' ';
+               write_real_float(s, im);
+       } else if (im_rationalp) {
+               s << 'J';
+               // real part is a FP number,
+               // imaginary part is a rational number
+               write_real_float(s, re);
+               s << ' ' << im;
+       } else  {
+               // both real and imaginary parts are floating point
+               s << 'C';
+               write_real_float(s, re);
+               s << ' ';
+               write_real_float(s, im);
        }
        n.add_string("number", s.str());
 }
 
-DEFAULT_UNARCHIVE(numeric)
-
 //////////
 // functions overriding virtual functions from base classes
 //////////
@@ -408,7 +468,7 @@ static inline bool coerce(T1& dst, const T2& arg);
  * @sa http://www.ginac.de/pipermail/cln-list/2006-October/000248.html
  */
 template<>
-static inline bool coerce<int, cln::cl_I>(int& dst, const cln::cl_I& arg)
+inline bool coerce<int, cln::cl_I>(int& dst, const cln::cl_I& arg)
 {
        static const cln::cl_I cl_max_int = 
                (cln::cl_I)(long)(std::numeric_limits<int>::max());
@@ -422,7 +482,7 @@ static inline bool coerce<int, cln::cl_I>(int& dst, const cln::cl_I& arg)
 }
 
 template<>
-static inline bool coerce<unsigned int, cln::cl_I>(unsigned int& dst, const cln::cl_I& arg)
+inline bool coerce<unsigned int, cln::cl_I>(unsigned int& dst, const cln::cl_I& arg)
 {
        static const cln::cl_I cl_max_uint = 
                (cln::cl_I)(unsigned long)(std::numeric_limits<unsigned int>::max());
@@ -640,7 +700,7 @@ bool numeric::info(unsigned inf) const
                case info_flags::negative:
                        return is_negative();
                case info_flags::nonnegative:
-                       return !is_negative();
+                       return is_zero() || is_positive();
                case info_flags::posint:
                        return is_pos_integer();
                case info_flags::negint:
@@ -1379,7 +1439,7 @@ const numeric I = numeric(cln::complex(cln::cl_I(0),cln::cl_I(1)));
  *  @return  arbitrary precision numerical exp(x). */
 const numeric exp(const numeric &x)
 {
-       return cln::exp(x.to_cl_N());
+       return numeric(cln::exp(x.to_cl_N()));
 }
 
 
@@ -1392,7 +1452,7 @@ const numeric log(const numeric &x)
 {
        if (x.is_zero())
                throw pole_error("log(): logarithmic pole",0);
-       return cln::log(x.to_cl_N());
+       return numeric(cln::log(x.to_cl_N()));
 }
 
 
@@ -1401,7 +1461,7 @@ const numeric log(const numeric &x)
  *  @return  arbitrary precision numerical sin(x). */
 const numeric sin(const numeric &x)
 {
-       return cln::sin(x.to_cl_N());
+       return numeric(cln::sin(x.to_cl_N()));
 }
 
 
@@ -1410,7 +1470,7 @@ const numeric sin(const numeric &x)
  *  @return  arbitrary precision numerical cos(x). */
 const numeric cos(const numeric &x)
 {
-       return cln::cos(x.to_cl_N());
+       return numeric(cln::cos(x.to_cl_N()));
 }
 
 
@@ -1419,7 +1479,7 @@ const numeric cos(const numeric &x)
  *  @return  arbitrary precision numerical tan(x). */
 const numeric tan(const numeric &x)
 {
-       return cln::tan(x.to_cl_N());
+       return numeric(cln::tan(x.to_cl_N()));
 }
        
 
@@ -1428,7 +1488,7 @@ const numeric tan(const numeric &x)
  *  @return  arbitrary precision numerical asin(x). */
 const numeric asin(const numeric &x)
 {
-       return cln::asin(x.to_cl_N());
+       return numeric(cln::asin(x.to_cl_N()));
 }
 
 
@@ -1437,7 +1497,7 @@ const numeric asin(const numeric &x)
  *  @return  arbitrary precision numerical acos(x). */
 const numeric acos(const numeric &x)
 {
-       return cln::acos(x.to_cl_N());
+       return numeric(cln::acos(x.to_cl_N()));
 }
        
 
@@ -1452,7 +1512,7 @@ const numeric atan(const numeric &x)
            x.real().is_zero() &&
            abs(x.imag()).is_equal(*_num1_p))
                throw pole_error("atan(): logarithmic pole",0);
-       return cln::atan(x.to_cl_N());
+       return numeric(cln::atan(x.to_cl_N()));
 }
 
 
@@ -1468,8 +1528,8 @@ const numeric atan(const numeric &y, const numeric &x)
        if (x.is_zero() && y.is_zero())
                return *_num0_p;
        if (x.is_real() && y.is_real())
-               return cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
-                                cln::the<cln::cl_R>(y.to_cl_N()));
+               return numeric(cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
+                                cln::the<cln::cl_R>(y.to_cl_N())));
 
        // Compute -I*log((x+I*y)/sqrt(x^2+y^2))
        //      == -I*log((x+I*y)/sqrt((x+I*y)*(x-I*y)))
@@ -1485,7 +1545,7 @@ const numeric atan(const numeric &y, const numeric &x)
                // x-I*y==0 => y/x==-I, so this is a pole (we have x!=0).
                throw pole_error("atan(): logarithmic pole",0);
        }
-       return cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m));
+       return numeric(cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m)));
 }
 
 
@@ -1494,7 +1554,7 @@ const numeric atan(const numeric &y, const numeric &x)
  *  @return  arbitrary precision numerical sinh(x). */
 const numeric sinh(const numeric &x)
 {
-       return cln::sinh(x.to_cl_N());
+       return numeric(cln::sinh(x.to_cl_N()));
 }
 
 
@@ -1503,7 +1563,7 @@ const numeric sinh(const numeric &x)
  *  @return  arbitrary precision numerical cosh(x). */
 const numeric cosh(const numeric &x)
 {
-       return cln::cosh(x.to_cl_N());
+       return numeric(cln::cosh(x.to_cl_N()));
 }
 
 
@@ -1512,7 +1572,7 @@ const numeric cosh(const numeric &x)
  *  @return  arbitrary precision numerical tanh(x). */
 const numeric tanh(const numeric &x)
 {
-       return cln::tanh(x.to_cl_N());
+       return numeric(cln::tanh(x.to_cl_N()));
 }
        
 
@@ -1521,7 +1581,7 @@ const numeric tanh(const numeric &x)
  *  @return  arbitrary precision numerical asinh(x). */
 const numeric asinh(const numeric &x)
 {
-       return cln::asinh(x.to_cl_N());
+       return numeric(cln::asinh(x.to_cl_N()));
 }
 
 
@@ -1530,7 +1590,7 @@ const numeric asinh(const numeric &x)
  *  @return  arbitrary precision numerical acosh(x). */
 const numeric acosh(const numeric &x)
 {
-       return cln::acosh(x.to_cl_N());
+       return numeric(cln::acosh(x.to_cl_N()));
 }
 
 
@@ -1539,7 +1599,7 @@ const numeric acosh(const numeric &x)
  *  @return  arbitrary precision numerical atanh(x). */
 const numeric atanh(const numeric &x)
 {
-       return cln::atanh(x.to_cl_N());
+       return numeric(cln::atanh(x.to_cl_N()));
 }
 
 
@@ -1611,24 +1671,24 @@ static cln::cl_N Li2_projection(const cln::cl_N &x,
        return Li2_series(x, prec);
 }
 
+
 /** Numeric evaluation of Dilogarithm.  The domain is the entire complex plane,
  *  the branch cut lies along the positive real axis, starting at 1 and
  *  continuous with quadrant IV.
  *
  *  @return  arbitrary precision numerical Li2(x). */
-const numeric Li2(const numeric &x)
+const cln::cl_N Li2_(const cln::cl_N& value)
 {
-       if (x.is_zero())
-               return *_num0_p;
+       if (zerop(value))
+               return 0;
        
        // what is the desired float format?
        // first guess: default format
        cln::float_format_t prec = cln::default_float_format;
-       const cln::cl_N value = x.to_cl_N();
        // second guess: the argument's format
-       if (!x.real().is_rational())
+       if (!instanceof(realpart(value), cln::cl_RA_ring))
                prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
-       else if (!x.imag().is_rational())
+       else if (!instanceof(imagpart(value), cln::cl_RA_ring))
                prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
        
        if (value==1)  // may cause trouble with log(1-x)
@@ -1640,7 +1700,16 @@ const numeric Li2(const numeric &x)
                       - cln::zeta(2, prec)
                       - Li2_projection(cln::recip(value), prec));
        else
-               return Li2_projection(x.to_cl_N(), prec);
+               return Li2_projection(value, prec);
+}
+
+const numeric Li2(const numeric &x)
+{
+       const cln::cl_N x_ = x.to_cl_N();
+       if (zerop(x_))
+               return *_num0_p;
+       const cln::cl_N result = Li2_(x_);
+       return numeric(result);
 }
 
 
@@ -1656,7 +1725,7 @@ const numeric zeta(const numeric &x)
        if (x.is_real()) {
                const int aux = (int)(cln::double_approx(cln::the<cln::cl_R>(x.to_cl_N())));
                if (cln::zerop(x.to_cl_N()-aux))
-                       return cln::zeta(aux);
+                       return numeric(cln::zeta(aux));
        }
        throw dunno();
 }
@@ -1953,24 +2022,34 @@ lanczos_coeffs::lanczos_coeffs()
        coeffs[3].swap(coeffs_120);
 }
 
+static const cln::float_format_t guess_precision(const cln::cl_N& x)
+{
+       cln::float_format_t prec = cln::default_float_format;
+       if (!instanceof(realpart(x), cln::cl_RA_ring))
+               prec = cln::float_format(cln::the<cln::cl_F>(realpart(x)));
+       if (!instanceof(imagpart(x), cln::cl_RA_ring))
+               prec = cln::float_format(cln::the<cln::cl_F>(imagpart(x)));
+       return prec;
+}
 
 /** The Gamma function.
  *  Use the Lanczos approximation. If the coefficients used here are not
  *  sufficiently many or sufficiently accurate, more can be calculated
  *  using the program doc/examples/lanczos.cpp. In that case, be sure to
  *  read the comments in that file. */
-const numeric lgamma(const numeric &x)
+const cln::cl_N lgamma(const cln::cl_N &x)
 {
+       cln::float_format_t prec = guess_precision(x);
        lanczos_coeffs lc;
-       if (lc.sufficiently_accurate(Digits)) {
-               cln::cl_N pi_val = cln::pi(cln::default_float_format);
-               if (x.real() < 0.5)
-                       return log(pi_val) - log(sin(pi_val*x.to_cl_N()))
-                               - lgamma(numeric(1).sub(x));
-               cln::cl_N A = lc.calc_lanczos_A(x.to_cl_N());
-               cln::cl_N temp = x.to_cl_N() + lc.get_order() - cln::cl_N(1)/2;
+       if (lc.sufficiently_accurate(prec)) {
+               cln::cl_N pi_val = cln::pi(prec);
+               if (realpart(x) < 0.5)
+                       return cln::log(pi_val) - cln::log(sin(pi_val*x))
+                               - lgamma(1 - x);
+               cln::cl_N A = lc.calc_lanczos_A(x);
+               cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2;
        cln::cl_N result = log(cln::cl_I(2)*pi_val)/2
-                             + (x.to_cl_N()-cln::cl_N(1)/2)*log(temp)
+                             + (x-cln::cl_N(1)/2)*log(temp)
                              - temp
                              + log(A);
        return result;
@@ -1979,17 +2058,25 @@ const numeric lgamma(const numeric &x)
                throw dunno();
 }
 
-const numeric tgamma(const numeric &x)
+const numeric lgamma(const numeric &x)
+{
+       const cln::cl_N x_ = x.to_cl_N();
+       const cln::cl_N result = lgamma(x_);
+       return numeric(result);
+}
+
+const cln::cl_N tgamma(const cln::cl_N &x)
 {
+       cln::float_format_t prec = guess_precision(x);
        lanczos_coeffs lc;
-       if (lc.sufficiently_accurate(Digits)) {
-               cln::cl_N pi_val = cln::pi(cln::default_float_format);
-               if (x.real() < 0.5)
-                       return pi_val/(sin(pi_val*x))/(tgamma(numeric(1).sub(x)).to_cl_N());
-               cln::cl_N A = lc.calc_lanczos_A(x.to_cl_N());
-               cln::cl_N temp = x.to_cl_N() + lc.get_order() - cln::cl_N(1)/2;
+       if (lc.sufficiently_accurate(prec)) {
+               cln::cl_N pi_val = cln::pi(prec);
+               if (realpart(x) < 0.5)
+                       return pi_val/(cln::sin(pi_val*x))/tgamma(1 - x);
+               cln::cl_N A = lc.calc_lanczos_A(x);
+               cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2;
        cln::cl_N result
-                       = sqrt(cln::cl_I(2)*pi_val) * expt(temp, x.to_cl_N()-cln::cl_N(1)/2)
+                       = sqrt(cln::cl_I(2)*pi_val) * expt(temp, x - cln::cl_N(1)/2)
                          * exp(-temp) * A;
        return result;
        }
@@ -1997,6 +2084,12 @@ const numeric tgamma(const numeric &x)
                throw dunno();
 }
 
+const numeric tgamma(const numeric &x)
+{
+       const cln::cl_N x_ = x.to_cl_N();
+       const cln::cl_N result = tgamma(x_);
+       return numeric(result);
+}
 
 /** The psi function (aka polygamma function).
  *  This is only a stub! */
@@ -2128,7 +2221,7 @@ const numeric bernoulli(const numeric &nn)
                next_r = 4;
        }
        if (n<next_r)
-               return results[n/2-1];
+               return numeric(results[n/2-1]);
 
        results.reserve(n/2);
        for (unsigned p=next_r; p<=n;  p+=2) {
@@ -2153,7 +2246,7 @@ const numeric bernoulli(const numeric &nn)
                results.push_back(-b/(p+1));
        }
        next_r = n+2;
-       return results[n/2-1];
+       return numeric(results[n/2-1]);
 }
 
 
@@ -2185,11 +2278,14 @@ const numeric fibonacci(const numeric &n)
        //      F(2n+2) = F(n+1)*(2*F(n) + F(n+1))
        if (n.is_zero())
                return *_num0_p;
-       if (n.is_negative())
-               if (n.is_even())
+       if (n.is_negative()) {
+               if (n.is_even()) {
                        return -fibonacci(-n);
-               else
+               }
+               else {
                        return fibonacci(-n);
+               }
+       }
        
        cln::cl_I u(0);
        cln::cl_I v(1);
@@ -2210,16 +2306,16 @@ const numeric fibonacci(const numeric &n)
        if (n.is_even())
                // Here we don't use the squaring formula because one multiplication
                // is cheaper than two squarings.
-               return u * ((v << 1) - u);
+               return numeric(u * ((v << 1) - u));
        else
-               return cln::square(u) + cln::square(v);    
+               return numeric(cln::square(u) + cln::square(v)); 
 }
 
 
 /** Absolute value. */
 const numeric abs(const numeric& x)
 {
-       return cln::abs(x.to_cl_N());
+       return numeric(cln::abs(x.to_cl_N()));
 }
 
 
@@ -2233,23 +2329,26 @@ const numeric abs(const numeric& x)
 const numeric mod(const numeric &a, const numeric &b)
 {
        if (a.is_integer() && b.is_integer())
-               return cln::mod(cln::the<cln::cl_I>(a.to_cl_N()),
-                               cln::the<cln::cl_I>(b.to_cl_N()));
+               return numeric(cln::mod(cln::the<cln::cl_I>(a.to_cl_N()),
+                               cln::the<cln::cl_I>(b.to_cl_N())));
        else
                return *_num0_p;
 }
 
 
 /** Modulus (in symmetric representation).
- *  Equivalent to Maple's mods.
  *
- *  @return a mod b in the range [-iquo(abs(b)-1,2), iquo(abs(b),2)]. */
-const numeric smod(const numeric &a, const numeric &b)
-{
-       if (a.is_integer() && b.is_integer()) {
-               const cln::cl_I b2 = cln::ceiling1(cln::the<cln::cl_I>(b.to_cl_N()) >> 1) - 1;
-               return cln::mod(cln::the<cln::cl_I>(a.to_cl_N()) + b2,
-                               cln::the<cln::cl_I>(b.to_cl_N())) - b2;
+ *  @return a mod b in the range [-iquo(abs(b),2), iquo(abs(b),2)]. */
+const numeric smod(const numeric &a_, const numeric &b_)
+{
+       if (a_.is_integer() && b_.is_integer()) {
+               const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N());
+               const cln::cl_I b = cln::the<cln::cl_I>(b_.to_cl_N());
+               const cln::cl_I b2 = b >> 1;
+               const cln::cl_I m = cln::mod(a, b);
+               const cln::cl_I m_b = m - b;
+               const cln::cl_I ret = m > b2 ? m_b : m;
+               return numeric(ret);
        } else
                return *_num0_p;
 }
@@ -2267,8 +2366,8 @@ 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()));
+               return numeric(cln::rem(cln::the<cln::cl_I>(a.to_cl_N()),
+                               cln::the<cln::cl_I>(b.to_cl_N())));
        else
                return *_num0_p;
 }
@@ -2289,8 +2388,8 @@ const numeric irem(const numeric &a, const numeric &b, numeric &q)
        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()));
-               q = rem_quo.quotient;
-               return rem_quo.remainder;
+               q = numeric(rem_quo.quotient);
+               return numeric(rem_quo.remainder);
        } else {
                q = *_num0_p;
                return *_num0_p;
@@ -2308,8 +2407,8 @@ 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()));
+               return numeric(cln::truncate1(cln::the<cln::cl_I>(a.to_cl_N()),
+                                 cln::the<cln::cl_I>(b.to_cl_N())));
        else
                return *_num0_p;
 }
@@ -2329,8 +2428,8 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r)
        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()));
-               r = rem_quo.remainder;
-               return rem_quo.quotient;
+               r = numeric(rem_quo.remainder);
+               return numeric(rem_quo.quotient);
        } else {
                r = *_num0_p;
                return *_num0_p;
@@ -2345,8 +2444,8 @@ const numeric iquo(const numeric &a, const numeric &b, numeric &r)
 const numeric gcd(const numeric &a, const numeric &b)
 {
        if (a.is_integer() && b.is_integer())
-               return cln::gcd(cln::the<cln::cl_I>(a.to_cl_N()),
-                               cln::the<cln::cl_I>(b.to_cl_N()));
+               return numeric(cln::gcd(cln::the<cln::cl_I>(a.to_cl_N()),
+                               cln::the<cln::cl_I>(b.to_cl_N())));
        else
                return *_num1_p;
 }
@@ -2359,8 +2458,8 @@ const numeric gcd(const numeric &a, const numeric &b)
 const numeric lcm(const numeric &a, const numeric &b)
 {
        if (a.is_integer() && b.is_integer())
-               return cln::lcm(cln::the<cln::cl_I>(a.to_cl_N()),
-                               cln::the<cln::cl_I>(b.to_cl_N()));
+               return numeric(cln::lcm(cln::the<cln::cl_I>(a.to_cl_N()),
+                               cln::the<cln::cl_I>(b.to_cl_N())));
        else
                return a.mul(b);
 }
@@ -2376,7 +2475,7 @@ const numeric lcm(const numeric &a, const numeric &b)
  *  where imag(x)>0. */
 const numeric sqrt(const numeric &x)
 {
-       return cln::sqrt(x.to_cl_N());
+       return numeric(cln::sqrt(x.to_cl_N()));
 }
 
 
@@ -2386,7 +2485,7 @@ const numeric isqrt(const numeric &x)
        if (x.is_integer()) {
                cln::cl_I root;
                cln::isqrt(cln::the<cln::cl_I>(x.to_cl_N()), &root);
-               return root;
+               return numeric(root);
        } else
                return *_num0_p;
 }