From: Richard Kreckel <Richard.Kreckel@uni-mainz.de>
Date: Sat, 25 Dec 1999 21:09:06 +0000 (+0000)
Subject: - speedup by declaring x_pt and so on const
X-Git-Tag: release_0-5-0~78
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=a6bb52b00bf185271774e7d56215923700a3ec40;hp=68acc9698d0f6cdff0c2c6d1ea3b622916009f71

- speedup by declaring x_pt and so on const
- added some code for series expansion of beta function (untested!)
---

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index aba77986..f69d4510 100644
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -1433,10 +1433,14 @@ When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
 @cindex Laurent expansion
 
 Expressions know how to expand themselves as a Taylor series or (more
-generally) a Laurent series.  Similar to most conventional Computer
-Algebra Systems, no distinction is made between those two.  There is a
-class of its own for storing such series as well as a class for storing
-the order of the series.  A sample program could read:
+generally) a Laurent series.  As in most conventional Computer Algebra
+Systems, no distinction is made between those two.  There is a class of
+its own for storing such series as well as a class for storing the order
+of the series.  As a consequence, if you want to work with series,
+i.e. multiply two series, you need to call the method @code{ex::series}
+again to convert it to a series object with the usual structure
+(expansion plus order term).  A sample application from special
+relativity could read:
 
 @example
 #include <ginac/ginac.h>
@@ -1444,25 +1448,28 @@ using namespace GiNaC;
 
 int main()
 @{
-    symbol x("x");
-    numeric point(0);
-    ex MyExpr1 = sin(x);
-    ex MyExpr2 = 1/(x - pow(x, 2) - pow(x, 3));
-    ex MyTailor, MySeries;
+    symbol v("v"), c("c");
+    
+    ex gamma = 1/sqrt(1 - pow(v/c,2));
+    ex mass_nonrel = gamma.series(v, 0, 10);
+    
+    cout << "the relativistic mass increase with v is " << endl
+         << mass_nonrel << endl;
+    
+    cout << "the inverse square of this series is " << endl
+         << pow(mass_nonrel,-2).series(v, 0, 10) << endl;
     
-    MyTailor = MyExpr1.series(x, point, 5);
-    cout << MyExpr1 << " == " << MyTailor
-         << " for small " << x << endl;
-    MySeries = MyExpr2.series(x, point, 7);
-    cout << MyExpr2 << " == " << MySeries
-         << " for small " << x << endl;
     // ...
 @}
 @end example
 
+Only calling the series method makes the last output simplify to
+@math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
+series raised to the power @math{-2}.
+
 @cindex M@'echain's formula
-As an instructive application, let us calculate the numerical value of
-Archimedes' constant
+As another instructive application, let us calculate the numerical 
+value of Archimedes' constant
 @tex
 $\pi$
 @end tex
diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index bcf2539f..8ed38515 100644
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -101,7 +101,7 @@ static ex gamma_diff(const ex & x, unsigned diff_param)
     return psi(x)*gamma(x);
 }
 
-static ex gamma_series(const ex & x, const symbol & s, const ex & point, int order)
+static ex gamma_series(const ex & x, const symbol & s, const ex & pt, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to psi function
@@ -110,17 +110,17 @@ static ex gamma_series(const ex & x, const symbol & s, const ex & point, int ord
     //   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);
-    if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive))
+    //   series(gamma(x+m+1)/(x*(x+1)*...*(x+m)),x,-m,order+1);
+    const ex x_pt = x.subs(s==pt);
+    if (!x_pt.info(info_flags::integer) || x_pt.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:
-    numeric m = -ex_to_numeric(xpoint);
+    numeric m = -ex_to_numeric(x_pt);
     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);
+    return (ser_numer/ser_denom).series(s, pt, order+1);
 }
 
 REGISTER_FUNCTION(gamma, gamma_eval, gamma_evalf, gamma_diff, gamma_series);
@@ -143,11 +143,11 @@ static ex beta_evalf(const ex & x, const ex & 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
         // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y)
+        numeric nx(ex_to_numeric(x));
+        numeric ny(ex_to_numeric(y));
         if (nx.is_real() && nx.is_integer() &&
             ny.is_real() && ny.is_integer()) {
             if (nx.is_negative()) {
@@ -190,23 +190,36 @@ static ex beta_diff(const ex & x, const ex & y, unsigned diff_param)
     return retval;
 }
 
-static ex beta_series(const ex & x, const ex & y, const symbol & s, const ex & point, int order)
+static ex beta_series(const ex & x, const ex & y, const symbol & s, const ex & pt, 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)))
+    // Taylor series where there is no pole of one of the gamma functions
+    // falls back to beta function evaluation.  Otherwise, fall back to
+    // gamma series directly.
+    // FIXME: this could need some testing, maybe it's wrong in some cases?
+    const ex x_pt = x.subs(s==pt);
+    const ex y_pt = y.subs(s==pt);
+    ex x_ser, y_ser, xy_ser;
+    if ((!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) &&
+        (!y_pt.info(info_flags::integer) || y_pt.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(): Mama, please code me!"));
+    // trap the case where x is on a pole directly:
+    if (x.info(info_flags::integer) && !x.info(info_flags::positive))
+        x_ser = gamma(x+s).series(s,pt,order);
+    else
+        x_ser = gamma(x).series(s,pt,order);
+    // trap the case where y is on a pole directly:
+    if (y.info(info_flags::integer) && !y.info(info_flags::positive))
+        y_ser = gamma(y+s).series(s,pt,order);
+    else
+        y_ser = gamma(y).series(s,pt,order);
+    // trap the case where y is on a pole directly:
+    if ((x+y).info(info_flags::integer) && !(x+y).info(info_flags::positive))
+        xy_ser = gamma(y+x+s).series(s,pt,order);
+    else
+        xy_ser = gamma(y+x).series(s,pt,order);
+    // compose the result:
+    return (x_ser*y_ser/xy_ser).series(s,pt,order);
 }
 
 REGISTER_FUNCTION(beta, beta_eval, beta_evalf, beta_diff, beta_series);
@@ -277,7 +290,7 @@ static ex psi1_diff(const ex & x, unsigned diff_param)
     return psi(_ex1(), x);
 }
 
-static ex psi1_series(const ex & x, const symbol & s, const ex & point, int order)
+static ex psi1_series(const ex & x, const symbol & s, const ex & pt, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
@@ -287,15 +300,15 @@ static ex psi1_series(const ex & x, const symbol & s, const ex & point, int orde
     // from which follows
     //   series(psi(x),x,-m,order) ==
     //   series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x,-m,order);
-    ex xpoint = x.subs(s==point);
-    if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive))
+    const ex x_pt = x.subs(s==pt);
+    if (!x_pt.info(info_flags::integer) || x_pt.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:
-    numeric m = -ex_to_numeric(xpoint);
+    numeric m = -ex_to_numeric(x_pt);
     ex recur;
     for (numeric p; p<=m; ++p)
         recur += power(x+p,_ex_1());
-    return (psi(x+m+_ex1())-recur).series(s, point, order);
+    return (psi(x+m+_ex1())-recur).series(s, pt, order);
 }
 
 const unsigned function_index_psi1 = function::register_new("psi", psi1_eval, psi1_evalf, psi1_diff, psi1_series);
@@ -391,7 +404,7 @@ static ex psi2_diff(const ex & n, const ex & x, unsigned diff_param)
     return psi(n+_ex1(), x);
 }
 
-static ex psi2_series(const ex & n, const ex & x, const symbol & s, const ex & point, int order)
+static ex psi2_series(const ex & n, const ex & x, const symbol & s, const ex & pt, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
@@ -402,16 +415,16 @@ static ex psi2_series(const ex & n, const ex & x, const symbol & s, const ex & p
     //   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);
-    ex xpoint = x.subs(s==point);
-    if (!xpoint.info(info_flags::integer) || xpoint.info(info_flags::positive))
+    const ex x_pt = x.subs(s==pt);
+    if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive))
         throw do_taylor();  // caught by function::series()
     // if we got here we have to care for a pole of order n+1 at -m:
-    numeric m = -ex_to_numeric(xpoint);
+    numeric m = -ex_to_numeric(x_pt);
     ex recur;
     for (numeric p; p<=m; ++p)
         recur += power(x+p,-n+_ex_1());
     recur *= factorial(n)*power(_ex_1(),n);
-    return (psi(n, x+m+_ex1())-recur).series(s, point, order);
+    return (psi(n, x+m+_ex1())-recur).series(s, pt, order);
 }
 
 const unsigned function_index_psi2 = function::register_new("psi", psi2_eval, psi2_evalf, psi2_diff, psi2_series);
diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index 9937e305..eae92385 100644
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -41,7 +41,7 @@ namespace GiNaC {
 // exponential function
 //////////
 
-static ex exp_evalf(ex const & x)
+static ex exp_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -50,7 +50,7 @@ static ex exp_evalf(ex const & x)
     return exp(ex_to_numeric(x)); // -> numeric exp(numeric)
 }
 
-static ex exp_eval(ex const & x)
+static ex exp_eval(const ex & x)
 {
     // exp(0) -> 1
     if (x.is_zero()) {
@@ -80,7 +80,7 @@ static ex exp_eval(ex const & x)
     return exp(x).hold();
 }
 
-static ex exp_diff(ex const & x, unsigned diff_param)
+static ex exp_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
 
@@ -94,7 +94,7 @@ REGISTER_FUNCTION(exp, exp_eval, exp_evalf, exp_diff, NULL);
 // natural logarithm
 //////////
 
-static ex log_evalf(ex const & x)
+static ex log_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -103,7 +103,7 @@ static ex log_evalf(ex const & x)
     return log(ex_to_numeric(x)); // -> numeric log(numeric)
 }
 
-static ex log_eval(ex const & x)
+static ex log_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         if (x.is_equal(_ex1()))  // log(1) -> 0
@@ -120,17 +120,20 @@ static ex log_eval(ex const & x)
         if (!x.info(info_flags::crational))
             return log_evalf(x);
     }
-    // log(exp(t)) -> t (for real-valued t):
+    // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
     if (is_ex_the_function(x, exp)) {
         ex t = x.op(0);
-        if (t.info(info_flags::real))
-            return t;
+        if (t.info(info_flags::numeric)) {
+            numeric nt = ex_to_numeric(t);
+            if (nt.is_real())
+                return t;
+        }
     }
     
     return log(x).hold();
 }
 
-static ex log_diff(ex const & x, unsigned diff_param)
+static ex log_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
 
@@ -144,7 +147,7 @@ REGISTER_FUNCTION(log, log_eval, log_evalf, log_diff, NULL);
 // sine (trigonometric function)
 //////////
 
-static ex sin_evalf(ex const & x)
+static ex sin_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -153,7 +156,7 @@ static ex sin_evalf(ex const & x)
     return sin(ex_to_numeric(x)); // -> numeric sin(numeric)
 }
 
-static ex sin_eval(ex const & x)
+static ex sin_eval(const ex & x)
 {
     // sin(n/d*Pi) -> { all known non-nested radicals }
     ex SixtyExOverPi = _ex60()*x/Pi;
@@ -209,7 +212,7 @@ static ex sin_eval(ex const & x)
     return sin(x).hold();
 }
 
-static ex sin_diff(ex const & x, unsigned diff_param)
+static ex sin_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -223,7 +226,7 @@ REGISTER_FUNCTION(sin, sin_eval, sin_evalf, sin_diff, NULL);
 // cosine (trigonometric function)
 //////////
 
-static ex cos_evalf(ex const & x)
+static ex cos_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -232,7 +235,7 @@ static ex cos_evalf(ex const & x)
     return cos(ex_to_numeric(x)); // -> numeric cos(numeric)
 }
 
-static ex cos_eval(ex const & x)
+static ex cos_eval(const ex & x)
 {
     // cos(n/d*Pi) -> { all known non-nested radicals }
     ex SixtyExOverPi = _ex60()*x/Pi;
@@ -288,7 +291,7 @@ static ex cos_eval(ex const & x)
     return cos(x).hold();
 }
 
-static ex cos_diff(ex const & x, unsigned diff_param)
+static ex cos_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
 
@@ -302,7 +305,7 @@ REGISTER_FUNCTION(cos, cos_eval, cos_evalf, cos_diff, NULL);
 // tangent (trigonometric function)
 //////////
 
-static ex tan_evalf(ex const & x)
+static ex tan_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -311,7 +314,7 @@ static ex tan_evalf(ex const & x)
     return tan(ex_to_numeric(x));
 }
 
-static ex tan_eval(ex const & x)
+static ex tan_eval(const ex & x)
 {
     // tan(n/d*Pi) -> { all known non-nested radicals }
     ex SixtyExOverPi = _ex60()*x/Pi;
@@ -364,7 +367,7 @@ static ex tan_eval(ex const & x)
     return tan(x).hold();
 }
 
-static ex tan_diff(ex const & x, unsigned diff_param)
+static ex tan_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -372,16 +375,16 @@ static ex tan_diff(ex const & x, unsigned diff_param)
     return (_ex1()+power(tan(x),_ex2()));
 }
 
-static ex tan_series(ex const & x, symbol const & s, ex const & point, int order)
+static ex tan_series(const ex & x, const symbol & s, const ex & pt, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to tan_diff.
     // On a pole simply expand sin(x)/cos(x).
-    ex xpoint = x.subs(s==point);
-    if (!(2*xpoint/Pi).info(info_flags::odd))
+    const ex x_pt = x.subs(s==pt);
+    if (!(2*x_pt/Pi).info(info_flags::odd))
         throw do_taylor();  // caught by function::series()
     // if we got here we have to care for a simple pole
-    return (sin(x)/cos(x)).series(s, point, order+2);
+    return (sin(x)/cos(x)).series(s, pt, order+2);
 }
 
 REGISTER_FUNCTION(tan, tan_eval, tan_evalf, tan_diff, tan_series);
@@ -390,7 +393,7 @@ REGISTER_FUNCTION(tan, tan_eval, tan_evalf, tan_diff, tan_series);
 // inverse sine (arc sine)
 //////////
 
-static ex asin_evalf(ex const & x)
+static ex asin_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -399,7 +402,7 @@ static ex asin_evalf(ex const & x)
     return asin(ex_to_numeric(x)); // -> numeric asin(numeric)
 }
 
-static ex asin_eval(ex const & x)
+static ex asin_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         // asin(0) -> 0
@@ -425,7 +428,7 @@ static ex asin_eval(ex const & x)
     return asin(x).hold();
 }
 
-static ex asin_diff(ex const & x, unsigned diff_param)
+static ex asin_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -439,7 +442,7 @@ REGISTER_FUNCTION(asin, asin_eval, asin_evalf, asin_diff, NULL);
 // inverse cosine (arc cosine)
 //////////
 
-static ex acos_evalf(ex const & x)
+static ex acos_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -448,7 +451,7 @@ static ex acos_evalf(ex const & x)
     return acos(ex_to_numeric(x)); // -> numeric acos(numeric)
 }
 
-static ex acos_eval(ex const & x)
+static ex acos_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         // acos(1) -> 0
@@ -474,7 +477,7 @@ static ex acos_eval(ex const & x)
     return acos(x).hold();
 }
 
-static ex acos_diff(ex const & x, unsigned diff_param)
+static ex acos_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -488,7 +491,7 @@ REGISTER_FUNCTION(acos, acos_eval, acos_evalf, acos_diff, NULL);
 // inverse tangent (arc tangent)
 //////////
 
-static ex atan_evalf(ex const & x)
+static ex atan_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -497,7 +500,7 @@ static ex atan_evalf(ex const & x)
     return atan(ex_to_numeric(x)); // -> numeric atan(numeric)
 }
 
-static ex atan_eval(ex const & x)
+static ex atan_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         // atan(0) -> 0
@@ -511,7 +514,7 @@ static ex atan_eval(ex const & x)
     return atan(x).hold();
 }    
 
-static ex atan_diff(ex const & x, unsigned diff_param)
+static ex atan_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
 
@@ -525,7 +528,7 @@ REGISTER_FUNCTION(atan, atan_eval, atan_evalf, atan_diff, NULL);
 // inverse tangent (atan2(y,x))
 //////////
 
-static ex atan2_evalf(ex const & y, ex const & x)
+static ex atan2_evalf(const ex & y, const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(y,numeric)
@@ -535,7 +538,7 @@ static ex atan2_evalf(ex const & y, ex const & x)
     return atan(ex_to_numeric(y),ex_to_numeric(x)); // -> numeric atan(numeric)
 }
 
-static ex atan2_eval(ex const & y, ex const & x)
+static ex atan2_eval(const ex & y, const ex & x)
 {
     if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
         x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
@@ -545,7 +548,7 @@ static ex atan2_eval(ex const & y, ex const & x)
     return atan2(y,x).hold();
 }    
 
-static ex atan2_diff(ex const & y, ex const & x, unsigned diff_param)
+static ex atan2_diff(const ex & y, const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param<2);
     
@@ -563,7 +566,7 @@ REGISTER_FUNCTION(atan2, atan2_eval, atan2_evalf, atan2_diff, NULL);
 // hyperbolic sine (trigonometric function)
 //////////
 
-static ex sinh_evalf(ex const & x)
+static ex sinh_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -572,7 +575,7 @@ static ex sinh_evalf(ex const & x)
     return sinh(ex_to_numeric(x)); // -> numeric sinh(numeric)
 }
 
-static ex sinh_eval(ex const & x)
+static ex sinh_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         if (x.is_zero())  // sinh(0) -> 0
@@ -601,7 +604,7 @@ static ex sinh_eval(ex const & x)
     return sinh(x).hold();
 }
 
-static ex sinh_diff(ex const & x, unsigned diff_param)
+static ex sinh_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -615,7 +618,7 @@ REGISTER_FUNCTION(sinh, sinh_eval, sinh_evalf, sinh_diff, NULL);
 // hyperbolic cosine (trigonometric function)
 //////////
 
-static ex cosh_evalf(ex const & x)
+static ex cosh_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -624,7 +627,7 @@ static ex cosh_evalf(ex const & x)
     return cosh(ex_to_numeric(x)); // -> numeric cosh(numeric)
 }
 
-static ex cosh_eval(ex const & x)
+static ex cosh_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         if (x.is_zero())  // cosh(0) -> 1
@@ -653,7 +656,7 @@ static ex cosh_eval(ex const & x)
     return cosh(x).hold();
 }
 
-static ex cosh_diff(ex const & x, unsigned diff_param)
+static ex cosh_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -667,7 +670,7 @@ REGISTER_FUNCTION(cosh, cosh_eval, cosh_evalf, cosh_diff, NULL);
 // hyperbolic tangent (trigonometric function)
 //////////
 
-static ex tanh_evalf(ex const & x)
+static ex tanh_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -676,7 +679,7 @@ static ex tanh_evalf(ex const & x)
     return tanh(ex_to_numeric(x)); // -> numeric tanh(numeric)
 }
 
-static ex tanh_eval(ex const & x)
+static ex tanh_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         if (x.is_zero())  // tanh(0) -> 0
@@ -705,7 +708,7 @@ static ex tanh_eval(ex const & x)
     return tanh(x).hold();
 }
 
-static ex tanh_diff(ex const & x, unsigned diff_param)
+static ex tanh_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -713,16 +716,16 @@ static ex tanh_diff(ex const & x, unsigned diff_param)
     return _ex1()-power(tanh(x),_ex2());
 }
 
-static ex tanh_series(ex const & x, symbol const & s, ex const & point, int order)
+static ex tanh_series(const ex & x, const symbol & s, const ex & pt, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to tanh_diff.
     // On a pole simply expand sinh(x)/cosh(x).
-    ex xpoint = x.subs(s==point);
-    if (!(2*I*xpoint/Pi).info(info_flags::odd))
+    const ex x_pt = x.subs(s==pt);
+    if (!(2*I*x_pt/Pi).info(info_flags::odd))
         throw do_taylor();  // caught by function::series()
     // if we got here we have to care for a simple pole
-    return (sinh(x)/cosh(x)).series(s, point, order+2);
+    return (sinh(x)/cosh(x)).series(s, pt, order+2);
 }
 
 REGISTER_FUNCTION(tanh, tanh_eval, tanh_evalf, tanh_diff, tanh_series);
@@ -731,7 +734,7 @@ REGISTER_FUNCTION(tanh, tanh_eval, tanh_evalf, tanh_diff, tanh_series);
 // inverse hyperbolic sine (trigonometric function)
 //////////
 
-static ex asinh_evalf(ex const & x)
+static ex asinh_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -740,7 +743,7 @@ static ex asinh_evalf(ex const & x)
     return asinh(ex_to_numeric(x)); // -> numeric asinh(numeric)
 }
 
-static ex asinh_eval(ex const & x)
+static ex asinh_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         // asinh(0) -> 0
@@ -754,7 +757,7 @@ static ex asinh_eval(ex const & x)
     return asinh(x).hold();
 }
 
-static ex asinh_diff(ex const & x, unsigned diff_param)
+static ex asinh_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -768,7 +771,7 @@ REGISTER_FUNCTION(asinh, asinh_eval, asinh_evalf, asinh_diff, NULL);
 // inverse hyperbolic cosine (trigonometric function)
 //////////
 
-static ex acosh_evalf(ex const & x)
+static ex acosh_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -777,7 +780,7 @@ static ex acosh_evalf(ex const & x)
     return acosh(ex_to_numeric(x)); // -> numeric acosh(numeric)
 }
 
-static ex acosh_eval(ex const & x)
+static ex acosh_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         // acosh(0) -> Pi*I/2
@@ -797,7 +800,7 @@ static ex acosh_eval(ex const & x)
     return acosh(x).hold();
 }
 
-static ex acosh_diff(ex const & x, unsigned diff_param)
+static ex acosh_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -811,7 +814,7 @@ REGISTER_FUNCTION(acosh, acosh_eval, acosh_evalf, acosh_diff, NULL);
 // inverse hyperbolic tangent (trigonometric function)
 //////////
 
-static ex atanh_evalf(ex const & x)
+static ex atanh_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
        TYPECHECK(x,numeric)
@@ -820,7 +823,7 @@ static ex atanh_evalf(ex const & x)
     return atanh(ex_to_numeric(x)); // -> numeric atanh(numeric)
 }
 
-static ex atanh_eval(ex const & x)
+static ex atanh_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         // atanh(0) -> 0
@@ -837,7 +840,7 @@ static ex atanh_eval(ex const & x)
     return atanh(x).hold();
 }
 
-static ex atanh_diff(ex const & x, unsigned diff_param)
+static ex atanh_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
diff --git a/ginac/inifcns_zeta.cpp b/ginac/inifcns_zeta.cpp
index dfd4638e..a8b23f15 100644
--- a/ginac/inifcns_zeta.cpp
+++ b/ginac/inifcns_zeta.cpp
@@ -39,7 +39,7 @@ namespace GiNaC {
 // Riemann's Zeta-function
 //////////
 
-static ex zeta1_evalf(ex const & x)
+static ex zeta1_evalf(const ex & x)
 {
     BEGIN_TYPECHECK
         TYPECHECK(x,numeric)
@@ -48,7 +48,7 @@ static ex zeta1_evalf(ex const & x)
     return zeta(ex_to_numeric(x));
 }
 
-static ex zeta1_eval(ex const & x)
+static ex zeta1_eval(const ex & x)
 {
     if (x.info(info_flags::numeric)) {
         numeric y = ex_to_numeric(x);
@@ -74,7 +74,7 @@ static ex zeta1_eval(ex const & x)
     return zeta(x).hold();
 }
 
-static ex zeta1_diff(ex const & x, unsigned diff_param)
+static ex zeta1_diff(const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param==0);
     
@@ -87,7 +87,7 @@ const unsigned function_index_zeta1 = function::register_new("zeta", zeta1_eval,
 // Derivatives of Riemann's Zeta-function  zeta(0,x)==zeta(x)
 //////////
 
-static ex zeta2_eval(ex const & n, ex const & x)
+static ex zeta2_eval(const ex & n, const ex & x)
 {
     if (n.info(info_flags::numeric)) {
         // zeta(0,x) -> zeta(x)
@@ -98,7 +98,7 @@ static ex zeta2_eval(ex const & n, ex const & x)
     return zeta(n, x).hold();
 }
 
-static ex zeta2_diff(ex const & n, ex const & x, unsigned diff_param)
+static ex zeta2_diff(const ex & n, const ex & x, unsigned diff_param)
 {
     GINAC_ASSERT(diff_param<2);
     
diff --git a/ginac/series.cpp b/ginac/series.cpp
index 4189b098..de849df1 100644
--- a/ginac/series.cpp
+++ b/ginac/series.cpp
@@ -264,7 +264,7 @@ ex basic::series(symbol const & s, ex const & point, int order) const
             // Series terminates
             return series::series(s, point, seq);
         }
-        coeff = power(fac, -1) * deriv.subs(s == point);
+        coeff = fac.inverse() * deriv.subs(s == point);
         if (!coeff.is_zero())
             seq.push_back(expair(coeff, numeric(n)));
     }