- more logic on the trigonometric function stuff.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Mon, 20 Dec 1999 20:12:36 +0000 (20:12 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Mon, 20 Dec 1999 20:12:36 +0000 (20:12 +0000)
- changed several occurences of numeric const & to const numeric &, which
  is the same, but doxygen kept being confused because declarations
  differed from implementations.

17 files changed:
ginac/add.cpp
ginac/add.h
ginac/basic.h
ginac/ex.cpp
ginac/ex.h
ginac/inifcns.cpp
ginac/inifcns_gamma.cpp
ginac/inifcns_trans.cpp
ginac/inifcns_zeta.cpp
ginac/mul.cpp
ginac/mul.h
ginac/numeric.cpp
ginac/numeric.h
ginac/operators.h
ginac/power.cpp
ginac/utils.cpp
ginac/utils.h

index ad0ac1f..d376b5d 100644 (file)
@@ -160,7 +160,7 @@ basic * add::duplicate() const
     return new add(*this);
 }
 
-void add::print(ostream & os, unsigned upper_precedence) const
+/*void add::print(ostream & os, unsigned upper_precedence) const
 {
     debugmsg("add print",LOGLEVEL_PRINT);
     if (precedence<=upper_precedence) os << "(";
@@ -176,10 +176,17 @@ void add::print(ostream & os, unsigned upper_precedence) const
         }
         if (!coeff.is_equal(_num1()) &&
             !coeff.is_equal(_num_1())) {
-            if (coeff.csgn()==-1)
-                (_num_1()*coeff).print(os, precedence);
-            else
-                coeff.print(os, precedence);
+            if (coeff.is_rational()) {
+                if (coeff.is_negative())
+                    os << -coeff;
+                else
+                    os << coeff;
+            } else {
+                if (coeff.csgn()==-1)
+                    (-coeff).print(os, precedence);
+                else
+                    coeff.print(os, precedence);
+            }
             os << '*';
         }
         os << cit->rest;
@@ -190,6 +197,46 @@ void add::print(ostream & os, unsigned upper_precedence) const
         os << overall_coeff;
     }
     if (precedence<=upper_precedence) os << ")";
+}*/
+
+void add::print(ostream & os, unsigned upper_precedence) const
+{
+    debugmsg("add print",LOGLEVEL_PRINT);
+    if (precedence<=upper_precedence) os << "(";
+    numeric coeff;
+    bool first = true;
+    // first print the overall numeric coefficient, if present:
+    if (!overall_coeff.is_zero()) {
+        os << overall_coeff;
+        first = false;
+    }
+    // then proceed with the remaining factors:
+    for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
+        coeff = ex_to_numeric(cit->coeff);
+        if (!first) {
+            if (coeff.csgn()==-1) os << '-'; else os << '+';
+        } else {
+            if (coeff.csgn()==-1) os << '-';
+            first = false;
+        }
+        if (!coeff.is_equal(_num1()) &&
+            !coeff.is_equal(_num_1())) {
+            if (coeff.is_rational()) {
+                if (coeff.is_negative())
+                    os << -coeff;
+                else
+                    os << coeff;
+            } else {
+                if (coeff.csgn()==-1)
+                    (-coeff).print(os, precedence);
+                else
+                    coeff.print(os, precedence);
+            }
+            os << '*';
+        }
+        os << cit->rest;
+    }
+    if (precedence<=upper_precedence) os << ")";
 }
 
 void add::printraw(ostream & os) const
index 1943805..40cc25a 100644 (file)
@@ -72,7 +72,7 @@ public:
     ex series(symbol const & s, ex const & point, int order) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
     numeric integer_content(void) const;
-    ex smod(numeric const &xi) const;
+    ex smod(const numeric &xi) const;
     numeric max_coefficient(void) const;
     exvector get_indices(void) const;
     ex simplify_ncmul(exvector const & v) const;
index 87dc1de..d087186 100644 (file)
@@ -138,7 +138,7 @@ public: // only const functions please (may break reference counting)
     virtual ex subs(lst const & ls, lst const & lr) const;
     virtual ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
     virtual numeric integer_content(void) const;
-    virtual ex smod(numeric const &xi) const;
+    virtual ex smod(const numeric &xi) const;
     virtual numeric max_coefficient(void) const;
     virtual exvector get_indices(void) const;
     virtual ex simplify_ncmul(exvector const & v) const;
index 6a7a0aa..b853ced 100644 (file)
@@ -105,7 +105,22 @@ ex::ex(basic const & other)
 ex::ex(int const i)
 {
     debugmsg("ex constructor from int",LOGLEVEL_CONSTRUCT);
-    construct_from_basic(numeric(i));
+    switch (i) {  // some tiny efficiency-hack (FIXME: is this ok?)
+        case -1:
+            bp = _ex_1().bp;
+            ++bp->refcount;
+            break;
+        case 0:
+            bp = _ex0().bp;
+            ++bp->refcount;
+            break;
+        case 1:
+            bp = _ex1().bp;
+            ++bp->refcount;
+            break;
+        default:
+            construct_from_basic(numeric(i));
+    }
 }
 
 ex::ex(unsigned int const i)
index f8d2970..a719103 100644 (file)
@@ -44,7 +44,7 @@ extern ex const & _ex0(void);  /* FIXME: should this pollute headers? */
 
 #define INLINE_EX_CONSTRUCTORS
 
-/** Lightweight interface to GiNaC's symbolic objects. Basically all it does is
+/** Lightweight wrapper for GiNaC's symbolic objects.  Basically all it does is
  *  to hold a pointer to the other objects, manage the reference counting and
  *  provide methods for manipulation of these objects. */
 class ex
index 0a38c4b..28cadfb 100644 (file)
@@ -50,9 +50,9 @@ static ex Li2_eval(ex const & x)
     if (x.is_zero())
         return x;
     if (x.is_equal(_ex1()))
-        return power(Pi, 2) / 6;
+        return power(Pi, _ex2()) / _ex6();
     if (x.is_equal(_ex_1()))
-        return -power(Pi, 2) / 12;
+        return -power(Pi, _ex2()) / _ex12();
     return Li2(x).hold();
 }
 
@@ -142,7 +142,10 @@ static ex Order_series(ex const & x, symbol const & s, ex const & point, int ord
 
 REGISTER_FUNCTION(Order, Order_eval, NULL, NULL, Order_series);
 
-/** linear solve. */
+//////////
+// Solve linear system
+//////////
+
 ex lsolve(ex const &eqns, ex const &symbols)
 {
     // solve a system of linear equations
@@ -180,7 +183,7 @@ ex lsolve(ex const &eqns, ex const &symbols)
     matrix sys(eqns.nops(),symbols.nops());
     matrix rhs(eqns.nops(),1);
     matrix vars(symbols.nops(),1);
-
+    
     for (int r=0; r<eqns.nops(); r++) {
         ex eq=eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
         ex linpart=eq;
index 3f7fc67..dbe8840 100644 (file)
@@ -42,7 +42,7 @@ namespace GiNaC {
 // Gamma-function
 //////////
 
-static ex gamma_evalf(ex const & x)
+static ex gamma_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -56,7 +56,7 @@ static ex gamma_evalf(ex const & x)
  *  evaluation some day...
  *
  *  @exception std::domain_error("gamma_eval(): simple pole") */
-static ex gamma_eval(ex const & x)
+static ex gamma_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         // trap integer arguments:
@@ -72,25 +72,25 @@ static ex gamma_eval(ex const & x)
         if ((x*2).info(info_flags::integer)) {
             // trap positive x==(n+1/2)
             // gamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n)
-            if ((x*2).info(info_flags::posint)) {
+            if ((x*_ex2()).info(info_flags::posint)) {
                 numeric n = ex_to_numeric(x).sub(_num1_2());
                 numeric coefficient = doublefactorial(n.mul(_num2()).sub(_num1()));
-                coefficient = coefficient.div(_num2().power(n));
-                return coefficient * pow(Pi,_num1_2());
+                coefficient = coefficient.div(pow(_num2(),n));
+                return coefficient * pow(Pi,_ex1_2());
             } else {
                 // trap negative x==(-n+1/2)
                 // gamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1))
                 numeric n = abs(ex_to_numeric(x).sub(_num1_2()));
-                numeric coefficient = numeric(-2).power(n);
+                numeric coefficient = pow(_num_2(), n);
                 coefficient = coefficient.div(doublefactorial(n.mul(_num2()).sub(_num1())));;
-                return coefficient*sqrt(Pi);
+                return coefficient*power(Pi,_ex1_2());
             }
         }
     }
     return gamma(x).hold();
 }    
 
-static ex gamma_diff(ex const & x, unsigned diff_param)
+static ex gamma_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -99,10 +99,11 @@ static ex gamma_diff(ex const & x, unsigned diff_param)
     return psi(x)*gamma(x);
 }
 
-static ex gamma_series(ex const & x, symbol const & s, ex const & point, int order)
+static ex gamma_series(const ex & x, const symbol & s, const ex & point, int order)
 {
     // method:
-    // Taylor series where there is no pole falls back to psi function evaluation.
+    // Taylor series where there is no pole falls back to psi function
+    // evaluation.
     // On a pole at -m use the recurrence relation
     //   gamma(x) == gamma(x+1) / x
     // from which follows
@@ -126,7 +127,7 @@ REGISTER_FUNCTION(gamma, gamma_eval, gamma_evalf, gamma_diff, gamma_series);
 // Beta-function
 //////////
 
-static ex beta_evalf(ex const & x, ex const & y)
+static ex beta_evalf(const ex & x, const ex & y)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -137,24 +138,25 @@ static ex beta_evalf(ex const & x, ex const & y)
         / gamma(ex_to_numeric(x+y));
 }
 
-static ex beta_eval(ex const & x, ex const & y)
+static ex beta_eval(const ex & x, const ex & y)
 {
     if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) {
         numeric nx(ex_to_numeric(x));
         numeric ny(ex_to_numeric(y));
         // treat all problematic x and y that may not be passed into gamma,
-        // because they would throw there although beta(x,y) is well-defined:
+        // because they would throw there although beta(x,y) is well-defined
+        // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y)
         if (nx.is_real() && nx.is_integer() &&
             ny.is_real() && ny.is_integer()) {
             if (nx.is_negative()) {
                 if (nx<=-ny)
-                    return _num_1().power(ny)*beta(1-x-y, y);
+                    return pow(_num_1(), ny)*beta(1-x-y, y);
                 else
                     throw (std::domain_error("beta_eval(): simple pole"));
             }
             if (ny.is_negative()) {
                 if (ny<=-nx)
-                    return _num_1().power(nx)*beta(1-y-x, x);
+                    return pow(_num_1(), nx)*beta(1-y-x, x);
                 else
                     throw (std::domain_error("beta_eval(): simple pole"));
             }
@@ -164,13 +166,15 @@ static ex beta_eval(ex const & x, ex const & y)
         if ((nx+ny).is_real() &&
             (nx+ny).is_integer() &&
             !(nx+ny).is_positive())
-            return _ex0();
+             return _ex0();
+        // everything is ok:
         return gamma(x)*gamma(y)/gamma(x+y);
     }
+    
     return beta(x,y).hold();
 }
 
-static ex beta_diff(ex const & x, ex const & y, unsigned diff_param)
+static ex beta_diff(const ex & x, const ex & y, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param<2);
     ex retval;
@@ -179,18 +183,43 @@ static ex beta_diff(ex const & x, ex const & y, unsigned diff_param)
     if (diff_param==0)
         retval = (psi(x)-psi(x+y))*beta(x,y);
     // d/dy beta(x,y) -> (psi(y)-psi(x+y)) * beta(x,y)
-    if (diff_param==1)  
+    if (diff_param==1)
         retval = (psi(y)-psi(x+y))*beta(x,y);
     return retval;
 }
 
-REGISTER_FUNCTION(beta, beta_eval, beta_evalf, beta_diff, NULL);
+static ex beta_series(const ex & x, const ex & y, const symbol & s, const ex & point, int order)
+{
+    // method:
+    // Taylor series where there is no pole falls back to beta function
+    // evaluation.
+    // On a pole at -m use the recurrence relation
+    //   gamma(x) == gamma(x+1) / x
+    // from which follows
+    //   series(gamma(x),x,-m,order) ==
+    //   series(gamma(x+m+1)/(x*(x+1)...*(x+m)),x,-m,order+1);
+    ex xpoint = x.subs(s==point);
+    ex ypoint = y.subs(s==point);
+    if ((!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive)) &&
+        (!ypoint.info(info_flags::integer) || ypoint.info(info_flags::positive)))
+        throw do_taylor();  // caught by function::series()
+    // if we got here we have to care for a simple pole at -m:
+    throw (std::domain_error("beta_series(): please code me"));
+    /*numeric m = -ex_to_numeric(xpoint);
+     *ex ser_numer = gamma(x+m+_ex1());
+     *ex ser_denom = _ex1();
+     *for (numeric p; p<=m; ++p)
+     *    ser_denom *= x+p;
+     *return (ser_numer/ser_denom).series(s, point, order+1);*/
+}
+
+REGISTER_FUNCTION(beta, beta_eval, beta_evalf, beta_diff, beta_series);
 
 //////////
 // Psi-function (aka digamma-function)
 //////////
 
-static ex psi1_evalf(ex const & x)
+static ex psi1_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -201,7 +230,7 @@ static ex psi1_evalf(ex const & x)
 
 /** Evaluation of digamma-function psi(x).
  *  Somebody ought to provide some good numerical evaluation some day... */
-static ex psi1_eval(ex const & x)
+static ex psi1_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         if (x.info(info_flags::integer) && !x.info(info_flags::positive))
@@ -229,7 +258,7 @@ static ex psi1_eval(ex const & x)
     return psi(x).hold();
 }
 
-static ex psi1_diff(ex const & x, unsigned diff_param)
+static ex psi1_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -237,7 +266,7 @@ static ex psi1_diff(ex const & x, unsigned diff_param)
     return psi(_ex1(), x);
 }
 
-static ex psi1_series(ex const & x, symbol const & s, ex const & point, int order)
+static ex psi1_series(const ex & x, const symbol & s, const ex & point, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
@@ -264,7 +293,7 @@ const unsigned function_index_psi1 = function::register_new("psi", psi1_eval, ps
 // Psi-functions (aka polygamma-functions)  psi(0,x)==psi(x)
 //////////
 
-static ex psi2_evalf(ex const & n, ex const & x)
+static ex psi2_evalf(const ex & n, const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(n,numeric)
@@ -276,7 +305,7 @@ static ex psi2_evalf(ex const & n, ex const & x)
 
 /** Evaluation of polygamma-function psi(n,x). 
  *  Somebody ought to provide some good numerical evaluation some day... */
-static ex psi2_eval(ex const & n, ex const & x)
+static ex psi2_eval(const ex & n, const ex & x)
 {
     // psi(0,x) -> psi(x)
     if (n.is_zero())
@@ -288,13 +317,29 @@ static ex psi2_eval(ex const & n, ex const & x)
         x.info(info_flags::numeric)) {
         numeric nn = ex_to_numeric(n);
         numeric nx = ex_to_numeric(x);
-        if (x.is_equal(_ex1()))
-            return _num_1().power(nn+_num1())*factorial(nn)*zeta(ex(nn+_num1()));
+        if (nx.is_integer()) {
+            if (nx.is_equal(_num1()))
+                return pow(_num_1(), nn+_num1())*factorial(nn)*zeta(ex(nn+_num1()));
+            if (nx.is_positive()) {
+                // use the recurrence relation
+                //   psi(n,m) == psi(n,m+1) - (-)^n * n! / m^(n+1)
+                // to relate psi(n,m) to psi(n,1):
+                //   psi(n,m) == psi(n,1) + r
+                // where r == (-)^n * n! * (1^(-n-1) + ... + (m-1)^(-n-1))
+                numeric recur;
+                for (numeric p(1); p<nx; ++p)
+                    recur += pow(p, -nn+_num_1());
+                recur *= factorial(nn)*pow(_num_1(), nn);
+                return recur+psi(n,_ex1());
+            }
+            // for non-positive integers there is a pole:
+            throw (std::domain_error("psi2_eval(): pole"));
+        }
     }
     return psi(n, x).hold();
 }    
 
-static ex psi2_diff(ex const & n, ex const & x, unsigned diff_param)
+static ex psi2_diff(const ex & n, const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param<2);
     
@@ -303,20 +348,20 @@ static ex psi2_diff(ex const & n, ex const & x, unsigned diff_param)
         throw(std::logic_error("cannot diff psi(n,x) with respect to n"));
     }
     // d/dx psi(n,x) -> psi(n+1,x)
-    return psi(n+1, x);
+    return psi(n+_ex1(), x);
 }
 
-static ex psi2_series(ex const & n, ex const & x, symbol const & s, ex const & point, int order)
+static ex psi2_series(const ex & n, const ex & x, const symbol & s, const ex & point, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
     // evaluation.
     // On a pole at -m use the recurrence relation
-    //   psi(n,x) == psi(n,x+1) - (-)^n * n! / z^(n+1)
+    //   psi(n,x) == psi(n,x+1) - (-)^n * n! / x^(n+1)
     // from which follows
     //   series(psi(x),x,-m,order) == 
-    //   series(psi(x+m+1) - (-1)^n * n!
-    //            * ((x)^(-n-1) + (x+1)^(-n-1) + (x+m)^(-n-1))),x,-m,order);
+    //   series(psi(x+m+1) - (-1)^n * n! * ((x)^(-n-1) + (x+1)^(-n-1) + ...
+    //                                      ... + (x+m)^(-n-1))),x,-m,order);
     ex xpoint = x.subs(s==point);
     if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()
index 4473fe0..9937e30 100644 (file)
@@ -122,7 +122,7 @@ static ex log_eval(ex const & x)
     }
     // log(exp(t)) -> t (for real-valued t):
     if (is_ex_the_function(x, exp)) {
-        ex t=x.op(0);
+        ex t = x.op(0);
         if (t.info(info_flags::real))
             return t;
     }
@@ -135,7 +135,7 @@ static ex log_diff(ex const & x, unsigned diff_param)
     GINAC_ASSERT(diff_param==0);
 
     // d/dx log(x) -> 1/x
-    return power(x, -1);
+    return power(x, _ex_1());
 }
 
 REGISTER_FUNCTION(log, log_eval, log_evalf, log_diff, NULL);
@@ -155,38 +155,42 @@ static ex sin_evalf(ex const & x)
 
 static ex sin_eval(ex const & x)
 {
-    // sin(n*Pi/6) -> { 0 | +/-1/2 | +/-sqrt(3)/2 | +/-1 }
-    ex SixExOverPi = _ex6()*x/Pi;
-    if (SixExOverPi.info(info_flags::integer)) {
-        numeric z = smod(ex_to_numeric(SixExOverPi),_num12());
-        if (z.is_equal(_num_5()))  // sin(7*Pi/6)  -> -1/2
-            return _ex_1_2();
-        if (z.is_equal(_num_4()))  // sin(8*Pi/6)  -> -sqrt(3)/2
-            return _ex_1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num_3()))  // sin(9*Pi/6)  -> -1
-            return _ex_1();
-        if (z.is_equal(_num_2()))  // sin(10*Pi/6) -> -sqrt(3)/2
-            return _ex_1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num_1()))  // sin(11*Pi/6) -> -1/2
-            return _ex_1_2();
-        if (z.is_equal(_num0()))   // sin(0)       ->  0
-            return _ex0();
-        if (z.is_equal(_num1()))   // sin(1*Pi/6)  ->  1/2
-            return _ex1_2();
-        if (z.is_equal(_num2()))   // sin(2*Pi/6)  ->  sqrt(3)/2
-            return _ex1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num3()))   // sin(3*Pi/6)  ->  1
-            return _ex1();
-        if (z.is_equal(_num4()))   // sin(4*Pi/6)  ->  sqrt(3)/2
-            return _ex1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num5()))   // sin(5*Pi/6)  ->  1/2
-            return _ex1_2();
-        if (z.is_equal(_num6()))   // sin(6*Pi/6)  ->  0
+    // sin(n/d*Pi) -> { all known non-nested radicals }
+    ex SixtyExOverPi = _ex60()*x/Pi;
+    ex sign = _ex1();
+    if (SixtyExOverPi.info(info_flags::integer)) {
+        numeric z = mod(ex_to_numeric(SixtyExOverPi),_num120());
+        if (z>=_num60()) {
+            // wrap to interval [0, Pi)
+            z -= _num60();
+            sign = _ex_1();
+        }
+        if (z>_num30()) {
+            // wrap to interval [0, Pi/2)
+            z = _num60()-z;
+        }
+        if (z.is_equal(_num0()))  // sin(0)       -> 0
             return _ex0();
+        if (z.is_equal(_num5()))  // sin(Pi/12)   -> sqrt(6)/4*(1-sqrt(3)/3)
+            return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
+        if (z.is_equal(_num6()))  // sin(Pi/10)   -> sqrt(5)/4-1/4
+            return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
+        if (z.is_equal(_num10())) // sin(Pi/6)    -> 1/2
+            return sign*_ex1_2();
+        if (z.is_equal(_num15())) // sin(Pi/4)    -> sqrt(2)/2
+            return sign*_ex1_2()*power(_ex2(),_ex1_2());
+        if (z.is_equal(_num18())) // sin(3/10*Pi) -> sqrt(5)/4+1/4
+            return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
+        if (z.is_equal(_num20())) // sin(Pi/3)    -> sqrt(3)/2
+            return sign*_ex1_2()*power(_ex3(),_ex1_2());
+        if (z.is_equal(_num25())) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
+            return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
+        if (z.is_equal(_num30())) // sin(Pi/2)    -> 1
+            return sign*_ex1();
     }
     
     if (is_ex_exactly_of_type(x, function)) {
-        ex t=x.op(0);
+        ex t = x.op(0);
         // sin(asin(x)) -> x
         if (is_ex_the_function(x, asin))
             return t;
@@ -230,38 +234,42 @@ static ex cos_evalf(ex const & x)
 
 static ex cos_eval(ex const & x)
 {
-    // cos(n*Pi/6) -> { 0 | +/-1/2 | +/-sqrt(3)/2 | +/-1 }
-    ex SixExOverPi = _ex6()*x/Pi;
-    if (SixExOverPi.info(info_flags::integer)) {
-        numeric z = smod(ex_to_numeric(SixExOverPi),_num12());
-        if (z.is_equal(_num_5()))  // cos(7*Pi/6)  -> -sqrt(3)/2
-            return _ex_1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num_4()))  // cos(8*Pi/6)  -> -1/2
-            return _ex_1_2();
-        if (z.is_equal(_num_3()))  // cos(9*Pi/6)  ->  0
-            return _ex0();
-        if (z.is_equal(_num_2()))  // cos(10*Pi/6) ->  1/2
-            return _ex1_2();
-        if (z.is_equal(_num_1()))  // cos(11*Pi/6) ->  sqrt(3)/2
-            return _ex1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num0()))   // cos(0)       ->  1
-            return _ex1();
-        if (z.is_equal(_num1()))   // cos(1*Pi/6)  ->  sqrt(3)/2
-            return _ex1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num2()))   // cos(2*Pi/6)  ->  1/2
-            return _ex1_2();
-        if (z.is_equal(_num3()))   // cos(3*Pi/6)  ->  0
-            return _ex0();
-        if (z.is_equal(_num4()))   // cos(4*Pi/6)  -> -1/2
-            return _ex_1_2();
-        if (z.is_equal(_num5()))   // cos(5*Pi/6)  -> -sqrt(3)/2
-            return _ex_1_2()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num6()))   // cos(6*Pi/6)  -> -1
-            return _ex_1();
+    // cos(n/d*Pi) -> { all known non-nested radicals }
+    ex SixtyExOverPi = _ex60()*x/Pi;
+    ex sign = _ex1();
+    if (SixtyExOverPi.info(info_flags::integer)) {
+        numeric z = mod(ex_to_numeric(SixtyExOverPi),_num120());
+        if (z>=_num60()) {
+            // wrap to interval [0, Pi)
+            z = _num120()-z;
+        }
+        if (z>=_num30()) {
+            // wrap to interval [0, Pi/2)
+            z = _num60()-z;
+            sign = _ex_1();
+        }
+        if (z.is_equal(_num0()))  // cos(0)       -> 1
+            return sign*_ex1();
+        if (z.is_equal(_num5()))  // cos(Pi/12)   -> sqrt(6)/4*(1+sqrt(3)/3)
+            return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
+        if (z.is_equal(_num10())) // cos(Pi/6)    -> sqrt(3)/2
+            return sign*_ex1_2()*power(_ex3(),_ex1_2());
+        if (z.is_equal(_num12())) // cos(Pi/5)    -> sqrt(5)/4+1/4
+            return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
+        if (z.is_equal(_num15())) // cos(Pi/4)    -> sqrt(2)/2
+            return sign*_ex1_2()*power(_ex2(),_ex1_2());
+        if (z.is_equal(_num20())) // cos(Pi/3)    -> 1/2
+            return sign*_ex1_2();
+        if (z.is_equal(_num24())) // cos(2/5*Pi)  -> sqrt(5)/4-1/4x
+            return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
+        if (z.is_equal(_num25())) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
+            return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
+        if (z.is_equal(_num30())) // cos(Pi/2)    -> 0
+            return sign*_ex0();
     }
     
     if (is_ex_exactly_of_type(x, function)) {
-        ex t=x.op(0);
+        ex t = x.op(0);
         // cos(acos(x)) -> x
         if (is_ex_the_function(x, acos))
             return t;
@@ -305,26 +313,38 @@ static ex tan_evalf(ex const & x)
 
 static ex tan_eval(ex const & x)
 {
-    // tan(n*Pi/6) -> { 0 | +/-sqrt(3) | +/-sqrt(3)/2 }
-    ex SixExOverPi = _ex6()*x/Pi;
-    if (SixExOverPi.info(info_flags::integer)) {
-        numeric z = smod(ex_to_numeric(SixExOverPi),_num6());
-        if (z.is_equal(_num_2()))  // tan(4*Pi/6) -> -sqrt(3)
-            return _ex_1()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num_1()))  // tan(5*Pi/6) -> -sqrt(3)/3
-            return _ex_1_3()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num0()))   // tan(0)      ->  0
+    // tan(n/d*Pi) -> { all known non-nested radicals }
+    ex SixtyExOverPi = _ex60()*x/Pi;
+    ex sign = _ex1();
+    if (SixtyExOverPi.info(info_flags::integer)) {
+        numeric z = mod(ex_to_numeric(SixtyExOverPi),_num60());
+        if (z>=_num60()) {
+            // wrap to interval [0, Pi)
+            z -= _num60();
+        }
+        if (z>=_num30()) {
+            // wrap to interval [0, Pi/2)
+            z = _num60()-z;
+            sign = _ex_1();
+        }
+        if (z.is_equal(_num0()))  // tan(0)       -> 0
             return _ex0();
-        if (z.is_equal(_num1()))   // tan(1*Pi/6) ->  sqrt(3)/3
-            return _ex1_3()*power(_ex3(),_ex1_2());
-        if (z.is_equal(_num2()))   // tan(2*Pi/6) ->  sqrt(3)
-            return power(_ex3(),_ex1_2());
-        if (z.is_equal(_num3()))   // tan(3*Pi/6) ->  infinity
+        if (z.is_equal(_num5()))  // tan(Pi/12)   -> 2-sqrt(3)
+            return sign*(_ex2()-power(_ex3(),_ex1_2()));
+        if (z.is_equal(_num10())) // tan(Pi/6)    -> sqrt(3)/3
+            return sign*_ex1_3()*power(_ex3(),_ex1_2());
+        if (z.is_equal(_num15())) // tan(Pi/4)    -> 1
+            return sign*_ex1();
+        if (z.is_equal(_num20())) // tan(Pi/3)    -> sqrt(3)
+            return sign*power(_ex3(),_ex1_2());
+        if (z.is_equal(_num25())) // tan(5/12*Pi) -> 2+sqrt(3)
+            return sign*(power(_ex3(),_ex1_2())+_ex2());
+        if (z.is_equal(_num30())) // tan(Pi/2)    -> infinity
             throw (std::domain_error("tan_eval(): infinity"));
     }
-        
+    
     if (is_ex_exactly_of_type(x, function)) {
-        ex t=x.op(0);
+        ex t = x.op(0);
         // tan(atan(x)) -> x
         if (is_ex_the_function(x, atan))
             return t;
@@ -349,7 +369,7 @@ static ex tan_diff(ex const & x, unsigned diff_param)
     GINAC_ASSERT(diff_param==0);
     
     // d/dx tan(x) -> 1+tan(x)^2;
-    return (1+power(tan(x),_ex2()));
+    return (_ex1()+power(tan(x),_ex2()));
 }
 
 static ex tan_series(ex const & x, symbol const & s, ex const & point, int order)
@@ -436,10 +456,10 @@ static ex acos_eval(ex const & x)
             return _ex0();
         // acos(1/2) -> Pi/3
         if (x.is_equal(_ex1_2()))
-            return numeric(1,3)*Pi;
+            return _ex1_3()*Pi;
         // acos(0) -> Pi/2
         if (x.is_zero())
-            return numeric(1,2)*Pi;
+            return _ex1_2()*Pi;
         // acos(-1/2) -> 2/3*Pi
         if (x.is_equal(_ex_1_2()))
             return numeric(2,3)*Pi;
@@ -495,7 +515,8 @@ static ex atan_diff(ex const & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
 
-    return power(1+x*x, -1);
+    // d/dx atan(x) -> 1/(1+x^2)
+    return power(_ex1()+power(x,_ex2()), _ex_1());
 }
 
 REGISTER_FUNCTION(atan, atan_eval, atan_evalf, atan_diff, NULL);
@@ -530,10 +551,10 @@ static ex atan2_diff(ex const & y, ex const & x, unsigned diff_param)
     
     if (diff_param==0) {
         // d/dy atan(y,x)
-        return x*pow(pow(x,2)+pow(y,2),-1);
+        return x*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
     }
     // d/dx atan(y,x)
-    return -y*pow(pow(x,2)+pow(y,2),-1);
+    return -y*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
 }
 
 REGISTER_FUNCTION(atan2, atan2_eval, atan2_evalf, atan2_diff, NULL);
@@ -565,7 +586,7 @@ static ex sinh_eval(ex const & x)
         return I*sin(x/I);
     
     if (is_ex_exactly_of_type(x, function)) {
-        ex t=x.op(0);
+        ex t = x.op(0);
         // sinh(asinh(x)) -> x
         if (is_ex_the_function(x, asinh))
             return t;
@@ -617,7 +638,7 @@ static ex cosh_eval(ex const & x)
         return cos(x/I);
     
     if (is_ex_exactly_of_type(x, function)) {
-        ex t=x.op(0);
+        ex t = x.op(0);
         // cosh(acosh(x)) -> x
         if (is_ex_the_function(x, acosh))
             return t;
@@ -669,7 +690,7 @@ static ex tanh_eval(ex const & x)
         return I*tan(x/I);
     
     if (is_ex_exactly_of_type(x, function)) {
-        ex t=x.op(0);
+        ex t = x.op(0);
         // tanh(atanh(x)) -> x
         if (is_ex_the_function(x, atanh))
             return t;
@@ -738,7 +759,7 @@ static ex asinh_diff(ex const & x, unsigned diff_param)
     GINAC_ASSERT(diff_param==0);
     
     // d/dx asinh(x) -> 1/sqrt(1+x^2)
-    return power(1+power(x,_ex2()),_ex_1_2());
+    return power(_ex1()+power(x,_ex2()),_ex_1_2());
 }
 
 REGISTER_FUNCTION(asinh, asinh_eval, asinh_evalf, asinh_diff, NULL);
index 9fdce1e..dfd4638 100644 (file)
@@ -62,7 +62,7 @@ static ex zeta1_eval(ex const & x)
                 if (x.info(info_flags::odd))
                     return zeta(x).hold();
                 else
-                    return abs(bernoulli(y))*pow(Pi,x)*_num2().power(y-_num1())/factorial(y);
+                    return abs(bernoulli(y))*pow(Pi,x)*pow(_num2(),y-_num1())/factorial(y);
             } else {
                 if (x.info(info_flags::odd))
                     return -bernoulli(_num1()-y)/(_num1()-y);
index 9bee46f..b444b1e 100644 (file)
@@ -181,13 +181,21 @@ void mul::print(ostream & os, unsigned upper_precedence) const
     if (precedence<=upper_precedence) os << "(";
     bool first=true;
     // first print the overall numeric coefficient:
-    if (ex_to_numeric(overall_coeff).csgn()==-1) os << '-';
-    if (!overall_coeff.is_equal(_ex1()) &&
-        !overall_coeff.is_equal(_ex_1())) {
-        if (ex_to_numeric(overall_coeff).csgn()==-1)
-            (_num_1()*overall_coeff).print(os, precedence);
-        else
-            overall_coeff.print(os, precedence);
+    numeric coeff = ex_to_numeric(overall_coeff);
+    if (coeff.csgn()==-1) os << '-';
+    if (!coeff.is_equal(_num1()) &&
+        !coeff.is_equal(_num_1())) {
+        if (coeff.is_rational()) {
+            if (coeff.is_negative())
+                os << -coeff;
+            else
+                os << coeff;
+        } else {
+            if (coeff.csgn()==-1)
+                (-coeff).print(os, precedence);
+            else
+                coeff.print(os, precedence);
+        }
         os << '*';
     }
     // then proceed with the remaining factors:
index 71b686b..1d95071 100644 (file)
@@ -73,7 +73,7 @@ public:
     ex series(symbol const & s, ex const & point, int order) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
     numeric integer_content(void) const;
-    ex smod(numeric const &xi) const;
+    ex smod(const numeric &xi) const;
     numeric max_coefficient(void) const;
     exvector get_indices(void) const;
     ex simplify_ncmul(exvector const & v) const;
index 472d9c3..4937448 100644 (file)
@@ -72,13 +72,13 @@ numeric::~numeric()
     destroy(0);
 }
 
-numeric::numeric(numeric const & other)
+numeric::numeric(const numeric & other)
 {
     debugmsg("numeric copy constructor", LOGLEVEL_CONSTRUCT);
     copy(other);
 }
 
-numeric const & numeric::operator=(numeric const & other)
+const numeric & numeric::operator=(const numeric & other)
 {
     debugmsg("numeric operator=", LOGLEVEL_ASSIGNMENT);
     if (this != &other) {
@@ -90,7 +90,7 @@ numeric const & numeric::operator=(numeric const & other)
 
 // protected
 
-void numeric::copy(numeric const & other)
+void numeric::copy(const numeric & other)
 {
     basic::copy(other);
     value = new cl_N(*other.value);
@@ -110,7 +110,7 @@ void numeric::destroy(bool call_parent)
 
 numeric::numeric(int i) : basic(TINFO_numeric)
 {
-    debugmsg("numeric constructor from int",LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from int",LOGLEVEL_CONSTRUCT);
     // 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:
@@ -122,7 +122,7 @@ numeric::numeric(int i) : basic(TINFO_numeric)
 
 numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 {
-    debugmsg("numeric constructor from uint",LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from uint",LOGLEVEL_CONSTRUCT);
     // 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:
@@ -134,7 +134,7 @@ numeric::numeric(unsigned int i) : basic(TINFO_numeric)
 
 numeric::numeric(long i) : basic(TINFO_numeric)
 {
-    debugmsg("numeric constructor from long",LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from long",LOGLEVEL_CONSTRUCT);
     value = new cl_I(i);
     calchash();
     setflag(status_flags::evaluated|
@@ -143,7 +143,7 @@ numeric::numeric(long i) : basic(TINFO_numeric)
 
 numeric::numeric(unsigned long i) : basic(TINFO_numeric)
 {
-    debugmsg("numeric constructor from ulong",LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from ulong",LOGLEVEL_CONSTRUCT);
     value = new cl_I(i);
     calchash();
     setflag(status_flags::evaluated|
@@ -155,7 +155,7 @@ numeric::numeric(unsigned long i) : basic(TINFO_numeric)
  *  @exception overflow_error (division by zero) */
 numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
 {
-    debugmsg("numeric constructor from long/long",LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from long/long",LOGLEVEL_CONSTRUCT);
     if (!denom)
         throw (std::overflow_error("division by zero"));
     value = new cl_I(numer);
@@ -167,7 +167,7 @@ numeric::numeric(long numer, long denom) : basic(TINFO_numeric)
 
 numeric::numeric(double d) : basic(TINFO_numeric)
 {
-    debugmsg("numeric constructor from double",LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from double",LOGLEVEL_CONSTRUCT);
     // 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
     // will not be promoted to cl_LF if overflow occurs:
@@ -180,7 +180,7 @@ numeric::numeric(double d) : basic(TINFO_numeric)
 
 numeric::numeric(char const *s) : basic(TINFO_numeric)
 {   // MISSING: treatment of complex and ints and rationals.
-    debugmsg("numeric constructor from string",LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from string",LOGLEVEL_CONSTRUCT);
     if (strchr(s, '.'))
         value = new cl_LF(s);
     else
@@ -194,7 +194,7 @@ numeric::numeric(char const *s) : basic(TINFO_numeric)
  *  only. */
 numeric::numeric(cl_N const & z) : basic(TINFO_numeric)
 {
-    debugmsg("numeric constructor from cl_N", LOGLEVEL_CONSTRUCT);
+    debugmsg("const numericructor from cl_N", LOGLEVEL_CONSTRUCT);
     value = new cl_N(z);
     calchash();
     setflag(status_flags::evaluated|
@@ -213,19 +213,11 @@ basic * numeric::duplicate() const
     return new numeric(*this);
 }
 
-// The method printraw doesn't do much, it simply uses CLN's operator<<() for
-// output, which is ugly but reliable. Examples:
-// 2+2i 
-void numeric::printraw(ostream & os) const
-{
-    debugmsg("numeric printraw", LOGLEVEL_PRINT);
-    os << "numeric(" << *value << ")";
-}
-
-// The method print adds to the output so it blends more consistently together
-// with the other routines and produces something compatible to Maple input.
 void numeric::print(ostream & os, unsigned upper_precedence) const
 {
+    // The method print adds to the output so it blends more consistently
+    // together with the other routines and produces something compatible to
+    // ginsh input.
     debugmsg("numeric print", LOGLEVEL_PRINT);
     if (is_real()) {
         // case 1, real:  x  or  -x
@@ -276,6 +268,14 @@ void numeric::print(ostream & os, unsigned upper_precedence) const
     }
 }
 
+
+void numeric::printraw(ostream & os) const
+{
+    // The method printraw doesn't do much, it simply uses CLN's operator<<()
+    // for output, which is ugly but reliable. e.g: 2+2i
+    debugmsg("numeric printraw", LOGLEVEL_PRINT);
+    os << "numeric(" << *value << ")";
+}
 void numeric::printtree(ostream & os, unsigned indent) const
 {
     debugmsg("numeric printtree", LOGLEVEL_PRINT);
@@ -379,7 +379,7 @@ ex numeric::evalf(int level) const
 int numeric::compare_same_type(basic const & other) const
 {
     GINAC_ASSERT(is_exactly_of_type(other, numeric));
-    numeric const & o = static_cast<numeric &>(const_cast<basic &>(other));
+    const numeric & o = static_cast<numeric &>(const_cast<basic &>(other));
 
     if (*value == *o.value) {
         return 0;
@@ -391,7 +391,7 @@ int numeric::compare_same_type(basic const & other) const
 bool numeric::is_equal_same_type(basic const & other) const
 {
     GINAC_ASSERT(is_exactly_of_type(other,numeric));
-    numeric const *o = static_cast<numeric const *>(&other);
+    const numeric *o = static_cast<const numeric *>(&other);
     
     return is_equal(*o);
 }
@@ -424,21 +424,21 @@ unsigned numeric::calchash(void) const
 
 /** Numerical addition method.  Adds argument to *this and returns result as
  *  a new numeric object. */
-numeric numeric::add(numeric const & other) const
+numeric numeric::add(const numeric & other) const
 {
     return numeric((*value)+(*other.value));
 }
 
 /** Numerical subtraction method.  Subtracts argument from *this and returns
  *  result as a new numeric object. */
-numeric numeric::sub(numeric const & other) const
+numeric numeric::sub(const numeric & other) const
 {
     return numeric((*value)-(*other.value));
 }
 
 /** Numerical multiplication method.  Multiplies *this and argument and returns
  *  result as a new numeric object. */
-numeric numeric::mul(numeric const & other) const
+numeric numeric::mul(const numeric & other) const
 {
     static const numeric * _num1p=&_num1();
     if (this==_num1p) {
@@ -453,14 +453,14 @@ numeric numeric::mul(numeric const & other) const
  *  a new numeric object.
  *
  *  @exception overflow_error (division by zero) */
-numeric numeric::div(numeric const & other) const
+numeric numeric::div(const numeric & other) const
 {
     if (::zerop(*other.value))
         throw (std::overflow_error("division by zero"));
     return numeric((*value)/(*other.value));
 }
 
-numeric numeric::power(numeric const & other) const
+numeric numeric::power(const numeric & other) const
 {
     static const numeric * _num1p=&_num1();
     if (&other==_num1p)
@@ -476,19 +476,19 @@ numeric numeric::inverse(void) const
     return numeric(::recip(*value));  // -> CLN
 }
 
-numeric const & numeric::add_dyn(numeric const & other) const
+const numeric & numeric::add_dyn(const numeric & other) const
 {
-    return static_cast<numeric const &>((new numeric((*value)+(*other.value)))->
+    return static_cast<const numeric &>((new numeric((*value)+(*other.value)))->
                                         setflag(status_flags::dynallocated));
 }
 
-numeric const & numeric::sub_dyn(numeric const & other) const
+const numeric & numeric::sub_dyn(const numeric & other) const
 {
-    return static_cast<numeric const &>((new numeric((*value)-(*other.value)))->
+    return static_cast<const numeric &>((new numeric((*value)-(*other.value)))->
                                         setflag(status_flags::dynallocated));
 }
 
-numeric const & numeric::mul_dyn(numeric const & other) const
+const numeric & numeric::mul_dyn(const numeric & other) const
 {
     static const numeric * _num1p=&_num1();
     if (this==_num1p) {
@@ -496,55 +496,55 @@ numeric const & numeric::mul_dyn(numeric const & other) const
     } else if (&other==_num1p) {
         return *this;
     }
-    return static_cast<numeric const &>((new numeric((*value)*(*other.value)))->
+    return static_cast<const numeric &>((new numeric((*value)*(*other.value)))->
                                         setflag(status_flags::dynallocated));
 }
 
-numeric const & numeric::div_dyn(numeric const & other) const
+const numeric & numeric::div_dyn(const numeric & other) const
 {
     if (::zerop(*other.value))
         throw (std::overflow_error("division by zero"));
-    return static_cast<numeric const &>((new numeric((*value)/(*other.value)))->
+    return static_cast<const numeric &>((new numeric((*value)/(*other.value)))->
                                         setflag(status_flags::dynallocated));
 }
 
-numeric const & numeric::power_dyn(numeric const & other) const
+const numeric & numeric::power_dyn(const numeric & other) const
 {
     static const numeric * _num1p=&_num1();
     if (&other==_num1p)
         return *this;
     if (::zerop(*value) && other.is_real() && ::minusp(realpart(*other.value)))
         throw (std::overflow_error("division by zero"));
-    return static_cast<numeric const &>((new numeric(::expt(*value,*other.value)))->
+    return static_cast<const numeric &>((new numeric(::expt(*value,*other.value)))->
                                         setflag(status_flags::dynallocated));
 }
 
-numeric const & numeric::operator=(int i)
+const numeric & numeric::operator=(int i)
 {
     return operator=(numeric(i));
 }
 
-numeric const & numeric::operator=(unsigned int i)
+const numeric & numeric::operator=(unsigned int i)
 {
     return operator=(numeric(i));
 }
 
-numeric const & numeric::operator=(long i)
+const numeric & numeric::operator=(long i)
 {
     return operator=(numeric(i));
 }
 
-numeric const & numeric::operator=(unsigned long i)
+const numeric & numeric::operator=(unsigned long i)
 {
     return operator=(numeric(i));
 }
 
-numeric const & numeric::operator=(double d)
+const numeric & numeric::operator=(double d)
 {
     return operator=(numeric(d));
 }
 
-numeric const & numeric::operator=(char const * s)
+const numeric & numeric::operator=(char const * s)
 {
     return operator=(numeric(s));
 }
@@ -553,7 +553,7 @@ numeric const & numeric::operator=(char const * s)
  *  csgn(x)==0 for x==0, csgn(x)==1 for Re(x)>0 or Re(x)=0 and Im(x)>0,
  *  csgn(x)==-1 for Re(x)<0 or Re(x)=0 and Im(x)<0.
  *
- *  @see numeric::compare(numeric const & other) */
+ *  @see numeric::compare(const numeric & other) */
 int numeric::csgn(void) const
 {
     if (is_zero())
@@ -578,7 +578,7 @@ int numeric::csgn(void) const
  *
  *  @return csgn(*this-other)
  *  @see numeric::csgn(void) */
-int numeric::compare(numeric const & other) const
+int numeric::compare(const numeric & other) const
 {
     // Comparing two real numbers?
     if (is_real() && other.is_real())
@@ -594,7 +594,7 @@ int numeric::compare(numeric const & other) const
     }
 }
 
-bool numeric::is_equal(numeric const & other) const
+bool numeric::is_equal(const numeric & other) const
 {
     return (*value == *other.value);
 }
@@ -672,12 +672,12 @@ bool numeric::is_real(void) const
     return ::instanceof(*value, cl_R_ring);  // -> CLN
 }
 
-bool numeric::operator==(numeric const & other) const
+bool numeric::operator==(const numeric & other) const
 {
     return (*value == *other.value);  // -> CLN
 }
 
-bool numeric::operator!=(numeric const & other) const
+bool numeric::operator!=(const numeric & other) const
 {
     return (*value != *other.value);  // -> CLN
 }
@@ -713,7 +713,7 @@ bool numeric::is_crational(void) const
 /** Numerical comparison: less.
  *
  *  @exception invalid_argument (complex inequality) */ 
-bool numeric::operator<(numeric const & other) const
+bool numeric::operator<(const numeric & other) const
 {
     if (is_real() && other.is_real())
         return (bool)(The(cl_R)(*value) < The(cl_R)(*other.value));  // -> CLN
@@ -724,7 +724,7 @@ bool numeric::operator<(numeric const & other) const
 /** Numerical comparison: less or equal.
  *
  *  @exception invalid_argument (complex inequality) */ 
-bool numeric::operator<=(numeric const & other) const
+bool numeric::operator<=(const numeric & other) const
 {
     if (is_real() && other.is_real())
         return (bool)(The(cl_R)(*value) <= The(cl_R)(*other.value));  // -> CLN
@@ -735,7 +735,7 @@ bool numeric::operator<=(numeric const & other) const
 /** Numerical comparison: greater.
  *
  *  @exception invalid_argument (complex inequality) */ 
-bool numeric::operator>(numeric const & other) const
+bool numeric::operator>(const numeric & other) const
 {
     if (is_real() && other.is_real())
         return (bool)(The(cl_R)(*value) > The(cl_R)(*other.value));  // -> CLN
@@ -746,7 +746,7 @@ bool numeric::operator>(numeric const & other) const
 /** Numerical comparison: greater or equal.
  *
  *  @exception invalid_argument (complex inequality) */  
-bool numeric::operator>=(numeric const & other) const
+bool numeric::operator>=(const numeric & other) const
 {
     if (is_real() && other.is_real())
         return (bool)(The(cl_R)(*value) >= The(cl_R)(*other.value));  // -> CLN
@@ -928,7 +928,7 @@ const numeric I = numeric(complex(cl_I(0),cl_I(1)));
 /** Exponential function.
  *
  *  @return  arbitrary precision numerical exp(x). */
-numeric exp(numeric const & x)
+numeric exp(const numeric & x)
 {
     return ::exp(*x.value);  // -> CLN
 }
@@ -938,7 +938,7 @@ numeric exp(numeric const & x)
  *  @param z complex number
  *  @return  arbitrary precision numerical log(x).
  *  @exception overflow_error (logarithmic singularity) */
-numeric log(numeric const & z)
+numeric log(const numeric & z)
 {
     if (z.is_zero())
         throw (std::overflow_error("log(): logarithmic singularity"));
@@ -948,7 +948,7 @@ numeric log(numeric const & z)
 /** Numeric sine (trigonometric function).
  *
  *  @return  arbitrary precision numerical sin(x). */
-numeric sin(numeric const & x)
+numeric sin(const numeric & x)
 {
     return ::sin(*x.value);  // -> CLN
 }
@@ -956,7 +956,7 @@ numeric sin(numeric const & x)
 /** Numeric cosine (trigonometric function).
  *
  *  @return  arbitrary precision numerical cos(x). */
-numeric cos(numeric const & x)
+numeric cos(const numeric & x)
 {
     return ::cos(*x.value);  // -> CLN
 }
@@ -964,7 +964,7 @@ numeric cos(numeric const & x)
 /** Numeric tangent (trigonometric function).
  *
  *  @return  arbitrary precision numerical tan(x). */
-numeric tan(numeric const & x)
+numeric tan(const numeric & x)
 {
     return ::tan(*x.value);  // -> CLN
 }
@@ -972,7 +972,7 @@ numeric tan(numeric const & x)
 /** Numeric inverse sine (trigonometric function).
  *
  *  @return  arbitrary precision numerical asin(x). */
-numeric asin(numeric const & x)
+numeric asin(const numeric & x)
 {
     return ::asin(*x.value);  // -> CLN
 }
@@ -980,7 +980,7 @@ numeric asin(numeric const & x)
 /** Numeric inverse cosine (trigonometric function).
  *
  *  @return  arbitrary precision numerical acos(x). */
-numeric acos(numeric const & x)
+numeric acos(const numeric & x)
 {
     return ::acos(*x.value);  // -> CLN
 }
@@ -990,7 +990,7 @@ numeric acos(numeric const & x)
  *  @param z complex number
  *  @return atan(z)
  *  @exception overflow_error (logarithmic singularity) */
-numeric atan(numeric const & x)
+numeric atan(const numeric & x)
 {
     if (!x.is_real() &&
         x.real().is_zero() &&
@@ -1004,7 +1004,7 @@ numeric atan(numeric const & x)
  *  @param x real number
  *  @param y real number
  *  @return atan(y/x) */
-numeric atan(numeric const & y, numeric const & x)
+numeric atan(const numeric & y, const numeric & x)
 {
     if (x.is_real() && y.is_real())
         return ::atan(realpart(*x.value), realpart(*y.value));  // -> CLN
@@ -1015,7 +1015,7 @@ numeric atan(numeric const & y, numeric const & x)
 /** Numeric hyperbolic sine (trigonometric function).
  *
  *  @return  arbitrary precision numerical sinh(x). */
-numeric sinh(numeric const & x)
+numeric sinh(const numeric & x)
 {
     return ::sinh(*x.value);  // -> CLN
 }
@@ -1023,7 +1023,7 @@ numeric sinh(numeric const & x)
 /** Numeric hyperbolic cosine (trigonometric function).
  *
  *  @return  arbitrary precision numerical cosh(x). */
-numeric cosh(numeric const & x)
+numeric cosh(const numeric & x)
 {
     return ::cosh(*x.value);  // -> CLN
 }
@@ -1031,7 +1031,7 @@ numeric cosh(numeric const & x)
 /** Numeric hyperbolic tangent (trigonometric function).
  *
  *  @return  arbitrary precision numerical tanh(x). */
-numeric tanh(numeric const & x)
+numeric tanh(const numeric & x)
 {
     return ::tanh(*x.value);  // -> CLN
 }
@@ -1039,7 +1039,7 @@ numeric tanh(numeric const & x)
 /** Numeric inverse hyperbolic sine (trigonometric function).
  *
  *  @return  arbitrary precision numerical asinh(x). */
-numeric asinh(numeric const & x)
+numeric asinh(const numeric & x)
 {
     return ::asinh(*x.value);  // -> CLN
 }
@@ -1047,7 +1047,7 @@ numeric asinh(numeric const & x)
 /** Numeric inverse hyperbolic cosine (trigonometric function).
  *
  *  @return  arbitrary precision numerical acosh(x). */
-numeric acosh(numeric const & x)
+numeric acosh(const numeric & x)
 {
     return ::acosh(*x.value);  // -> CLN
 }
@@ -1055,14 +1055,14 @@ numeric acosh(numeric const & x)
 /** Numeric inverse hyperbolic tangent (trigonometric function).
  *
  *  @return  arbitrary precision numerical atanh(x). */
-numeric atanh(numeric const & x)
+numeric atanh(const numeric & x)
 {
     return ::atanh(*x.value);  // -> CLN
 }
 
 /** Numeric evaluation of Riemann's Zeta function.  Currently works only for
  *  integer arguments. */
-numeric zeta(numeric const & x)
+numeric zeta(const numeric & x)
 {
     // A dirty hack to allow for things like zeta(3.0), since CLN currently
     // only knows about integer arguments and zeta(3).evalf() automatically
@@ -1082,7 +1082,7 @@ numeric zeta(numeric const & x)
 
 /** The gamma function.
  *  This is only a stub! */
-numeric gamma(numeric const & x)
+numeric gamma(const numeric & x)
 {
     clog << "gamma(" << x
          << "): Does anybody know good way to calculate this numerically?"
@@ -1092,7 +1092,7 @@ numeric gamma(numeric const & x)
 
 /** The psi function (aka polygamma function).
  *  This is only a stub! */
-numeric psi(numeric const & x)
+numeric psi(const numeric & x)
 {
     clog << "psi(" << x
          << "): Does anybody know good way to calculate this numerically?"
@@ -1102,7 +1102,7 @@ numeric psi(numeric const & x)
 
 /** The psi functions (aka polygamma functions).
  *  This is only a stub! */
-numeric psi(numeric const & n, numeric const & x)
+numeric psi(const numeric & n, const numeric & x)
 {
     clog << "psi(" << n << "," << x
          << "): Does anybody know good way to calculate this numerically?"
@@ -1113,7 +1113,7 @@ numeric psi(numeric const & n, numeric const & x)
 /** Factorial combinatorial function.
  *
  *  @exception range_error (argument must be integer >= 0) */
-numeric factorial(numeric const & nn)
+numeric factorial(const numeric & nn)
 {
     if (!nn.is_nonneg_integer())
         throw (std::range_error("numeric::factorial(): argument must be integer >= 0"));
@@ -1126,7 +1126,7 @@ numeric factorial(numeric const & nn)
  *  @param n  integer argument >= -1
  *  @return n!! == n * (n-2) * (n-4) * ... * ({1|2}) with 0!! == 1 == (-1)!!
  *  @exception range_error (argument must be integer >= -1) */
-numeric doublefactorial(numeric const & nn)
+numeric doublefactorial(const numeric & nn)
 {
     // META-NOTE:  The whole shit here will become obsolete and may be moved
     // out once CLN learns about double factorial, which should be as soon as
@@ -1190,7 +1190,7 @@ numeric doublefactorial(numeric const & nn)
  *  integer n and k and positive n this is the number of ways of choosing k
  *  objects from n distinct objects.  If n is negative, the formula
  *  binomial(n,k) == (-1)^k*binomial(k-n-1,k) is used to compute the result. */
-numeric binomial(numeric const & n, numeric const & k)
+numeric binomial(const numeric & n, const numeric & k)
 {
     if (n.is_integer() && k.is_integer()) {
         if (n.is_nonneg_integer()) {
@@ -1212,7 +1212,7 @@ numeric binomial(numeric const & n, numeric const & k)
  *
  *  @return the nth Bernoulli number (a rational number).
  *  @exception range_error (argument must be integer >= 0) */
-numeric bernoulli(numeric const & nn)
+numeric bernoulli(const numeric & nn)
 {
     if (!nn.is_integer() || nn.is_negative())
         throw (std::range_error("numeric::bernoulli(): argument must be integer >= 0"));
@@ -1250,7 +1250,7 @@ numeric bernoulli(numeric const & nn)
 }
 
 /** Absolute value. */
-numeric abs(numeric const & x)
+numeric abs(const numeric & x)
 {
     return ::abs(*x.value);  // -> CLN
 }
@@ -1262,7 +1262,7 @@ numeric abs(numeric const & x)
  *
  *  @return a mod b in the range [0,abs(b)-1] with sign of b if both are
  *  integer, 0 otherwise. */
-numeric mod(numeric const & a, numeric const & b)
+numeric mod(const numeric & a, const numeric & b)
 {
     if (a.is_integer() && b.is_integer())
         return ::mod(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
@@ -1274,7 +1274,7 @@ numeric mod(numeric const & a, numeric const & b)
  *  Equivalent to Maple's mods.
  *
  *  @return a mod b in the range [-iquo(abs(m)-1,2), iquo(abs(m),2)]. */
-numeric smod(numeric const & a, numeric const & b)
+numeric smod(const numeric & a, const numeric & b)
 {
     //  FIXME: Should this become a member function?
     if (a.is_integer() && b.is_integer()) {
@@ -1290,7 +1290,7 @@ numeric smod(numeric const & a, numeric const & b)
  *  sign of a or is zero.
  *
  *  @return remainder of a/b if both are integer, 0 otherwise. */
-numeric irem(numeric const & a, numeric const & b)
+numeric irem(const numeric & a, const numeric & b)
 {
     if (a.is_integer() && b.is_integer())
         return ::rem(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
@@ -1305,7 +1305,7 @@ numeric irem(numeric const & a, numeric const & b)
  *
  *  @return remainder of a/b and quotient stored in q if both are integer,
  *  0 otherwise. */
-numeric irem(numeric const & a, numeric const & b, numeric & q)
+numeric irem(const numeric & a, const numeric & b, numeric & q)
 {
     if (a.is_integer() && b.is_integer()) {  // -> CLN
         cl_I_div_t rem_quo = truncate2(The(cl_I)(*a.value), The(cl_I)(*b.value));
@@ -1322,7 +1322,7 @@ numeric irem(numeric const & a, numeric const & b, numeric & q)
  *  Equivalent to Maple's iquo as far as sign conventions are concerned.
  *  
  *  @return truncated quotient of a/b if both are integer, 0 otherwise. */
-numeric iquo(numeric const & a, numeric const & b)
+numeric iquo(const numeric & a, const numeric & b)
 {
     if (a.is_integer() && b.is_integer())
         return truncate1(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
@@ -1336,7 +1336,7 @@ numeric iquo(numeric const & a, numeric const & b)
  *
  *  @return truncated quotient of a/b and remainder stored in r if both are
  *  integer, 0 otherwise. */
-numeric iquo(numeric const & a, numeric const & b, numeric & r)
+numeric iquo(const numeric & a, const numeric & b, numeric & r)
 {
     if (a.is_integer() && b.is_integer()) {  // -> CLN
         cl_I_div_t rem_quo = truncate2(The(cl_I)(*a.value), The(cl_I)(*b.value));
@@ -1356,13 +1356,13 @@ numeric iquo(numeric const & a, numeric const & b, numeric & r)
  *  @return square root of z. Branch cut along negative real axis, the negative
  *  real axis itself where imag(z)==0 and real(z)<0 belongs to the upper part
  *  where imag(z)>0. */
-numeric sqrt(numeric const & z)
+numeric sqrt(const numeric & z)
 {
     return ::sqrt(*z.value);  // -> CLN
 }
 
 /** Integer numeric square root. */
-numeric isqrt(numeric const & x)
+numeric isqrt(const numeric & x)
 {
     if (x.is_integer()) {
         cl_I root;
@@ -1376,7 +1376,7 @@ numeric isqrt(numeric const & x)
  *   
  *  @return  The GCD of two numbers if both are integer, a numerical 1
  *  if they are not. */
-numeric gcd(numeric const & a, numeric const & b)
+numeric gcd(const numeric & a, const numeric & b)
 {
     if (a.is_integer() && b.is_integer())
         return ::gcd(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
@@ -1388,7 +1388,7 @@ numeric gcd(numeric const & a, numeric const & b)
  *   
  *  @return  The LCM of two numbers if both are integer, the product of those
  *  two numbers if they are not. */
-numeric lcm(numeric const & a, numeric const & b)
+numeric lcm(const numeric & a, const numeric & b)
 {
     if (a.is_integer() && b.is_integer())
         return ::lcm(The(cl_I)(*a.value), The(cl_I)(*b.value));  // -> CLN
index 27153df..30d4edf 100644 (file)
@@ -64,34 +64,34 @@ private:
 class numeric : public basic
 {
 // friends
-    friend numeric exp(numeric const & x);
-    friend numeric log(numeric const & x);
-    friend numeric sin(numeric const & x);
-    friend numeric cos(numeric const & x);
-    friend numeric tan(numeric const & x);
-    friend numeric asin(numeric const & x);
-    friend numeric acos(numeric const & x);
-    friend numeric atan(numeric const & x);
-    friend numeric atan(numeric const & y, numeric const & x);
-    friend numeric sinh(numeric const & x);
-    friend numeric cosh(numeric const & x);
-    friend numeric tanh(numeric const & x);
-    friend numeric asinh(numeric const & x);
-    friend numeric acosh(numeric const & x);
-    friend numeric atanh(numeric const & x);
-    friend numeric zeta(numeric const & x);
-    friend numeric bernoulli(numeric const & n);
-    friend numeric abs(numeric const & x);
-    friend numeric mod(numeric const & a, numeric const & b);
-    friend numeric smod(numeric const & a, numeric const & b);
-    friend numeric irem(numeric const & a, numeric const & b);
-    friend numeric irem(numeric const & a, numeric const & b, numeric & q);
-    friend numeric iquo(numeric const & a, numeric const & b);
-    friend numeric iquo(numeric const & a, numeric const & b, numeric & r);
-    friend numeric sqrt(numeric const & x);
-    friend numeric isqrt(numeric const & x);
-    friend numeric gcd(numeric const & a, numeric const & b);
-    friend numeric lcm(numeric const & a, numeric const & b);
+    friend numeric exp(const numeric & x);
+    friend numeric log(const numeric & x);
+    friend numeric sin(const numeric & x);
+    friend numeric cos(const numeric & x);
+    friend numeric tan(const numeric & x);
+    friend numeric asin(const numeric & x);
+    friend numeric acos(const numeric & x);
+    friend numeric atan(const numeric & x);
+    friend numeric atan(const numeric & y, const numeric & x);
+    friend numeric sinh(const numeric & x);
+    friend numeric cosh(const numeric & x);
+    friend numeric tanh(const numeric & x);
+    friend numeric asinh(const numeric & x);
+    friend numeric acosh(const numeric & x);
+    friend numeric atanh(const numeric & x);
+    friend numeric zeta(const numeric & x);
+    friend numeric bernoulli(const numeric & n);
+    friend numeric abs(const numeric & x);
+    friend numeric mod(const numeric & a, const numeric & b);
+    friend numeric smod(const numeric & a, const numeric & b);
+    friend numeric irem(const numeric & a, const numeric & b);
+    friend numeric irem(const numeric & a, const numeric & b, numeric & q);
+    friend numeric iquo(const numeric & a, const numeric & b);
+    friend numeric iquo(const numeric & a, const numeric & b, numeric & r);
+    friend numeric sqrt(const numeric & x);
+    friend numeric isqrt(const numeric & x);
+    friend numeric gcd(const numeric & a, const numeric & b);
+    friend numeric lcm(const numeric & a, const numeric & b);
 
 // member functions
 
@@ -100,10 +100,10 @@ class numeric : public basic
 public:
     numeric();
     ~numeric();
-    numeric(numeric const & other);
-    numeric const & operator=(numeric const & other);
+    numeric(const numeric & other);
+    const numeric & operator=(const numeric & other);
 protected:
-    void copy(numeric const & other);
+    void copy(const numeric & other);
     void destroy(bool call_parent);
 
     // other constructors
@@ -129,7 +129,7 @@ public:
     ex diff(symbol const & s) const;
     ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
     numeric integer_content(void) const;
-    ex smod(numeric const &xi) const;
+    ex smod(const numeric &xi) const;
     numeric max_coefficient(void) const;
 protected:
     int compare_same_type(basic const & other) const;
@@ -144,33 +144,26 @@ protected:
 
     // non-virtual functions in this class
 public:
-    numeric add(numeric const & other) const;
-    numeric sub(numeric const & other) const;
-    numeric mul(numeric const & other) const;
-    numeric div(numeric const & other) const;
-    numeric power(numeric const & other) const;
-    numeric const & add_dyn(numeric const & other) const;
-    numeric const & sub_dyn(numeric const & other) const;
-    numeric const & mul_dyn(numeric const & other) const;
-    numeric const & div_dyn(numeric const & other) const;
-    numeric const & power_dyn(numeric const & other) const;
-    numeric const & operator=(int i);
-    numeric const & operator=(unsigned int i);
-    numeric const & operator=(long i);
-    numeric const & operator=(unsigned long i);
-    numeric const & operator=(double d);
-    numeric const & operator=(char const * s);
-    /*
-    numeric add_dyn(numeric const & other) const   { return add(other);   }
-    numeric sub_dyn(numeric const & other) const   { return sub(other);   }
-    numeric mul_dyn(numeric const & other) const   { return mul(other);   }
-    numeric div_dyn(numeric const & other) const   { return div(other);   }
-    numeric power_dyn(numeric const & other) const { return power(other); }
-    */
+    numeric add(const numeric & other) const;
+    numeric sub(const numeric & other) const;
+    numeric mul(const numeric & other) const;
+    numeric div(const numeric & other) const;
+    numeric power(const numeric & other) const;
+    const numeric & add_dyn(const numeric & other) const;
+    const numeric & sub_dyn(const numeric & other) const;
+    const numeric & mul_dyn(const numeric & other) const;
+    const numeric & div_dyn(const numeric & other) const;
+    const numeric & power_dyn(const numeric & other) const;
+    const numeric & operator=(int i);
+    const numeric & operator=(unsigned int i);
+    const numeric & operator=(long i);
+    const numeric & operator=(unsigned long i);
+    const numeric & operator=(double d);
+    const numeric & operator=(char const * s);
     numeric inverse(void) const;
     int csgn(void) const;
-    int compare(numeric const & other) const;
-    bool is_equal(numeric const & other) const;
+    int compare(const numeric & other) const;
+    bool is_equal(const numeric & other) const;
     bool is_zero(void) const;
     bool is_positive(void) const;
     bool is_negative(void) const;
@@ -184,12 +177,12 @@ public:
     bool is_real(void) const;
     bool is_cinteger(void) const;
     bool is_crational(void) const;
-    bool operator==(numeric const & other) const;
-    bool operator!=(numeric const & other) const;
-    bool operator<(numeric const & other) const;
-    bool operator<=(numeric const & other) const;
-    bool operator>(numeric const & other) const;
-    bool operator>=(numeric const & other) const;
+    bool operator==(const numeric & other) const;
+    bool operator!=(const numeric & other) const;
+    bool operator<(const numeric & other) const;
+    bool operator<=(const numeric & other) const;
+    bool operator>(const numeric & other) const;
+    bool operator>=(const numeric & other) const;
     int to_int(void) const;
     double to_double(void) const;
     numeric real(void) const;
@@ -217,42 +210,42 @@ extern _numeric_digits Digits;
 
 // global functions
 
-numeric exp(numeric const & x);
-numeric log(numeric const & x);
-numeric sin(numeric const & x);
-numeric cos(numeric const & x);
-numeric tan(numeric const & x);
-numeric asin(numeric const & x);
-numeric acos(numeric const & x);
-numeric atan(numeric const & x);
-numeric atan(numeric const & y, numeric const & x);
-numeric sinh(numeric const & x);
-numeric cosh(numeric const & x);
-numeric tanh(numeric const & x);
-numeric asinh(numeric const & x);
-numeric acosh(numeric const & x);
-numeric atanh(numeric const & x);
-numeric zeta(numeric const & x);
-numeric gamma(numeric const & x);
-numeric psi(numeric const & x);
-numeric psi(numeric const & n, numeric const & x);
-numeric factorial(numeric const & n);
-numeric doublefactorial(numeric const & n);
-numeric binomial(numeric const & n, numeric const & k);
-numeric bernoulli(numeric const & n);
-
-numeric abs(numeric const & x);
-numeric mod(numeric const & a, numeric const & b);
-numeric smod(numeric const & a, numeric const & b);
-numeric irem(numeric const & a, numeric const & b);
-numeric irem(numeric const & a, numeric const & b, numeric & q);
-numeric iquo(numeric const & a, numeric const & b);
-numeric iquo(numeric const & a, numeric const & b, numeric & r);
-numeric sqrt(numeric const & x);
-numeric isqrt(numeric const & x);
-
-numeric gcd(numeric const & a, numeric const & b);
-numeric lcm(numeric const & a, numeric const & b);
+numeric exp(const numeric & x);
+numeric log(const numeric & x);
+numeric sin(const numeric & x);
+numeric cos(const numeric & x);
+numeric tan(const numeric & x);
+numeric asin(const numeric & x);
+numeric acos(const numeric & x);
+numeric atan(const numeric & x);
+numeric atan(const numeric & y, const numeric & x);
+numeric sinh(const numeric & x);
+numeric cosh(const numeric & x);
+numeric tanh(const numeric & x);
+numeric asinh(const numeric & x);
+numeric acosh(const numeric & x);
+numeric atanh(const numeric & x);
+numeric zeta(const numeric & x);
+numeric gamma(const numeric & x);
+numeric psi(const numeric & x);
+numeric psi(const numeric & n, const numeric & x);
+numeric factorial(const numeric & n);
+numeric doublefactorial(const numeric & n);
+numeric binomial(const numeric & n, const numeric & k);
+numeric bernoulli(const numeric & n);
+
+numeric abs(const numeric & x);
+numeric mod(const numeric & a, const numeric & b);
+numeric smod(const numeric & a, const numeric & b);
+numeric irem(const numeric & a, const numeric & b);
+numeric irem(const numeric & a, const numeric & b, numeric & q);
+numeric iquo(const numeric & a, const numeric & b);
+numeric iquo(const numeric & a, const numeric & b, numeric & r);
+numeric sqrt(const numeric & x);
+numeric isqrt(const numeric & x);
+
+numeric gcd(const numeric & a, const numeric & b);
+numeric lcm(const numeric & a, const numeric & b);
 
 /** Exception thrown by numeric members to signal failure */
 struct numeric_fail
@@ -262,58 +255,61 @@ struct numeric_fail
 };
 
 // wrapper functions around member functions
-inline numeric inverse(numeric const & x)
+inline numeric pow(const numeric & x, const numeric & y)
+{ return x.power(y); }
+
+inline numeric inverse(const numeric & x)
 { return x.inverse(); }
 
-inline bool csgn(numeric const & x)
+inline bool csgn(const numeric & x)
 { return x.csgn(); }
 
-inline bool is_zero(numeric const & x)
+inline bool is_zero(const numeric & x)
 { return x.is_zero(); }
 
-inline bool is_positive(numeric const & x)
+inline bool is_positive(const numeric & x)
 { return x.is_positive(); }
 
-inline bool is_integer(numeric const & x)
+inline bool is_integer(const numeric & x)
 { return x.is_integer(); }
 
-inline bool is_pos_integer(numeric const & x)
+inline bool is_pos_integer(const numeric & x)
 { return x.is_pos_integer(); }
 
-inline bool is_nonneg_integer(numeric const & x)
+inline bool is_nonneg_integer(const numeric & x)
 { return x.is_nonneg_integer(); }
 
-inline bool is_even(numeric const & x)
+inline bool is_even(const numeric & x)
 { return x.is_even(); }
 
-inline bool is_odd(numeric const & x)
+inline bool is_odd(const numeric & x)
 { return x.is_odd(); }
 
-inline bool is_prime(numeric const & x)
+inline bool is_prime(const numeric & x)
 { return x.is_prime(); }
 
-inline bool is_rational(numeric const & x)
+inline bool is_rational(const numeric & x)
 { return x.is_rational(); }
 
-inline bool is_real(numeric const & x)
+inline bool is_real(const numeric & x)
 { return x.is_real(); }
 
-inline bool is_cinteger(numeric const & x)
+inline bool is_cinteger(const numeric & x)
 { return x.is_cinteger(); }
 
-inline bool is_crational(numeric const & x)
+inline bool is_crational(const numeric & x)
 { return x.is_crational(); }
 
-inline numeric real(numeric const & x)
+inline numeric real(const numeric & x)
 { return x.real(); }
 
-inline numeric imag(numeric const & x)
+inline numeric imag(const numeric & x)
 { return x.imag(); }
 
-inline numeric numer(numeric const & x)
+inline numeric numer(const numeric & x)
 { return x.numer(); }
 
-inline numeric denom(numeric const & x)
+inline numeric denom(const numeric & x)
 { return x.denom(); }
 
 // numeric evaluation functions for class constant objects:
index a60da5d..f226ff1 100644 (file)
@@ -41,10 +41,10 @@ ex operator/(ex const & lh, ex const & rh);
 ex operator%(ex const & lh, ex const & rh); // non-commutative multiplication
 
 // binary arithmetic operators numeric with numeric
-numeric operator+(numeric const & lh, numeric const & rh);
-numeric operator-(numeric const & lh, numeric const & rh);
-numeric operator*(numeric const & lh, numeric const & rh);
-numeric operator/(numeric const & lh, numeric const & rh);
+numeric operator+(const numeric & lh, const numeric & rh);
+numeric operator-(const numeric & lh, const numeric & rh);
+numeric operator*(const numeric & lh, const numeric & rh);
+numeric operator/(const numeric & lh, const numeric & rh);
 
 // binary arithmetic assignment operators with ex
 ex const & operator+=(ex & lh, ex const & rh);
@@ -54,17 +54,17 @@ ex const & operator/=(ex & lh, ex const & rh);
 ex const & operator%=(ex & lh, ex const & rh); // non-commutative multiplication
 
 // binary arithmetic assignment operators with numeric
-numeric const & operator+=(numeric & lh, numeric const & rh);
-numeric const & operator-=(numeric & lh, numeric const & rh);
-numeric const & operator*=(numeric & lh, numeric const & rh);
-numeric const & operator/=(numeric & lh, numeric const & rh);
+const numeric & operator+=(numeric & lh, const numeric & rh);
+const numeric & operator-=(numeric & lh, const numeric & rh);
+const numeric & operator*=(numeric & lh, const numeric & rh);
+const numeric & operator/=(numeric & lh, const numeric & rh);
 
 // unary operators
 ex operator+(ex const & lh);
 ex operator-(ex const & lh);
 
-numeric operator+(numeric const & lh);
-numeric operator-(numeric const & lh);
+numeric operator+(const numeric & lh);
+numeric operator-(const numeric & lh);
 numeric& operator++(numeric & rh);
 numeric& operator--(numeric & rh);
 numeric operator++(numeric & lh, int);
index 51e5bfc..dfd1a47 100644 (file)
@@ -120,11 +120,15 @@ basic * power::duplicate() const
 void power::print(ostream & os, unsigned upper_precedence) const
 {
     debugmsg("power print",LOGLEVEL_PRINT);
-    if (precedence<=upper_precedence) os << "(";
-    basis.print(os,precedence);
-    os << "^";
-    exponent.print(os,precedence);
-    if (precedence<=upper_precedence) os << ")";
+    if (exponent.is_equal(_ex1_2())) {
+        os << "sqrt(" << basis << ")";
+    } else {
+        if (precedence<=upper_precedence) os << "(";
+        basis.print(os,precedence);
+        os << "^";
+        exponent.print(os,precedence);
+        if (precedence<=upper_precedence) os << ")";
+    }
 }
 
 void power::printraw(ostream & os) const
index e113c04..1a52936 100644 (file)
@@ -53,6 +53,118 @@ int compare_pointers(void const * a, void const * b)
 // `construct on first use' chest of numbers
 //////////
 
+// numeric -60
+numeric const & _num_60(void)
+{
+    const static ex e = ex(numeric(-60));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_60(void)
+{
+    static ex * e = new ex(_num_60());
+    return *e;
+}
+
+// numeric -120
+numeric const & _num_120(void)
+{
+    const static ex e = ex(numeric(-120));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_120(void)
+{
+    static ex * e = new ex(_num_120());
+    return *e;
+}
+
+// numeric -30
+numeric const & _num_30(void)
+{
+    const static ex e = ex(numeric(-30));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_30(void)
+{
+    static ex * e = new ex(_num_30());
+    return *e;
+}
+
+// numeric -25
+numeric const & _num_25(void)
+{
+    const static ex e = ex(numeric(-25));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_25(void)
+{
+    static ex * e = new ex(_num_25());
+    return *e;
+}
+
+// numeric -24
+numeric const & _num_24(void)
+{
+    const static ex e = ex(numeric(-24));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_24(void)
+{
+    static ex * e = new ex(_num_24());
+    return *e;
+}
+
+// numeric -20
+numeric const & _num_20(void)
+{
+    const static ex e = ex(numeric(-20));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_20(void)
+{
+    static ex * e = new ex(_num_20());
+    return *e;
+}
+
+// numeric -18
+numeric const & _num_18(void)
+{
+    const static ex e = ex(numeric(-18));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_18(void)
+{
+    static ex * e = new ex(_num_18());
+    return *e;
+}
+
+// numeric -15
+numeric const & _num_15(void)
+{
+    const static ex e = ex(numeric(-15));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_15(void)
+{
+    static ex * e = new ex(_num_15());
+    return *e;
+}
+
 // numeric -12
 numeric const & _num_12(void)
 {
@@ -249,6 +361,20 @@ ex const & _ex_1_3(void)
     return *e;
 }    
 
+// numeric -1/4
+numeric const & _num_1_4(void)
+{
+    const static ex e = ex(numeric(-1,4));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex_1_4(void)
+{
+    static ex * e = new ex(_num_1_4());
+    return *e;
+}    
+
 // numeric  0
 numeric const & _num0(void)
 {
@@ -263,6 +389,20 @@ ex const & _ex0(void)
     return *e;
 }
 
+// numeric  1/4
+numeric const & _num1_4(void)
+{
+    const static ex e = ex(numeric(1,4));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex1_4(void)
+{
+    static ex * e = new ex(_num1_4());
+    return *e;
+}    
+
 // numeric  1/3
 numeric const & _num1_3(void)
 {
@@ -459,6 +599,118 @@ ex const & _ex12(void)
     return *e;
 }
 
+// numeric  15
+numeric const & _num15(void)
+{
+    const static ex e = ex(numeric(15));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex15(void)
+{
+    static ex * e = new ex(_num15());
+    return *e;
+}
+
+// numeric  18
+numeric const & _num18(void)
+{
+    const static ex e = ex(numeric(18));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex18(void)
+{
+    static ex * e = new ex(_num18());
+    return *e;
+}
+
+// numeric  20
+numeric const & _num20(void)
+{
+    const static ex e = ex(numeric(20));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex20(void)
+{
+    static ex * e = new ex(_num20());
+    return *e;
+}
+
+// numeric  24
+numeric const & _num24(void)
+{
+    const static ex e = ex(numeric(24));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex24(void)
+{
+    static ex * e = new ex(_num24());
+    return *e;
+}
+
+// numeric  25
+numeric const & _num25(void)
+{
+    const static ex e = ex(numeric(25));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex25(void)
+{
+    static ex * e = new ex(_num25());
+    return *e;
+}
+
+// numeric  30
+numeric const & _num30(void)
+{
+    const static ex e = ex(numeric(30));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex30(void)
+{
+    static ex * e = new ex(_num30());
+    return *e;
+}
+
+// numeric  60
+numeric const & _num60(void)
+{
+    const static ex e = ex(numeric(60));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex60(void)
+{
+    static ex * e = new ex(_num60());
+    return *e;
+}
+
+// numeric  120
+numeric const & _num120(void)
+{
+    const static ex e = ex(numeric(120));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+ex const & _ex120(void)
+{
+    static ex * e = new ex(_num120());
+    return *e;
+}
+
 // comment skeleton for header files
 
 
index 0f6082a..133e9f6 100644 (file)
@@ -146,65 +146,101 @@ OutputIterator mymerge3(InputIterator1 first1, InputIterator1 last1,
 
 class numeric;
 class ex;
-                       
+
+numeric const & _num_120(void);   // -120
+ex const & _ex_120(void);
+numeric const & _num_60(void);    // -60
+ex const & _ex_60(void);
+numeric const & _num_30(void);    // -30
+ex const & _ex_30(void);
+numeric const & _num_25(void);    // -25
+ex const & _ex_25(void);
+numeric const & _num_24(void);    // -24
+ex const & _ex_24(void);
+numeric const & _num_20(void);    // -20
+ex const & _ex_20(void);
+numeric const & _num_18(void);    // -18
+ex const & _ex_18(void);
+numeric const & _num_15(void);    // -15
+ex const & _ex_15(void);
 numeric const & _num_12(void);    // -12
+ex const & _ex_12(void);
 numeric const & _num_11(void);    // -11
+ex const & _ex_11(void);
 numeric const & _num_10(void);    // -10
+ex const & _ex_10(void);
 numeric const & _num_9(void);     // -9
+ex const & _ex_9(void);
 numeric const & _num_8(void);     // -8
+ex const & _ex_8(void);
 numeric const & _num_7(void);     // -7
+ex const & _ex_7(void);
 numeric const & _num_6(void);     // -6
+ex const & _ex_6(void);
 numeric const & _num_5(void);     // -5
+ex const & _ex_5(void);
 numeric const & _num_4(void);     // -4
+ex const & _ex_4(void);
 numeric const & _num_3(void);     // -3
+ex const & _ex_3(void);
 numeric const & _num_2(void);     // -2
+ex const & _ex_2(void);
 numeric const & _num_1(void);     // -1
+ex const & _ex_1(void);
 numeric const & _num_1_2(void);   // -1/2
+ex const & _ex_1_2(void);
 numeric const & _num_1_3(void);   // -1/3
+ex const & _ex_1_3(void);
+numeric const & _num_1_4(void);   // -1/4
+ex const & _ex_1_4(void);
 numeric const & _num0(void);      //  0
+ex const & _ex0(void);
+numeric const & _num1_4(void);    //  1/4
+ex const & _ex1_4(void);
 numeric const & _num1_3(void);    //  1/3
+ex const & _ex1_3(void);
 numeric const & _num1_2(void);    //  1/2
+ex const & _ex1_2(void);
 numeric const & _num1(void);      //  1
+ex const & _ex1(void);
 numeric const & _num2(void);      //  2
+ex const & _ex2(void);
 numeric const & _num3(void);      //  3
+ex const & _ex3(void);
 numeric const & _num4(void);      //  4
+ex const & _ex4(void);
 numeric const & _num5(void);      //  5
+ex const & _ex5(void);
 numeric const & _num6(void);      //  6
+ex const & _ex6(void);
 numeric const & _num7(void);      //  7
+ex const & _ex7(void);
 numeric const & _num8(void);      //  8
+ex const & _ex8(void);
 numeric const & _num9(void);      //  9
+ex const & _ex9(void);
 numeric const & _num10(void);     //  10
+ex const & _ex10(void);
 numeric const & _num11(void);     //  11
+ex const & _ex11(void);
 numeric const & _num12(void);     //  12
-ex const & _ex_12(void);          // -12
-ex const & _ex_11(void);          // -11
-ex const & _ex_10(void);          // -10
-ex const & _ex_9(void);           // -9
-ex const & _ex_8(void);           // -8
-ex const & _ex_7(void);           // -7
-ex const & _ex_6(void);           // -6
-ex const & _ex_5(void);           // -5
-ex const & _ex_4(void);           // -4
-ex const & _ex_3(void);           // -3
-ex const & _ex_2(void);           // -2
-ex const & _ex_1(void);           // -1
-ex const & _ex_1_2(void);         // -1/2
-ex const & _ex_1_3(void);         // -1/3
-ex const & _ex0(void);            //  0
-ex const & _ex1_3(void);          //  1/3
-ex const & _ex1_2(void);          //  1/2
-ex const & _ex1(void);            //  1
-ex const & _ex2(void);            //  2
-ex const & _ex3(void);            //  3
-ex const & _ex4(void);            //  4
-ex const & _ex5(void);            //  5
-ex const & _ex6(void);            //  6
-ex const & _ex7(void);            //  7
-ex const & _ex8(void);            //  8
-ex const & _ex9(void);            //  9
-ex const & _ex10(void);           //  10
-ex const & _ex11(void);           //  11
-ex const & _ex12(void);           //  12
+ex const & _ex12(void);
+numeric const & _num15(void);     //  15
+ex const & _ex15(void);
+numeric const & _num18(void);     //  18
+ex const & _ex18(void);
+numeric const & _num20(void);     //  20
+ex const & _ex20(void);
+numeric const & _num24(void);     //  24
+ex const & _ex24(void);
+numeric const & _num25(void);     //  25
+ex const & _ex25(void);
+numeric const & _num30(void);     //  30
+ex const & _ex30(void);
+numeric const & _num60(void);     //  60
+ex const & _ex60(void);
+numeric const & _num120(void);    //  120
+ex const & _ex120(void);
 
 #ifndef NO_GINAC_NAMESPACE
 } // namespace GiNaC