]> www.ginac.de Git - ginac.git/blobdiff - ginac/inifcns_trans.cpp
- ginac-config.1: removed (duplicate).
[ginac.git] / ginac / inifcns_trans.cpp
index c7dbfbe6b0508240d0666b1770c18f782568ab0c..a626ce1c29332370ad13e471554e3cf4f5260ba0 100644 (file)
@@ -34,9 +34,7 @@
 #include "pseries.h"
 #include "utils.h"
 
-#ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
-#endif // ndef NO_NAMESPACE_GINAC
 
 //////////
 // exponential function
@@ -44,11 +42,10 @@ namespace GiNaC {
 
 static ex exp_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-               TYPECHECK(x,numeric)
-       END_TYPECHECK(exp(x))
+       if (is_exactly_a<numeric>(x))
+               return exp(ex_to<numeric>(x));
        
-       return exp(ex_to_numeric(x)); // -> numeric exp(numeric)
+       return exp(x).hold();
 }
 
 static ex exp_eval(const ex & x)
@@ -58,9 +55,9 @@ static ex exp_eval(const ex & x)
                return _ex1();
        }
        // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
-       ex TwoExOverPiI=(_ex2()*x)/(Pi*I);
+       const ex TwoExOverPiI=(_ex2()*x)/(Pi*I);
        if (TwoExOverPiI.info(info_flags::integer)) {
-               numeric z=mod(ex_to_numeric(TwoExOverPiI),_num4());
+               numeric z=mod(ex_to<numeric>(TwoExOverPiI),_num4());
                if (z.is_equal(_num0()))
                        return _ex1();
                if (z.is_equal(_num1()))
@@ -76,7 +73,7 @@ static ex exp_eval(const ex & x)
        
        // exp(float)
        if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
-               return exp_evalf(x);
+               return exp(ex_to<numeric>(x));
        
        return exp(x).hold();
 }
@@ -91,7 +88,8 @@ static ex exp_deriv(const ex & x, unsigned deriv_param)
 
 REGISTER_FUNCTION(exp, eval_func(exp_eval).
                        evalf_func(exp_evalf).
-                       derivative_func(exp_deriv));
+                       derivative_func(exp_deriv).
+                       latex_name("\\exp"));
 
 //////////
 // natural logarithm
@@ -99,17 +97,16 @@ REGISTER_FUNCTION(exp, eval_func(exp_eval).
 
 static ex log_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-               TYPECHECK(x,numeric)
-       END_TYPECHECK(log(x))
+       if (is_exactly_a<numeric>(x))
+               return log(ex_to<numeric>(x));
        
-       return log(ex_to_numeric(x)); // -> numeric log(numeric)
+       return log(x).hold();
 }
 
 static ex log_eval(const ex & x)
 {
        if (x.info(info_flags::numeric)) {
-               if (x.is_equal(_ex0()))  // log(0) -> infinity
+               if (x.is_zero())         // log(0) -> infinity
                        throw(pole_error("log_eval(): log(0)",0));
                if (x.info(info_flags::real) && x.info(info_flags::negative))
                        return (log(-x)+I*Pi);
@@ -121,13 +118,13 @@ static ex log_eval(const ex & x)
                        return (Pi*I*_num_1_2());
                // log(float)
                if (!x.info(info_flags::crational))
-                       return log_evalf(x);
+                       return log(ex_to<numeric>(x));
        }
        // 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::numeric)) {
-                       numeric nt = ex_to_numeric(t);
+                       numeric nt = ex_to<numeric>(t);
                        if (nt.is_real())
                                return t;
                }
@@ -158,7 +155,7 @@ static ex log_series(const ex &arg,
        } catch (pole_error) {
                must_expand_arg = true;
        }
-       // or we are at the branch cut anyways
+       // or we are at the branch point anyways
        if (arg_pt.is_zero())
                must_expand_arg = true;
        
@@ -171,28 +168,36 @@ static ex log_series(const ex &arg,
                // Return a plain n*log(x) for the x^n part and series expand the
                // other part.  Add them together and reexpand again in order to have
                // one unnested pseries object.  All this also works for negative n.
-               const pseries argser = ex_to_pseries(arg.series(rel, order, options));
+               const pseries argser = ex_to<pseries>(arg.series(rel, order, options));
                const symbol *s = static_cast<symbol *>(rel.lhs().bp);
                const ex point = rel.rhs();
                const int n = argser.ldegree(*s);
                epvector seq;
-               seq.push_back(expair(n*log(*s-point), _ex0()));
+               // construct what we carelessly called the n*log(x) term above
+               ex coeff = argser.coeff(*s, n);
+               // expand the log, but only if coeff is real and > 0, since otherwise
+               // it would make the branch cut run into the wrong direction
+               if (coeff.info(info_flags::positive))
+                       seq.push_back(expair(n*log(*s-point)+log(coeff), _ex0()));
+               else
+                       seq.push_back(expair(log(coeff*pow(*s-point, n)), _ex0()));
                if (!argser.is_terminating() || argser.nops()!=1) {
                        // in this case n more terms are needed
-                       ex newarg = ex_to_pseries(arg.series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
-                       return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, options)));
+                       // (sadly, to generate them, we have to start from the beginning)
+                       ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
+                       return pseries(rel, seq).add_series(ex_to<pseries>(log(newarg).series(rel, order, options)));
                } else  // it was a monomial
                        return pseries(rel, seq);
        }
        if (!(options & series_options::suppress_branchcut) &&
-                arg_pt.info(info_flags::negative)) {
+            arg_pt.info(info_flags::negative)) {
                // method:
                // This is the branch cut: assemble the primitive series manually and
                // then add the corresponding complex step function.
                const symbol *s = static_cast<symbol *>(rel.lhs().bp);
                const ex point = rel.rhs();
                const symbol foo;
-               ex replarg = series(log(arg), *s==foo, order, false).subs(foo==point);
+               const ex replarg = series(log(arg), *s==foo, order).subs(foo==point);
                epvector seq;
                seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
                seq.push_back(expair(Order(_ex1()), order));
@@ -204,7 +209,8 @@ static ex log_series(const ex &arg,
 REGISTER_FUNCTION(log, eval_func(log_eval).
                        evalf_func(log_evalf).
                        derivative_func(log_deriv).
-                       series_func(log_series));
+                       series_func(log_series).
+                       latex_name("\\ln"));
 
 //////////
 // sine (trigonometric function)
@@ -212,20 +218,19 @@ REGISTER_FUNCTION(log, eval_func(log_eval).
 
 static ex sin_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(sin(x))
+       if (is_exactly_a<numeric>(x))
+               return sin(ex_to<numeric>(x));
        
-       return sin(ex_to_numeric(x)); // -> numeric sin(numeric)
+       return sin(x).hold();
 }
 
 static ex sin_eval(const ex & x)
 {
        // sin(n/d*Pi) -> { all known non-nested radicals }
-       ex SixtyExOverPi = _ex60()*x/Pi;
+       const ex SixtyExOverPi = _ex60()*x/Pi;
        ex sign = _ex1();
        if (SixtyExOverPi.info(info_flags::integer)) {
-               numeric z = mod(ex_to_numeric(SixtyExOverPi),_num120());
+               numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
                if (z>=_num60()) {
                        // wrap to interval [0, Pi)
                        z -= _num60();
@@ -255,7 +260,7 @@ static ex sin_eval(const ex & x)
                        return sign*_ex1();
        }
        
-       if (is_ex_exactly_of_type(x, function)) {
+       if (is_exactly_a<function>(x)) {
                ex t = x.op(0);
                // sin(asin(x)) -> x
                if (is_ex_the_function(x, asin))
@@ -270,7 +275,7 @@ static ex sin_eval(const ex & x)
        
        // sin(float) -> float
        if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
-               return sin_evalf(x);
+               return sin(ex_to<numeric>(x));
        
        return sin(x).hold();
 }
@@ -285,7 +290,8 @@ static ex sin_deriv(const ex & x, unsigned deriv_param)
 
 REGISTER_FUNCTION(sin, eval_func(sin_eval).
                        evalf_func(sin_evalf).
-                       derivative_func(sin_deriv));
+                       derivative_func(sin_deriv).
+                       latex_name("\\sin"));
 
 //////////
 // cosine (trigonometric function)
@@ -293,20 +299,19 @@ REGISTER_FUNCTION(sin, eval_func(sin_eval).
 
 static ex cos_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-               TYPECHECK(x,numeric)
-       END_TYPECHECK(cos(x))
+       if (is_exactly_a<numeric>(x))
+               return cos(ex_to<numeric>(x));
        
-       return cos(ex_to_numeric(x)); // -> numeric cos(numeric)
+       return cos(x).hold();
 }
 
 static ex cos_eval(const ex & x)
 {
        // cos(n/d*Pi) -> { all known non-nested radicals }
-       ex SixtyExOverPi = _ex60()*x/Pi;
+       const ex SixtyExOverPi = _ex60()*x/Pi;
        ex sign = _ex1();
        if (SixtyExOverPi.info(info_flags::integer)) {
-               numeric z = mod(ex_to_numeric(SixtyExOverPi),_num120());
+               numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
                if (z>=_num60()) {
                        // wrap to interval [0, Pi)
                        z = _num120()-z;
@@ -336,7 +341,7 @@ static ex cos_eval(const ex & x)
                        return sign*_ex0();
        }
        
-       if (is_ex_exactly_of_type(x, function)) {
+       if (is_exactly_a<function>(x)) {
                ex t = x.op(0);
                // cos(acos(x)) -> x
                if (is_ex_the_function(x, acos))
@@ -351,7 +356,7 @@ static ex cos_eval(const ex & x)
        
        // cos(float) -> float
        if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
-               return cos_evalf(x);
+               return cos(ex_to<numeric>(x));
        
        return cos(x).hold();
 }
@@ -366,7 +371,8 @@ static ex cos_deriv(const ex & x, unsigned deriv_param)
 
 REGISTER_FUNCTION(cos, eval_func(cos_eval).
                        evalf_func(cos_evalf).
-                       derivative_func(cos_deriv));
+                       derivative_func(cos_deriv).
+                       latex_name("\\cos"));
 
 //////////
 // tangent (trigonometric function)
@@ -374,20 +380,19 @@ REGISTER_FUNCTION(cos, eval_func(cos_eval).
 
 static ex tan_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(tan(x)) // -> numeric tan(numeric)
+       if (is_exactly_a<numeric>(x))
+               return tan(ex_to<numeric>(x));
        
-       return tan(ex_to_numeric(x));
+       return tan(x).hold();
 }
 
 static ex tan_eval(const ex & x)
 {
        // tan(n/d*Pi) -> { all known non-nested radicals }
-       ex SixtyExOverPi = _ex60()*x/Pi;
+       const ex SixtyExOverPi = _ex60()*x/Pi;
        ex sign = _ex1();
        if (SixtyExOverPi.info(info_flags::integer)) {
-               numeric z = mod(ex_to_numeric(SixtyExOverPi),_num60());
+               numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num60());
                if (z>=_num60()) {
                        // wrap to interval [0, Pi)
                        z -= _num60();
@@ -413,7 +418,7 @@ static ex tan_eval(const ex & x)
                        throw (pole_error("tan_eval(): simple pole",1));
        }
        
-       if (is_ex_exactly_of_type(x, function)) {
+       if (is_exactly_a<function>(x)) {
                ex t = x.op(0);
                // tan(atan(x)) -> x
                if (is_ex_the_function(x, atan))
@@ -428,7 +433,7 @@ static ex tan_eval(const ex & x)
        
        // tan(float) -> float
        if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
-               return tan_evalf(x);
+               return tan(ex_to<numeric>(x));
        }
        
        return tan(x).hold();
@@ -461,7 +466,8 @@ static ex tan_series(const ex &x,
 REGISTER_FUNCTION(tan, eval_func(tan_eval).
                        evalf_func(tan_evalf).
                        derivative_func(tan_deriv).
-                       series_func(tan_series));
+                       series_func(tan_series).
+                       latex_name("\\tan"));
 
 //////////
 // inverse sine (arc sine)
@@ -469,11 +475,10 @@ REGISTER_FUNCTION(tan, eval_func(tan_eval).
 
 static ex asin_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(asin(x))
+       if (is_exactly_a<numeric>(x))
+               return asin(ex_to<numeric>(x));
        
-       return asin(ex_to_numeric(x)); // -> numeric asin(numeric)
+       return asin(x).hold();
 }
 
 static ex asin_eval(const ex & x)
@@ -496,7 +501,7 @@ static ex asin_eval(const ex & x)
                        return _num_1_2()*Pi;
                // asin(float) -> float
                if (!x.info(info_flags::crational))
-                       return asin_evalf(x);
+                       return asin(ex_to<numeric>(x));
        }
        
        return asin(x).hold();
@@ -512,7 +517,8 @@ static ex asin_deriv(const ex & x, unsigned deriv_param)
 
 REGISTER_FUNCTION(asin, eval_func(asin_eval).
                         evalf_func(asin_evalf).
-                        derivative_func(asin_deriv));
+                        derivative_func(asin_deriv).
+                        latex_name("\\arcsin"));
 
 //////////
 // inverse cosine (arc cosine)
@@ -520,11 +526,10 @@ REGISTER_FUNCTION(asin, eval_func(asin_eval).
 
 static ex acos_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(acos(x))
+       if (is_exactly_a<numeric>(x))
+               return acos(ex_to<numeric>(x));
        
-       return acos(ex_to_numeric(x)); // -> numeric acos(numeric)
+       return acos(x).hold();
 }
 
 static ex acos_eval(const ex & x)
@@ -547,7 +552,7 @@ static ex acos_eval(const ex & x)
                        return Pi;
                // acos(float) -> float
                if (!x.info(info_flags::crational))
-                       return acos_evalf(x);
+                       return acos(ex_to<numeric>(x));
        }
        
        return acos(x).hold();
@@ -563,7 +568,8 @@ static ex acos_deriv(const ex & x, unsigned deriv_param)
 
 REGISTER_FUNCTION(acos, eval_func(acos_eval).
                         evalf_func(acos_evalf).
-                        derivative_func(acos_deriv));
+                        derivative_func(acos_deriv).
+                        latex_name("\\arccos"));
 
 //////////
 // inverse tangent (arc tangent)
@@ -571,18 +577,17 @@ REGISTER_FUNCTION(acos, eval_func(acos_eval).
 
 static ex atan_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-               TYPECHECK(x,numeric)
-       END_TYPECHECK(atan(x))
+       if (is_exactly_a<numeric>(x))
+               return atan(ex_to<numeric>(x));
        
-       return atan(ex_to_numeric(x)); // -> numeric atan(numeric)
+       return atan(x).hold();
 }
 
 static ex atan_eval(const ex & x)
 {
        if (x.info(info_flags::numeric)) {
                // atan(0) -> 0
-               if (x.is_equal(_ex0()))
+               if (x.is_zero())
                        return _ex0();
                // atan(1) -> Pi/4
                if (x.is_equal(_ex1()))
@@ -594,7 +599,7 @@ static ex atan_eval(const ex & x)
                        throw (pole_error("atan_eval(): logarithmic pole",0));
                // atan(float) -> float
                if (!x.info(info_flags::crational))
-                       return atan_evalf(x);
+                       return atan(ex_to<numeric>(x));
        }
        
        return atan(x).hold();
@@ -608,7 +613,7 @@ static ex atan_deriv(const ex & x, unsigned deriv_param)
        return power(_ex1()+power(x,_ex2()), _ex_1());
 }
 
-static ex atan_series(const ex &x,
+static ex atan_series(const ex &arg,
                       const relational &rel,
                       int order,
                       unsigned options)
@@ -620,35 +625,53 @@ static ex atan_series(const ex &x,
        // one running from -I down the imaginary axis.  The points I and -I are
        // poles.
        // On the branch cuts and the poles series expand
-       //     log((1+I*x)/(1-I*x))/(2*I)
+       //     (log(1+I*x)-log(1-I*x))/(2*I)
        // instead.
-       // (The constant term on the cut itself could be made simpler.)
-       const ex x_pt = x.subs(rel);
-       if (!(I*x_pt).info(info_flags::real))
+       const ex arg_pt = arg.subs(rel);
+       if (!(I*arg_pt).info(info_flags::real))
                throw do_taylor();     // Re(x) != 0
-       if ((I*x_pt).info(info_flags::real) && abs(I*x_pt)<_ex1())
+       if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1())
                throw do_taylor();     // Re(x) == 0, but abs(x)<1
-       // if we got here we have to care for cuts and poles
-       return (log((1+I*x)/(1-I*x))/(2*I)).series(rel, order, options);
+       // care for the poles, using the defining formula for atan()...
+       if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
+               return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
+       if (!(options & series_options::suppress_branchcut)) {
+               // method:
+               // This is the branch cut: assemble the primitive series manually and
+               // then add the corresponding complex step function.
+               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               const ex point = rel.rhs();
+               const symbol foo;
+               const ex replarg = series(atan(arg), *s==foo, order).subs(foo==point);
+               ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2();
+               if ((I*arg_pt)<_ex0())
+                       Order0correction += log((I*arg_pt+_ex_1())/(I*arg_pt+_ex1()))*I*_ex_1_2();
+               else
+                       Order0correction += log((I*arg_pt+_ex1())/(I*arg_pt+_ex_1()))*I*_ex1_2();
+               epvector seq;
+               seq.push_back(expair(Order0correction, _ex0()));
+               seq.push_back(expair(Order(_ex1()), order));
+               return series(replarg - pseries(rel, seq), rel, order);
+       }
+       throw do_taylor();
 }
 
 REGISTER_FUNCTION(atan, eval_func(atan_eval).
                         evalf_func(atan_evalf).
                         derivative_func(atan_deriv).
-                        series_func(atan_series));
+                        series_func(atan_series).
+                        latex_name("\\arctan"));
 
 //////////
 // inverse tangent (atan2(y,x))
 //////////
 
-static ex atan2_evalf(const ex & y, const ex & x)
+static ex atan2_evalf(const ex &y, const ex &x)
 {
-       BEGIN_TYPECHECK
-               TYPECHECK(y,numeric)
-               TYPECHECK(x,numeric)
-       END_TYPECHECK(atan2(y,x))
+       if (is_exactly_a<numeric>(y) && is_exactly_a<numeric>(x))
+               return atan2(ex_to<numeric>(y), ex_to<numeric>(x));
        
-       return atan(ex_to_numeric(y),ex_to_numeric(x)); // -> numeric atan(numeric)
+       return atan2(y, x).hold();
 }
 
 static ex atan2_eval(const ex & y, const ex & x)
@@ -683,11 +706,10 @@ REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
 
 static ex sinh_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(sinh(x))
+       if (is_exactly_a<numeric>(x))
+               return sinh(ex_to<numeric>(x));
        
-       return sinh(ex_to_numeric(x)); // -> numeric sinh(numeric)
+       return sinh(x).hold();
 }
 
 static ex sinh_eval(const ex & x)
@@ -696,14 +718,14 @@ static ex sinh_eval(const ex & x)
                if (x.is_zero())  // sinh(0) -> 0
                        return _ex0();        
                if (!x.info(info_flags::crational))  // sinh(float) -> float
-                       return sinh_evalf(x);
+                       return sinh(ex_to<numeric>(x));
        }
        
        if ((x/Pi).info(info_flags::numeric) &&
-               ex_to_numeric(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
+               ex_to<numeric>(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
                return I*sin(x/I);
        
-       if (is_ex_exactly_of_type(x, function)) {
+       if (is_exactly_a<function>(x)) {
                ex t = x.op(0);
                // sinh(asinh(x)) -> x
                if (is_ex_the_function(x, asinh))
@@ -729,7 +751,8 @@ static ex sinh_deriv(const ex & x, unsigned deriv_param)
 
 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
                         evalf_func(sinh_evalf).
-                        derivative_func(sinh_deriv));
+                        derivative_func(sinh_deriv).
+                        latex_name("\\sinh"));
 
 //////////
 // hyperbolic cosine (trigonometric function)
@@ -737,11 +760,10 @@ REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
 
 static ex cosh_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(cosh(x))
+       if (is_exactly_a<numeric>(x))
+               return cosh(ex_to<numeric>(x));
        
-       return cosh(ex_to_numeric(x)); // -> numeric cosh(numeric)
+       return cosh(x).hold();
 }
 
 static ex cosh_eval(const ex & x)
@@ -750,14 +772,14 @@ static ex cosh_eval(const ex & x)
                if (x.is_zero())  // cosh(0) -> 1
                        return _ex1();
                if (!x.info(info_flags::crational))  // cosh(float) -> float
-                       return cosh_evalf(x);
+                       return cosh(ex_to<numeric>(x));
        }
        
        if ((x/Pi).info(info_flags::numeric) &&
-               ex_to_numeric(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
+               ex_to<numeric>(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
                return cos(x/I);
        
-       if (is_ex_exactly_of_type(x, function)) {
+       if (is_exactly_a<function>(x)) {
                ex t = x.op(0);
                // cosh(acosh(x)) -> x
                if (is_ex_the_function(x, acosh))
@@ -783,8 +805,8 @@ static ex cosh_deriv(const ex & x, unsigned deriv_param)
 
 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
                         evalf_func(cosh_evalf).
-                        derivative_func(cosh_deriv));
-
+                        derivative_func(cosh_deriv).
+                        latex_name("\\cosh"));
 
 //////////
 // hyperbolic tangent (trigonometric function)
@@ -792,11 +814,10 @@ REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
 
 static ex tanh_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(tanh(x))
+       if (is_exactly_a<numeric>(x))
+               return tanh(ex_to<numeric>(x));
        
-       return tanh(ex_to_numeric(x)); // -> numeric tanh(numeric)
+       return tanh(x).hold();
 }
 
 static ex tanh_eval(const ex & x)
@@ -805,14 +826,14 @@ static ex tanh_eval(const ex & x)
                if (x.is_zero())  // tanh(0) -> 0
                        return _ex0();
                if (!x.info(info_flags::crational))  // tanh(float) -> float
-                       return tanh_evalf(x);
+                       return tanh(ex_to<numeric>(x));
        }
        
        if ((x/Pi).info(info_flags::numeric) &&
-               ex_to_numeric(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
+               ex_to<numeric>(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
                return I*tan(x/I);
        
-       if (is_ex_exactly_of_type(x, function)) {
+       if (is_exactly_a<function>(x)) {
                ex t = x.op(0);
                // tanh(atanh(x)) -> x
                if (is_ex_the_function(x, atanh))
@@ -855,7 +876,8 @@ static ex tanh_series(const ex &x,
 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
                         evalf_func(tanh_evalf).
                         derivative_func(tanh_deriv).
-                        series_func(tanh_series));
+                        series_func(tanh_series).
+                        latex_name("\\tanh"));
 
 //////////
 // inverse hyperbolic sine (trigonometric function)
@@ -863,11 +885,10 @@ REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
 
 static ex asinh_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(asinh(x))
+       if (is_exactly_a<numeric>(x))
+               return asinh(ex_to<numeric>(x));
        
-       return asinh(ex_to_numeric(x)); // -> numeric asinh(numeric)
+       return asinh(x).hold();
 }
 
 static ex asinh_eval(const ex & x)
@@ -878,7 +899,7 @@ static ex asinh_eval(const ex & x)
                        return _ex0();
                // asinh(float) -> float
                if (!x.info(info_flags::crational))
-                       return asinh_evalf(x);
+                       return asinh(ex_to<numeric>(x));
        }
        
        return asinh(x).hold();
@@ -902,11 +923,10 @@ REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
 
 static ex acosh_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(acosh(x))
+       if (is_exactly_a<numeric>(x))
+               return acosh(ex_to<numeric>(x));
        
-       return acosh(ex_to_numeric(x)); // -> numeric acosh(numeric)
+       return acosh(x).hold();
 }
 
 static ex acosh_eval(const ex & x)
@@ -923,7 +943,7 @@ static ex acosh_eval(const ex & x)
                        return Pi*I;
                // acosh(float) -> float
                if (!x.info(info_flags::crational))
-                       return acosh_evalf(x);
+                       return acosh(ex_to<numeric>(x));
        }
        
        return acosh(x).hold();
@@ -947,11 +967,10 @@ REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
 
 static ex atanh_evalf(const ex & x)
 {
-       BEGIN_TYPECHECK
-          TYPECHECK(x,numeric)
-       END_TYPECHECK(atanh(x))
+       if (is_exactly_a<numeric>(x))
+               return atanh(ex_to<numeric>(x));
        
-       return atanh(ex_to_numeric(x)); // -> numeric atanh(numeric)
+       return atanh(x).hold();
 }
 
 static ex atanh_eval(const ex & x)
@@ -965,7 +984,7 @@ static ex atanh_eval(const ex & x)
                        throw (pole_error("atanh_eval(): logarithmic pole",0));
                // atanh(float) -> float
                if (!x.info(info_flags::crational))
-                       return atanh_evalf(x);
+                       return atanh(ex_to<numeric>(x));
        }
        
        return atanh(x).hold();
@@ -979,27 +998,47 @@ static ex atanh_deriv(const ex & x, unsigned deriv_param)
        return power(_ex1()-power(x,_ex2()),_ex_1());
 }
 
-static ex atanh_series(const ex &x,
+static ex atanh_series(const ex &arg,
                        const relational &rel,
                        int order,
                        unsigned options)
 {
        GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
        // method:
-       // Taylor series where there is no pole or cut falls back to atan_deriv.
+       // Taylor series where there is no pole or cut falls back to atanh_deriv.
        // There are two branch cuts, one runnig from 1 up the real axis and one
        // one running from -1 down the real axis.  The points 1 and -1 are poles
        // On the branch cuts and the poles series expand
-       //     log((1+x)/(1-x))/(2*I)
+       //     (log(1+x)-log(1-x))/2
        // instead.
-       // (The constant term on the cut itself could be made simpler.)
-       const ex x_pt = x.subs(rel);
-       if (!(x_pt).info(info_flags::real))
+       const ex arg_pt = arg.subs(rel);
+       if (!(arg_pt).info(info_flags::real))
                throw do_taylor();     // Im(x) != 0
-       if ((x_pt).info(info_flags::real) && abs(x_pt)<_ex1())
+       if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1())
                throw do_taylor();     // Im(x) == 0, but abs(x)<1
-       // if we got here we have to care for cuts and poles
-       return (log((1+x)/(1-x))/2).series(rel, order, options);
+       // care for the poles, using the defining formula for atanh()...
+       if (arg_pt.is_equal(_ex1()) || arg_pt.is_equal(_ex_1()))
+               return ((log(_ex1()+arg)-log(_ex1()-arg))*_ex1_2()).series(rel, order, options);
+       // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
+       if (!(options & series_options::suppress_branchcut)) {
+               // method:
+               // This is the branch cut: assemble the primitive series manually and
+               // then add the corresponding complex step function.
+               const symbol *s = static_cast<symbol *>(rel.lhs().bp);
+               const ex point = rel.rhs();
+               const symbol foo;
+               const ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point);
+               ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2();
+               if (arg_pt<_ex0())
+                       Order0correction += log((arg_pt+_ex_1())/(arg_pt+_ex1()))*_ex1_2();
+               else
+                       Order0correction += log((arg_pt+_ex1())/(arg_pt+_ex_1()))*_ex_1_2();
+               epvector seq;
+               seq.push_back(expair(Order0correction, _ex0()));
+               seq.push_back(expair(Order(_ex1()), order));
+               return series(replarg - pseries(rel, seq), rel, order);
+       }
+       throw do_taylor();
 }
 
 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
@@ -1008,6 +1047,4 @@ REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
                          series_func(atanh_series));
 
 
-#ifndef NO_NAMESPACE_GINAC
 } // namespace GiNaC
-#endif // ndef NO_NAMESPACE_GINAC