- Some corrections in series expansion of several functions.

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 96f4c9a596869ab47886ed8eac5ab61dc693fd17..c03f4a15a6e29b39d4dea7d5269711fc5ce5f3a4 100644 (file)
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -65,6 +65,51 @@ static ex abs_eval(const ex & x)
 REGISTER_FUNCTION(abs, eval_func(abs_eval).
                        evalf_func(abs_evalf));
 
+
+//////////
+// Complex sign
+//////////
+
+static ex csgn_evalf(const ex & x)
+{
+    BEGIN_TYPECHECK
+        TYPECHECK(x,numeric)
+    END_TYPECHECK(csgn(x))
+    
+    return csgn(ex_to_numeric(x));
+}
+
+static ex csgn_eval(const ex & x)
+{
+    if (is_ex_exactly_of_type(x, numeric))
+        return csgn(ex_to_numeric(x));
+    
+    if (is_ex_exactly_of_type(x, mul)) {
+        numeric oc = ex_to_numeric(x.op(x.nops()-1));
+        if (oc.is_real()) {
+            if (oc > 0)
+                // csgn(42*x) -> csgn(x)
+                return csgn(x/oc).hold();
+            else
+                // csgn(-42*x) -> -csgn(x)
+                return -csgn(x/oc).hold();
+        }
+        if (oc.real().is_zero()) {
+            if (oc.imag() > 0)
+                // csgn(42*I*x) -> csgn(I*x)
+                return csgn(I*x/oc).hold();
+            else
+                // csgn(-42*I*x) -> -csgn(I*x)
+                return -csgn(I*x/oc).hold();
+        }
+    }
+    
+    return csgn(x).hold();
+}
+
+REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
+                        evalf_func(csgn_evalf));
+
 //////////
 // dilogarithm
 //////////

diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 75599ea9bfe5e759dedff6ec631bb819b87bfad6..9781d935c1a2e777a2dfbf9326fba1da8eae8d22 100644 (file)
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -32,6 +32,9 @@ namespace GiNaC {
 
 /** Absolute value. */
 DECLARE_FUNCTION_1P(abs)
+    
+/** Complex sign. */
+DECLARE_FUNCTION_1P(csgn)
 
 /** Sine. */
 DECLARE_FUNCTION_1P(sin)

diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp
index b28c246253258e0980a1fbbdc7ad50eefdc769ef..85275f62aabec258d301d794bcc764a36f92ce21 100644 (file)
--- a/ginac/inifcns_gamma.cpp
+++ b/ginac/inifcns_gamma.cpp
@@ -84,11 +84,33 @@ static ex lgamma_deriv(const ex & x, unsigned deriv_param)
     return psi(x);
 }
 
-// need to implement lgamma_series.
+
+static ex lgamma_series(const ex & x, const relational & rel, int order)
+{
+    // method:
+    // Taylor series where there is no pole falls back to psi function
+    // evaluation.
+    // On a pole at -m use the recurrence relation
+    //   lgamma(x) == lgamma(x+1)-log(x)
+    // from which follows
+    //   series(lgamma(x),x==-m,order) ==
+    //   series(lgamma(x+m+1)-log(x)...-log(x+m)),x==-m,order+1);
+    const ex x_pt = x.subs(rel);
+    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(x_pt);
+    ex ser_sub = _ex0();
+    for (numeric p; p<=m; ++p)
+        ser_sub += log(x+p);
+    return (lgamma(x+m+_ex1())-ser_sub).series(rel, order);
+}
+
 
 REGISTER_FUNCTION(lgamma, eval_func(lgamma_eval).
                           evalf_func(lgamma_evalf).
-                          derivative_func(lgamma_deriv));
+                          derivative_func(lgamma_deriv).
+                          series_func(lgamma_series));
 
 
 //////////
@@ -156,7 +178,7 @@ static ex tgamma_deriv(const ex & x, unsigned deriv_param)
 }
 
 
-static ex tgamma_series(const ex & x, const relational & r, int order)
+static ex tgamma_series(const ex & x, const relational & rel, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to psi function
@@ -164,9 +186,9 @@ static ex tgamma_series(const ex & x, const relational & r, int order)
     // On a pole at -m use the recurrence relation
     //   tgamma(x) == tgamma(x+1) / x
     // from which follows
-    //   series(tgamma(x),x,-m,order) ==
-    //   series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x,-m,order+1);
-    const ex x_pt = x.subs(r);
+    //   series(tgamma(x),x==-m,order) ==
+    //   series(tgamma(x+m+1)/(x*(x+1)*...*(x+m)),x==-m,order+1);
+    const ex x_pt = x.subs(rel);
     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:
@@ -174,7 +196,7 @@ static ex tgamma_series(const ex & x, const relational & r, int order)
     ex ser_denom = _ex1();
     for (numeric p; p<=m; ++p)
         ser_denom *= x+p;
-    return (tgamma(x+m+_ex1())/ser_denom).series(r, order+1);
+    return (tgamma(x+m+_ex1())/ser_denom).series(rel, order+1);
 }
 
 
@@ -251,38 +273,37 @@ static ex beta_deriv(const ex & x, const ex & y, unsigned deriv_param)
 }
 
 
-static ex beta_series(const ex & x, const ex & y, const relational & r, int order)
+static ex beta_series(const ex & x, const ex & y, const relational & rel, int order)
 {
     // method:
     // Taylor series where there is no pole of one of the tgamma functions
     // falls back to beta function evaluation.  Otherwise, fall back to
     // tgamma series directly.
-    // FIXME: this could need some testing, maybe it's wrong in some cases?
-    const ex x_pt = x.subs(r);
-    const ex y_pt = y.subs(r);
-    GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
-    const symbol *s = static_cast<symbol *>(r.lhs().bp);
+    const ex x_pt = x.subs(rel);
+    const ex y_pt = y.subs(rel);
+    GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
+    const symbol *s = static_cast<symbol *>(rel.lhs().bp);
     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()
     // trap the case where x is on a pole directly:
     if (x.info(info_flags::integer) && !x.info(info_flags::positive))
-        x_ser = tgamma(x+*s).series(r,order);
+        x_ser = tgamma(x+*s).series(rel,order);
     else
-        x_ser = tgamma(x).series(r,order);
+        x_ser = tgamma(x).series(rel,order);
     // trap the case where y is on a pole directly:
     if (y.info(info_flags::integer) && !y.info(info_flags::positive))
-        y_ser = tgamma(y+*s).series(r,order);
+        y_ser = tgamma(y+*s).series(rel,order);
     else
-        y_ser = tgamma(y).series(r,order);
+        y_ser = tgamma(y).series(rel,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 = tgamma(y+x+*s).series(r,order);
+        xy_ser = tgamma(y+x+*s).series(rel,order);
     else
-        xy_ser = tgamma(y+x).series(r,order);
+        xy_ser = tgamma(y+x).series(rel,order);
     // compose the result:
-    return (x_ser*y_ser/xy_ser).series(r,order);
+    return (x_ser*y_ser/xy_ser).series(rel,order);
 }
 
 
@@ -358,7 +379,7 @@ static ex psi1_deriv(const ex & x, unsigned deriv_param)
     return psi(_ex1(), x);
 }
 
-static ex psi1_series(const ex & x, const relational & r, int order)
+static ex psi1_series(const ex & x, const relational & rel, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
@@ -366,9 +387,9 @@ static ex psi1_series(const ex & x, const relational & r, int order)
     // On a pole at -m use the recurrence relation
     //   psi(x) == psi(x+1) - 1/z
     // 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);
-    const ex x_pt = x.subs(r);
+    //   series(psi(x),x==-m,order) ==
+    //   series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x==-m,order);
+    const ex x_pt = x.subs(rel);
     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:
@@ -376,7 +397,7 @@ static ex psi1_series(const ex & x, const relational & r, int order)
     ex recur;
     for (numeric p; p<=m; ++p)
         recur += power(x+p,_ex_1());
-    return (psi(x+m+_ex1())-recur).series(r, order);
+    return (psi(x+m+_ex1())-recur).series(rel, order);
 }
 
 const unsigned function_index_psi1 =
@@ -478,7 +499,7 @@ static ex psi2_deriv(const ex & n, const ex & x, unsigned deriv_param)
     return psi(n+_ex1(), x);
 }
 
-static ex psi2_series(const ex & n, const ex & x, const relational & r, int order)
+static ex psi2_series(const ex & n, const ex & x, const relational & rel, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to polygamma function
@@ -486,10 +507,10 @@ static ex psi2_series(const ex & n, const ex & x, const relational & r, int orde
     // On a pole at -m use the recurrence relation
     //   psi(n,x) == psi(n,x+1) - (-)^n * n! / x^(n+1)
     // from which follows
-    //   series(psi(x),x,-m,order) == 
+    //   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);
-    const ex x_pt = x.subs(r);
+    //                                      ... + (x+m)^(-n-1))),x==-m,order);
+    const ex x_pt = x.subs(rel);
     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:
@@ -498,7 +519,7 @@ static ex psi2_series(const ex & n, const ex & x, const relational & r, int orde
     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(r, order);
+    return (psi(n, x+m+_ex1())-recur).series(rel, order);
 }
 
 const unsigned function_index_psi2 =

diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp
index 12fabc0e98499ea5e4c11ff196f753b6d32b3872..8705901464af978c08a8042a5e754ff5903ff921 100644 (file)
--- a/ginac/inifcns_trans.cpp
+++ b/ginac/inifcns_trans.cpp
@@ -144,6 +144,29 @@ static ex log_deriv(const ex & x, unsigned deriv_param)
     return power(x, _ex_1());
 }
 
+/*static ex log_series(const ex &x, const relational &rel, int order)
+{
+    const ex x_pt = x.subs(rel);
+    if (!x_pt.info(info_flags::negative) && !x_pt.is_zero())
+        throw do_taylor();  // caught by function::series()
+    // now we either have to care for the branch cut or the branch point:
+    if (x_pt.is_zero()) {  // branch point: return a plain log(x).
+        epvector seq;
+        seq.push_back(expair(log(x), _ex0()));
+        return pseries(rel, seq);
+    } // the branch cut:
+    const ex point = rel.rhs();
+    const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+    const symbol foo;
+    // compute the formal series:
+    ex replx = series(log(x),*s==foo,order).subs(foo==point);
+    epvector seq;
+    // FIXME: this is probably off by 2 or so:
+    seq.push_back(expair(-I*csgn(x*I)*Pi,_ex0()));
+    seq.push_back(expair(Order(_ex1()),order));
+    return series(replx + pseries(rel, seq),rel,order);
+}*/
+
 static ex log_series(const ex &x, const relational &r, int order)
 {
        if (x.subs(r).is_zero()) {
@@ -395,16 +418,16 @@ static ex tan_deriv(const ex & x, unsigned deriv_param)
     return (_ex1()+power(tan(x),_ex2()));
 }
 
-static ex tan_series(const ex &x, const relational &r, int order)
+static ex tan_series(const ex &x, const relational &rel, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to tan_deriv.
     // On a pole simply expand sin(x)/cos(x).
-    const ex x_pt = x.subs(r);
+    const ex x_pt = x.subs(rel);
     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(r, order+2);
+    return (sin(x)/cos(x)).series(rel, order+2);
 }
 
 REGISTER_FUNCTION(tan, eval_func(tan_eval).
@@ -752,16 +775,16 @@ static ex tanh_deriv(const ex & x, unsigned deriv_param)
     return _ex1()-power(tanh(x),_ex2());
 }
 
-static ex tanh_series(const ex &x, const relational &r, int order)
+static ex tanh_series(const ex &x, const relational &rel, int order)
 {
     // method:
     // Taylor series where there is no pole falls back to tanh_deriv.
     // On a pole simply expand sinh(x)/cosh(x).
-    const ex x_pt = x.subs(r);
+    const ex x_pt = x.subs(rel);
     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(r, order+2);
+    return (sinh(x)/cosh(x)).series(rel, order+2);
 }
 
 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).