+static ex atanh_series(const ex &arg,
+ const relational &rel,
+ int order,
+ unsigned options)
+{
+ GINAC_ASSERT(is_exactly_a<symbol>(rel.lhs()));
+ // method:
+ // 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)-log(1-x))/2
+ // instead.
+ const ex arg_pt = arg.subs(rel);
+ if (!(arg_pt).info(info_flags::real))
+ throw do_taylor(); // Im(x) != 0
+ if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1())
+ throw do_taylor(); // Im(x) == 0, but abs(x)<1
+ // 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 = ex_to<symbol>(rel.lhs());
+ 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();
+}
+