]> www.ginac.de Git - ginac.git/blobdiff - ginac/numeric.cpp
Fix MSC compiler warning.
[ginac.git] / ginac / numeric.cpp
index 5649353f308f68a0c26e715034b84977c5702979..948409e334bf6245c2384f915b1f6543f268dbd2 100644 (file)
@@ -7,7 +7,7 @@
  *  of special functions or implement the interface to the bignum package. */
 
 /*
- *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2016 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
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
+#ifdef HAVE_CONFIG_H
 #include "config.h"
-
-#include <vector>
-#include <stdexcept>
-#include <string>
-#include <sstream>
-#include <limits>
+#endif
 
 #include "numeric.h"
 #include "ex.h"
 #include "operators.h"
 #include "archive.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
@@ -88,7 +89,7 @@ numeric::numeric()
 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
+       // first.  This is due to the behavior of the cl_I-ctor, which
        // emphasizes efficiency.  However, if the integer is small enough
        // we save space and dereferences by using an immediate type.
        // (C.f. <cln/object.h>)
@@ -109,7 +110,7 @@ numeric::numeric(int i)
 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
+       // first.  This is due to the behavior of the cl_I-ctor, which
        // emphasizes efficiency.  However, if the integer is small enough
        // we save space and dereferences by using an immediate type.
        // (C.f. <cln/object.h>)
@@ -223,7 +224,7 @@ numeric::numeric(const char *s)
                        // E to lower case
                        term = term.replace(term.find("E"),1,"e");
                        // append _<Digits> to term
-                       term += "_" + ToString((unsigned)Digits);
+                       term += "_" + std::to_string((unsigned)Digits);
                        // construct float using cln::cl_F(const char *) ctor.
                        if (imaginary)
                                ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_F(term.c_str()));
@@ -278,8 +279,9 @@ static const cln::cl_F read_real_float(std::istream& s)
        return x;
 }
 
-numeric::numeric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+void numeric::read_archive(const archive_node &n, lst &sym_lst)
 {
+       inherited::read_archive(n, sym_lst);
        value = 0;
        
        // Read number as string
@@ -324,6 +326,7 @@ numeric::numeric(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
        }
        setflag(status_flags::evaluated | status_flags::expanded);
 }
+GINAC_BIND_UNARCHIVER(numeric);
 
 static void write_real_float(std::ostream& s, const cln::cl_R& n)
 {
@@ -373,8 +376,6 @@ void numeric::archive(archive_node &n) const
        n.add_string("number", s.str());
 }
 
-DEFAULT_UNARCHIVE(numeric)
-
 //////////
 // functions overriding virtual functions from base classes
 //////////
@@ -441,7 +442,7 @@ static void print_real_csrc(const print_context & c, const cln::cl_R & x)
                // 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) {
+               if (cln::plusp(x)) {
                        c.s << "(";
                        print_integer_csrc(c, numer);
                } else {
@@ -500,20 +501,20 @@ static void print_real_cl_N(const print_context & c, const cln::cl_R & x)
 {
        if (cln::instanceof(x, cln::cl_I_ring)) {
 
-    int dst;
-    // fixnum 
-    if (coerce(dst, cln::the<cln::cl_I>(x))) {
-      // can be converted to native int
-      if (dst < 0)
-        c.s << "(-" << dst << ")";
-      else
-        c.s << dst;
-    } else {
-      // bignum
-      c.s << "cln::cl_I(\"";
-      print_real_number(c, x);
-      c.s << "\")";
-    }
+               int dst;
+               // fixnum 
+               if (coerce(dst, cln::the<cln::cl_I>(x))) {
+                       // can be converted to native int
+                       if (dst < 0)
+                               c.s << "(-" << dst << ")";
+                       else
+                               c.s << dst;
+               } else {
+                       // bignum
+                       c.s << "cln::cl_I(\"";
+                       print_real_number(c, x);
+                       c.s << "\")";
+               }
        } else if (cln::instanceof(x, cln::cl_RA_ring)) {
 
                // Rational number
@@ -700,7 +701,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:
@@ -713,8 +714,6 @@ bool numeric::info(unsigned inf) const
                        return is_odd();
                case info_flags::prime:
                        return is_prime();
-               case info_flags::algebraic:
-                       return !is_real();
        }
        return false;
 }
@@ -774,10 +773,8 @@ bool numeric::has(const ex &other, unsigned options) const
 
 
 /** Evaluation of numbers doesn't do anything at all. */
-ex numeric::eval(int level) const
+ex numeric::eval() const
 {
-       // Warning: if this is ever gonna do something, the ex ctors from all kinds
-       // of numbers should be checking for status_flags::evaluated.
        return this->hold();
 }
 
@@ -787,11 +784,9 @@ ex numeric::eval(int level) const
  *  currently set.  In case the object already was a floating point number the
  *  precision is trimmed to match the currently set default.
  *
- *  @param level  ignored, only needed for overriding basic::evalf.
  *  @return  an ex-handle to a numeric. */
-ex numeric::evalf(int level) const
+ex numeric::evalf() const
 {
-       // level can safely be discarded for numeric objects.
        return numeric(cln::cl_float(1.0, cln::default_float_format) * value);
 }
 
@@ -928,9 +923,8 @@ const numeric &numeric::add_dyn(const numeric &other) const
                return other;
        else if (&other==_num0_p)
                return *this;
-       
-       return static_cast<const numeric &>((new numeric(value + other.value))->
-                                           setflag(status_flags::dynallocated));
+
+       return dynallocate<numeric>(value + other.value);
 }
 
 
@@ -944,9 +938,8 @@ const numeric &numeric::sub_dyn(const numeric &other) const
        // hack is supposed to keep the number of distinct numeric objects low.
        if (&other==_num0_p || cln::zerop(other.value))
                return *this;
-       
-       return static_cast<const numeric &>((new numeric(value - other.value))->
-                                           setflag(status_flags::dynallocated));
+
+       return dynallocate<numeric>(value - other.value);
 }
 
 
@@ -963,8 +956,7 @@ const numeric &numeric::mul_dyn(const numeric &other) const
        else if (&other==_num1_p)
                return *this;
        
-       return static_cast<const numeric &>((new numeric(value * other.value))->
-                                           setflag(status_flags::dynallocated));
+       return dynallocate<numeric>(value * other.value);
 }
 
 
@@ -982,8 +974,8 @@ const numeric &numeric::div_dyn(const numeric &other) const
                return *this;
        if (cln::zerop(cln::the<cln::cl_N>(other.value)))
                throw std::overflow_error("division by zero");
-       return static_cast<const numeric &>((new numeric(value / other.value))->
-                                           setflag(status_flags::dynallocated));
+
+       return dynallocate<numeric>(value / other.value);
 }
 
 
@@ -1009,8 +1001,8 @@ const numeric &numeric::power_dyn(const numeric &other) const
                else
                        return *_num0_p;
        }
-       return static_cast<const numeric &>((new numeric(cln::expt(value, other.value)))->
-                                            setflag(status_flags::dynallocated));
+
+       return dynallocate<numeric>(cln::expt(value, other.value));
 }
 
 
@@ -1374,7 +1366,7 @@ const numeric numeric::numer() const
                if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_RA_ring)) {
                        const cln::cl_I s = cln::lcm(cln::denominator(r), cln::denominator(i));
                        return numeric(cln::complex(cln::numerator(r)*(cln::exquo(s,cln::denominator(r))),
-                                                           cln::numerator(i)*(cln::exquo(s,cln::denominator(i)))));
+                                                   cln::numerator(i)*(cln::exquo(s,cln::denominator(i)))));
                }
        }
        // at least one float encountered
@@ -1529,7 +1521,7 @@ const numeric atan(const numeric &y, const numeric &x)
                return *_num0_p;
        if (x.is_real() && y.is_real())
                return numeric(cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
-                                cln::the<cln::cl_R>(y.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)))
@@ -1747,7 +1739,7 @@ class lanczos_coeffs
                std::vector<cln::cl_N> *current_vector;
 };
 
-std::vector<cln::cl_N>* lanczos_coeffs::coeffs = 0;
+std::vector<cln::cl_N>* lanczos_coeffs::coeffs = nullptr;
 
 bool lanczos_coeffs::sufficiently_accurate(int digits)
 {      if (digits<=20) {
@@ -1773,8 +1765,8 @@ cln::cl_N lanczos_coeffs::calc_lanczos_A(const cln::cl_N &x) const
 {
        cln::cl_N A = (*current_vector)[0];
        int size = current_vector->size();
-       for (int i=1; i<size; ++i)
-       A = A + (*current_vector)[i]/(x+cln::cl_I(-1+i));
+       for (int i=1; i<size; ++i)
+               A = A + (*current_vector)[i]/(x+cln::cl_I(-1+i));
        return A;
 }
 
@@ -2022,7 +2014,7 @@ lanczos_coeffs::lanczos_coeffs()
        coeffs[3].swap(coeffs_120);
 }
 
-static const cln::float_format_t guess_precision(const cln::cl_N& x)
+static 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))
@@ -2048,11 +2040,11 @@ const cln::cl_N lgamma(const cln::cl_N &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-cln::cl_N(1)/2)*log(temp)
-                             - temp
-                             + log(A);
-       return result;
+               cln::cl_N result = log(cln::cl_I(2)*pi_val)/2
+                                + (x-cln::cl_N(1)/2)*log(temp)
+                                - temp
+                                + log(A);
+               return result;
        }
        else 
                throw dunno();
@@ -2075,10 +2067,10 @@ const cln::cl_N tgamma(const cln::cl_N &x)
                        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 - cln::cl_N(1)/2)
-                         * exp(-temp) * A;
-       return result;
+               cln::cl_N result = sqrt(cln::cl_I(2)*pi_val)
+                                * expt(temp, x - cln::cl_N(1)/2)
+                                * exp(-temp) * A;
+               return result;
        }
        else
                throw dunno();
@@ -2225,7 +2217,7 @@ const numeric bernoulli(const numeric &nn)
 
        results.reserve(n/2);
        for (unsigned p=next_r; p<=n;  p+=2) {
-               cln::cl_I  c = 1;  // seed for binonmial coefficients
+               cln::cl_I  c = 1;  // seed for binomial coefficients
                cln::cl_RA b = cln::cl_RA(p-1)/-2;
                // The CLN manual says: "The conversion from `unsigned int' works only
                // if the argument is < 2^29" (This is for 32 Bit machines. More
@@ -2242,7 +2234,7 @@ const numeric bernoulli(const numeric &nn)
                                c = cln::exquo((c * (p+3-2*k)) * (p/2-k+1), cln::cl_I(2*k-1)*k);
                                b = b + c*results[k-1];
                        }
-               }
+               }
                results.push_back(-b/(p+1));
        }
        next_r = n+2;
@@ -2278,11 +2270,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);
@@ -2334,15 +2329,18 @@ const numeric mod(const numeric &a, const numeric &b)
 
 
 /** 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 numeric(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;
 }
@@ -2531,9 +2529,8 @@ _numeric_digits& _numeric_digits::operator=(long prec)
        cln::default_float_format = cln::float_format(prec);
 
        // call registered callbacks
-       std::vector<digits_changed_callback>::const_iterator it = callbacklist.begin(), end = callbacklist.end();
-       for (; it != end; ++it) {
-               (*it)(digitsdiff);
+       for (auto it : callbacklist) {
+               (it)(digitsdiff);
        }
 
        return *this;