* functions. */
/*
- * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
}
REGISTER_FUNCTION(exp, eval_func(exp_eval).
- evalf_func(exp_evalf).
- derivative_func(exp_deriv));
+ evalf_func(exp_evalf).
+ derivative_func(exp_deriv));
//////////
// natural logarithm
}
static ex log_series(const ex &arg,
- const relational &rel,
- int order,
- unsigned options)
+ const relational &rel,
+ int order,
+ unsigned options)
{
GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
ex arg_pt;
} 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;
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);
+ // (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);
+ 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));
}
REGISTER_FUNCTION(log, eval_func(log_eval).
- evalf_func(log_evalf).
- derivative_func(log_deriv).
- series_func(log_series));
+ evalf_func(log_evalf).
+ derivative_func(log_deriv).
+ series_func(log_series));
//////////
// sine (trigonometric function)
}
REGISTER_FUNCTION(sin, eval_func(sin_eval).
- evalf_func(sin_evalf).
- derivative_func(sin_deriv));
+ evalf_func(sin_evalf).
+ derivative_func(sin_deriv));
//////////
// cosine (trigonometric function)
}
REGISTER_FUNCTION(cos, eval_func(cos_eval).
- evalf_func(cos_evalf).
- derivative_func(cos_deriv));
+ evalf_func(cos_evalf).
+ derivative_func(cos_deriv));
//////////
// tangent (trigonometric function)
}
static ex tan_series(const ex &x,
- const relational &rel,
- int order,
- unsigned options)
+ 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 falls back to tan_deriv.
// On a pole simply expand sin(x)/cos(x).
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(rel, order+2);
+ return (sin(x)/cos(x)).series(rel, order+2, options);
}
REGISTER_FUNCTION(tan, eval_func(tan_eval).
- evalf_func(tan_evalf).
- derivative_func(tan_deriv).
- series_func(tan_series));
+ evalf_func(tan_evalf).
+ derivative_func(tan_deriv).
+ series_func(tan_series));
//////////
// inverse sine (arc sine)
}
REGISTER_FUNCTION(asin, eval_func(asin_eval).
- evalf_func(asin_evalf).
- derivative_func(asin_deriv));
+ evalf_func(asin_evalf).
+ derivative_func(asin_deriv));
//////////
// inverse cosine (arc cosine)
}
REGISTER_FUNCTION(acos, eval_func(acos_eval).
- evalf_func(acos_evalf).
- derivative_func(acos_deriv));
+ evalf_func(acos_evalf).
+ derivative_func(acos_deriv));
//////////
// inverse tangent (arc tangent)
// atan(0) -> 0
if (x.is_equal(_ex0()))
return _ex0();
+ // atan(1) -> Pi/4
+ if (x.is_equal(_ex1()))
+ return _ex1_4()*Pi;
+ // atan(-1) -> -Pi/4
+ if (x.is_equal(_ex_1()))
+ return _ex_1_4()*Pi;
+ if (x.is_equal(I) || x.is_equal(-I))
+ throw (pole_error("atan_eval(): logarithmic pole",0));
// atan(float) -> float
if (!x.info(info_flags::crational))
return atan_evalf(x);
}
return atan(x).hold();
-}
+}
static ex atan_deriv(const ex & x, unsigned deriv_param)
{
return power(_ex1()+power(x,_ex2()), _ex_1());
}
+static ex atan_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.
+ // There are two branch cuts, one runnig from I up the imaginary axis and
+ // 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)-log(1-I*x))/(2*I)
+ // instead.
+ const ex arg_pt = arg.subs(rel);
+ if (!(I*arg_pt).info(info_flags::real))
+ throw do_taylor(); // Re(x) != 0
+ if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1())
+ throw do_taylor(); // Re(x) == 0, but abs(x)<1
+ // 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;
+ 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));
+ evalf_func(atan_evalf).
+ derivative_func(atan_deriv).
+ series_func(atan_series));
//////////
// inverse tangent (atan2(y,x))
}
REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
- evalf_func(atan2_evalf).
- derivative_func(atan2_deriv));
+ evalf_func(atan2_evalf).
+ derivative_func(atan2_deriv));
//////////
// hyperbolic sine (trigonometric function)
}
REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
- evalf_func(sinh_evalf).
- derivative_func(sinh_deriv));
+ evalf_func(sinh_evalf).
+ derivative_func(sinh_deriv));
//////////
// hyperbolic cosine (trigonometric function)
}
REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
- evalf_func(cosh_evalf).
- derivative_func(cosh_deriv));
+ evalf_func(cosh_evalf).
+ derivative_func(cosh_deriv));
//////////
}
static ex tanh_series(const ex &x,
- const relational &rel,
- int order,
- unsigned options)
+ 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 falls back to tanh_deriv.
// On a pole simply expand sinh(x)/cosh(x).
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(rel, order+2);
+ return (sinh(x)/cosh(x)).series(rel, order+2, options);
}
REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
- evalf_func(tanh_evalf).
- derivative_func(tanh_deriv).
- series_func(tanh_series));
+ evalf_func(tanh_evalf).
+ derivative_func(tanh_deriv).
+ series_func(tanh_series));
//////////
// inverse hyperbolic sine (trigonometric function)
}
REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
- evalf_func(asinh_evalf).
- derivative_func(asinh_deriv));
+ evalf_func(asinh_evalf).
+ derivative_func(asinh_deriv));
//////////
// inverse hyperbolic cosine (trigonometric function)
}
REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
- evalf_func(acosh_evalf).
- derivative_func(acosh_deriv));
+ evalf_func(acosh_evalf).
+ derivative_func(acosh_deriv));
//////////
// inverse hyperbolic tangent (trigonometric function)
return power(_ex1()-power(x,_ex2()),_ex_1());
}
+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 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 = static_cast<symbol *>(rel.lhs().bp);
+ const ex point = rel.rhs();
+ const symbol foo;
+ 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).
- evalf_func(atanh_evalf).
- derivative_func(atanh_deriv));
+ evalf_func(atanh_evalf).
+ derivative_func(atanh_deriv).
+ series_func(atanh_series));
+
#ifndef NO_NAMESPACE_GINAC
} // namespace GiNaC
git://www.ginac.de / ginac.git / blobdiff