1 /** @file inifcns_trans.cpp
3 * Implementation of transcendental (and trigonometric and hyperbolic)
7 * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
32 #include "relational.h"
40 // exponential function
43 static ex exp_evalf(const ex & x)
49 return exp(ex_to<numeric>(x)); // -> numeric exp(numeric)
52 static ex exp_eval(const ex & x)
58 // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
59 const ex TwoExOverPiI=(_ex2()*x)/(Pi*I);
60 if (TwoExOverPiI.info(info_flags::integer)) {
61 numeric z=mod(ex_to<numeric>(TwoExOverPiI),_num4());
62 if (z.is_equal(_num0()))
64 if (z.is_equal(_num1()))
66 if (z.is_equal(_num2()))
68 if (z.is_equal(_num3()))
72 if (is_ex_the_function(x, log))
76 if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
82 static ex exp_deriv(const ex & x, unsigned deriv_param)
84 GINAC_ASSERT(deriv_param==0);
86 // d/dx exp(x) -> exp(x)
90 REGISTER_FUNCTION(exp, eval_func(exp_eval).
91 evalf_func(exp_evalf).
92 derivative_func(exp_deriv).
99 static ex log_evalf(const ex & x)
103 END_TYPECHECK(log(x))
105 return log(ex_to<numeric>(x)); // -> numeric log(numeric)
108 static ex log_eval(const ex & x)
110 if (x.info(info_flags::numeric)) {
111 if (x.is_zero()) // log(0) -> infinity
112 throw(pole_error("log_eval(): log(0)",0));
113 if (x.info(info_flags::real) && x.info(info_flags::negative))
114 return (log(-x)+I*Pi);
115 if (x.is_equal(_ex1())) // log(1) -> 0
117 if (x.is_equal(I)) // log(I) -> Pi*I/2
118 return (Pi*I*_num1_2());
119 if (x.is_equal(-I)) // log(-I) -> -Pi*I/2
120 return (Pi*I*_num_1_2());
122 if (!x.info(info_flags::crational))
125 // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
126 if (is_ex_the_function(x, exp)) {
128 if (t.info(info_flags::numeric)) {
129 numeric nt = ex_to<numeric>(t);
135 return log(x).hold();
138 static ex log_deriv(const ex & x, unsigned deriv_param)
140 GINAC_ASSERT(deriv_param==0);
142 // d/dx log(x) -> 1/x
143 return power(x, _ex_1());
146 static ex log_series(const ex &arg,
147 const relational &rel,
151 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
153 bool must_expand_arg = false;
154 // maybe substitution of rel into arg fails because of a pole
156 arg_pt = arg.subs(rel);
157 } catch (pole_error) {
158 must_expand_arg = true;
160 // or we are at the branch point anyways
161 if (arg_pt.is_zero())
162 must_expand_arg = true;
164 if (must_expand_arg) {
166 // This is the branch point: Series expand the argument first, then
167 // trivially factorize it to isolate that part which has constant
168 // leading coefficient in this fashion:
169 // x^n + Order(x^(n+m)) -> x^n * (1 + Order(x^m)).
170 // Return a plain n*log(x) for the x^n part and series expand the
171 // other part. Add them together and reexpand again in order to have
172 // one unnested pseries object. All this also works for negative n.
173 const pseries argser = ex_to<pseries>(arg.series(rel, order, options));
174 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
175 const ex point = rel.rhs();
176 const int n = argser.ldegree(*s);
178 // construct what we carelessly called the n*log(x) term above
179 ex coeff = argser.coeff(*s, n);
180 // expand the log, but only if coeff is real and > 0, since otherwise
181 // it would make the branch cut run into the wrong direction
182 if (coeff.info(info_flags::positive))
183 seq.push_back(expair(n*log(*s-point)+log(coeff), _ex0()));
185 seq.push_back(expair(log(coeff*pow(*s-point, n)), _ex0()));
186 if (!argser.is_terminating() || argser.nops()!=1) {
187 // in this case n more terms are needed
188 // (sadly, to generate them, we have to start from the beginning)
189 ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
190 return pseries(rel, seq).add_series(ex_to<pseries>(log(newarg).series(rel, order, options)));
191 } else // it was a monomial
192 return pseries(rel, seq);
194 if (!(options & series_options::suppress_branchcut) &&
195 arg_pt.info(info_flags::negative)) {
197 // This is the branch cut: assemble the primitive series manually and
198 // then add the corresponding complex step function.
199 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
200 const ex point = rel.rhs();
202 const ex replarg = series(log(arg), *s==foo, order).subs(foo==point);
204 seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
205 seq.push_back(expair(Order(_ex1()), order));
206 return series(replarg - I*Pi + pseries(rel, seq), rel, order);
208 throw do_taylor(); // caught by function::series()
211 REGISTER_FUNCTION(log, eval_func(log_eval).
212 evalf_func(log_evalf).
213 derivative_func(log_deriv).
214 series_func(log_series).
218 // sine (trigonometric function)
221 static ex sin_evalf(const ex & x)
225 END_TYPECHECK(sin(x))
227 return sin(ex_to<numeric>(x)); // -> numeric sin(numeric)
230 static ex sin_eval(const ex & x)
232 // sin(n/d*Pi) -> { all known non-nested radicals }
233 const ex SixtyExOverPi = _ex60()*x/Pi;
235 if (SixtyExOverPi.info(info_flags::integer)) {
236 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
238 // wrap to interval [0, Pi)
243 // wrap to interval [0, Pi/2)
246 if (z.is_equal(_num0())) // sin(0) -> 0
248 if (z.is_equal(_num5())) // sin(Pi/12) -> sqrt(6)/4*(1-sqrt(3)/3)
249 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
250 if (z.is_equal(_num6())) // sin(Pi/10) -> sqrt(5)/4-1/4
251 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
252 if (z.is_equal(_num10())) // sin(Pi/6) -> 1/2
253 return sign*_ex1_2();
254 if (z.is_equal(_num15())) // sin(Pi/4) -> sqrt(2)/2
255 return sign*_ex1_2()*power(_ex2(),_ex1_2());
256 if (z.is_equal(_num18())) // sin(3/10*Pi) -> sqrt(5)/4+1/4
257 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
258 if (z.is_equal(_num20())) // sin(Pi/3) -> sqrt(3)/2
259 return sign*_ex1_2()*power(_ex3(),_ex1_2());
260 if (z.is_equal(_num25())) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
261 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
262 if (z.is_equal(_num30())) // sin(Pi/2) -> 1
266 if (is_exactly_a<function>(x)) {
269 if (is_ex_the_function(x, asin))
271 // sin(acos(x)) -> sqrt(1-x^2)
272 if (is_ex_the_function(x, acos))
273 return power(_ex1()-power(t,_ex2()),_ex1_2());
274 // sin(atan(x)) -> x*(1+x^2)^(-1/2)
275 if (is_ex_the_function(x, atan))
276 return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
279 // sin(float) -> float
280 if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
283 return sin(x).hold();
286 static ex sin_deriv(const ex & x, unsigned deriv_param)
288 GINAC_ASSERT(deriv_param==0);
290 // d/dx sin(x) -> cos(x)
294 REGISTER_FUNCTION(sin, eval_func(sin_eval).
295 evalf_func(sin_evalf).
296 derivative_func(sin_deriv).
297 latex_name("\\sin"));
300 // cosine (trigonometric function)
303 static ex cos_evalf(const ex & x)
307 END_TYPECHECK(cos(x))
309 return cos(ex_to<numeric>(x)); // -> numeric cos(numeric)
312 static ex cos_eval(const ex & x)
314 // cos(n/d*Pi) -> { all known non-nested radicals }
315 const ex SixtyExOverPi = _ex60()*x/Pi;
317 if (SixtyExOverPi.info(info_flags::integer)) {
318 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
320 // wrap to interval [0, Pi)
324 // wrap to interval [0, Pi/2)
328 if (z.is_equal(_num0())) // cos(0) -> 1
330 if (z.is_equal(_num5())) // cos(Pi/12) -> sqrt(6)/4*(1+sqrt(3)/3)
331 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
332 if (z.is_equal(_num10())) // cos(Pi/6) -> sqrt(3)/2
333 return sign*_ex1_2()*power(_ex3(),_ex1_2());
334 if (z.is_equal(_num12())) // cos(Pi/5) -> sqrt(5)/4+1/4
335 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
336 if (z.is_equal(_num15())) // cos(Pi/4) -> sqrt(2)/2
337 return sign*_ex1_2()*power(_ex2(),_ex1_2());
338 if (z.is_equal(_num20())) // cos(Pi/3) -> 1/2
339 return sign*_ex1_2();
340 if (z.is_equal(_num24())) // cos(2/5*Pi) -> sqrt(5)/4-1/4x
341 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
342 if (z.is_equal(_num25())) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
343 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
344 if (z.is_equal(_num30())) // cos(Pi/2) -> 0
348 if (is_exactly_a<function>(x)) {
351 if (is_ex_the_function(x, acos))
353 // cos(asin(x)) -> (1-x^2)^(1/2)
354 if (is_ex_the_function(x, asin))
355 return power(_ex1()-power(t,_ex2()),_ex1_2());
356 // cos(atan(x)) -> (1+x^2)^(-1/2)
357 if (is_ex_the_function(x, atan))
358 return power(_ex1()+power(t,_ex2()),_ex_1_2());
361 // cos(float) -> float
362 if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
365 return cos(x).hold();
368 static ex cos_deriv(const ex & x, unsigned deriv_param)
370 GINAC_ASSERT(deriv_param==0);
372 // d/dx cos(x) -> -sin(x)
373 return _ex_1()*sin(x);
376 REGISTER_FUNCTION(cos, eval_func(cos_eval).
377 evalf_func(cos_evalf).
378 derivative_func(cos_deriv).
379 latex_name("\\cos"));
382 // tangent (trigonometric function)
385 static ex tan_evalf(const ex & x)
389 END_TYPECHECK(tan(x)) // -> numeric tan(numeric)
391 return tan(ex_to<numeric>(x));
394 static ex tan_eval(const ex & x)
396 // tan(n/d*Pi) -> { all known non-nested radicals }
397 const ex SixtyExOverPi = _ex60()*x/Pi;
399 if (SixtyExOverPi.info(info_flags::integer)) {
400 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num60());
402 // wrap to interval [0, Pi)
406 // wrap to interval [0, Pi/2)
410 if (z.is_equal(_num0())) // tan(0) -> 0
412 if (z.is_equal(_num5())) // tan(Pi/12) -> 2-sqrt(3)
413 return sign*(_ex2()-power(_ex3(),_ex1_2()));
414 if (z.is_equal(_num10())) // tan(Pi/6) -> sqrt(3)/3
415 return sign*_ex1_3()*power(_ex3(),_ex1_2());
416 if (z.is_equal(_num15())) // tan(Pi/4) -> 1
418 if (z.is_equal(_num20())) // tan(Pi/3) -> sqrt(3)
419 return sign*power(_ex3(),_ex1_2());
420 if (z.is_equal(_num25())) // tan(5/12*Pi) -> 2+sqrt(3)
421 return sign*(power(_ex3(),_ex1_2())+_ex2());
422 if (z.is_equal(_num30())) // tan(Pi/2) -> infinity
423 throw (pole_error("tan_eval(): simple pole",1));
426 if (is_exactly_a<function>(x)) {
429 if (is_ex_the_function(x, atan))
431 // tan(asin(x)) -> x*(1+x^2)^(-1/2)
432 if (is_ex_the_function(x, asin))
433 return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
434 // tan(acos(x)) -> (1-x^2)^(1/2)/x
435 if (is_ex_the_function(x, acos))
436 return power(t,_ex_1())*power(_ex1()-power(t,_ex2()),_ex1_2());
439 // tan(float) -> float
440 if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
444 return tan(x).hold();
447 static ex tan_deriv(const ex & x, unsigned deriv_param)
449 GINAC_ASSERT(deriv_param==0);
451 // d/dx tan(x) -> 1+tan(x)^2;
452 return (_ex1()+power(tan(x),_ex2()));
455 static ex tan_series(const ex &x,
456 const relational &rel,
460 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
462 // Taylor series where there is no pole falls back to tan_deriv.
463 // On a pole simply expand sin(x)/cos(x).
464 const ex x_pt = x.subs(rel);
465 if (!(2*x_pt/Pi).info(info_flags::odd))
466 throw do_taylor(); // caught by function::series()
467 // if we got here we have to care for a simple pole
468 return (sin(x)/cos(x)).series(rel, order+2, options);
471 REGISTER_FUNCTION(tan, eval_func(tan_eval).
472 evalf_func(tan_evalf).
473 derivative_func(tan_deriv).
474 series_func(tan_series).
475 latex_name("\\tan"));
478 // inverse sine (arc sine)
481 static ex asin_evalf(const ex & x)
485 END_TYPECHECK(asin(x))
487 return asin(ex_to<numeric>(x)); // -> numeric asin(numeric)
490 static ex asin_eval(const ex & x)
492 if (x.info(info_flags::numeric)) {
497 if (x.is_equal(_ex1_2()))
498 return numeric(1,6)*Pi;
500 if (x.is_equal(_ex1()))
502 // asin(-1/2) -> -Pi/6
503 if (x.is_equal(_ex_1_2()))
504 return numeric(-1,6)*Pi;
506 if (x.is_equal(_ex_1()))
507 return _num_1_2()*Pi;
508 // asin(float) -> float
509 if (!x.info(info_flags::crational))
510 return asin_evalf(x);
513 return asin(x).hold();
516 static ex asin_deriv(const ex & x, unsigned deriv_param)
518 GINAC_ASSERT(deriv_param==0);
520 // d/dx asin(x) -> 1/sqrt(1-x^2)
521 return power(1-power(x,_ex2()),_ex_1_2());
524 REGISTER_FUNCTION(asin, eval_func(asin_eval).
525 evalf_func(asin_evalf).
526 derivative_func(asin_deriv).
527 latex_name("\\arcsin"));
530 // inverse cosine (arc cosine)
533 static ex acos_evalf(const ex & x)
537 END_TYPECHECK(acos(x))
539 return acos(ex_to<numeric>(x)); // -> numeric acos(numeric)
542 static ex acos_eval(const ex & x)
544 if (x.info(info_flags::numeric)) {
546 if (x.is_equal(_ex1()))
549 if (x.is_equal(_ex1_2()))
554 // acos(-1/2) -> 2/3*Pi
555 if (x.is_equal(_ex_1_2()))
556 return numeric(2,3)*Pi;
558 if (x.is_equal(_ex_1()))
560 // acos(float) -> float
561 if (!x.info(info_flags::crational))
562 return acos_evalf(x);
565 return acos(x).hold();
568 static ex acos_deriv(const ex & x, unsigned deriv_param)
570 GINAC_ASSERT(deriv_param==0);
572 // d/dx acos(x) -> -1/sqrt(1-x^2)
573 return _ex_1()*power(1-power(x,_ex2()),_ex_1_2());
576 REGISTER_FUNCTION(acos, eval_func(acos_eval).
577 evalf_func(acos_evalf).
578 derivative_func(acos_deriv).
579 latex_name("\\arccos"));
582 // inverse tangent (arc tangent)
585 static ex atan_evalf(const ex & x)
589 END_TYPECHECK(atan(x))
591 return atan(ex_to<numeric>(x)); // -> numeric atan(numeric)
594 static ex atan_eval(const ex & x)
596 if (x.info(info_flags::numeric)) {
601 if (x.is_equal(_ex1()))
604 if (x.is_equal(_ex_1()))
606 if (x.is_equal(I) || x.is_equal(-I))
607 throw (pole_error("atan_eval(): logarithmic pole",0));
608 // atan(float) -> float
609 if (!x.info(info_flags::crational))
610 return atan_evalf(x);
613 return atan(x).hold();
616 static ex atan_deriv(const ex & x, unsigned deriv_param)
618 GINAC_ASSERT(deriv_param==0);
620 // d/dx atan(x) -> 1/(1+x^2)
621 return power(_ex1()+power(x,_ex2()), _ex_1());
624 static ex atan_series(const ex &arg,
625 const relational &rel,
629 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
631 // Taylor series where there is no pole or cut falls back to atan_deriv.
632 // There are two branch cuts, one runnig from I up the imaginary axis and
633 // one running from -I down the imaginary axis. The points I and -I are
635 // On the branch cuts and the poles series expand
636 // (log(1+I*x)-log(1-I*x))/(2*I)
638 const ex arg_pt = arg.subs(rel);
639 if (!(I*arg_pt).info(info_flags::real))
640 throw do_taylor(); // Re(x) != 0
641 if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1())
642 throw do_taylor(); // Re(x) == 0, but abs(x)<1
643 // care for the poles, using the defining formula for atan()...
644 if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
645 return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
646 if (!(options & series_options::suppress_branchcut)) {
648 // This is the branch cut: assemble the primitive series manually and
649 // then add the corresponding complex step function.
650 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
651 const ex point = rel.rhs();
653 const ex replarg = series(atan(arg), *s==foo, order).subs(foo==point);
654 ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2();
655 if ((I*arg_pt)<_ex0())
656 Order0correction += log((I*arg_pt+_ex_1())/(I*arg_pt+_ex1()))*I*_ex_1_2();
658 Order0correction += log((I*arg_pt+_ex1())/(I*arg_pt+_ex_1()))*I*_ex1_2();
660 seq.push_back(expair(Order0correction, _ex0()));
661 seq.push_back(expair(Order(_ex1()), order));
662 return series(replarg - pseries(rel, seq), rel, order);
667 REGISTER_FUNCTION(atan, eval_func(atan_eval).
668 evalf_func(atan_evalf).
669 derivative_func(atan_deriv).
670 series_func(atan_series).
671 latex_name("\\arctan"));
674 // inverse tangent (atan2(y,x))
677 static ex atan2_evalf(const ex & y, const ex & x)
682 END_TYPECHECK(atan2(y,x))
684 return atan(ex_to<numeric>(y),ex_to<numeric>(x)); // -> numeric atan(numeric)
687 static ex atan2_eval(const ex & y, const ex & x)
689 if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
690 x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
691 return atan2_evalf(y,x);
694 return atan2(y,x).hold();
697 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
699 GINAC_ASSERT(deriv_param<2);
701 if (deriv_param==0) {
703 return x*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
706 return -y*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
709 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
710 evalf_func(atan2_evalf).
711 derivative_func(atan2_deriv));
714 // hyperbolic sine (trigonometric function)
717 static ex sinh_evalf(const ex & x)
721 END_TYPECHECK(sinh(x))
723 return sinh(ex_to<numeric>(x)); // -> numeric sinh(numeric)
726 static ex sinh_eval(const ex & x)
728 if (x.info(info_flags::numeric)) {
729 if (x.is_zero()) // sinh(0) -> 0
731 if (!x.info(info_flags::crational)) // sinh(float) -> float
732 return sinh_evalf(x);
735 if ((x/Pi).info(info_flags::numeric) &&
736 ex_to<numeric>(x/Pi).real().is_zero()) // sinh(I*x) -> I*sin(x)
739 if (is_exactly_a<function>(x)) {
741 // sinh(asinh(x)) -> x
742 if (is_ex_the_function(x, asinh))
744 // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2)
745 if (is_ex_the_function(x, acosh))
746 return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2());
747 // sinh(atanh(x)) -> x*(1-x^2)^(-1/2)
748 if (is_ex_the_function(x, atanh))
749 return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
752 return sinh(x).hold();
755 static ex sinh_deriv(const ex & x, unsigned deriv_param)
757 GINAC_ASSERT(deriv_param==0);
759 // d/dx sinh(x) -> cosh(x)
763 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
764 evalf_func(sinh_evalf).
765 derivative_func(sinh_deriv).
766 latex_name("\\sinh"));
769 // hyperbolic cosine (trigonometric function)
772 static ex cosh_evalf(const ex & x)
776 END_TYPECHECK(cosh(x))
778 return cosh(ex_to<numeric>(x)); // -> numeric cosh(numeric)
781 static ex cosh_eval(const ex & x)
783 if (x.info(info_flags::numeric)) {
784 if (x.is_zero()) // cosh(0) -> 1
786 if (!x.info(info_flags::crational)) // cosh(float) -> float
787 return cosh_evalf(x);
790 if ((x/Pi).info(info_flags::numeric) &&
791 ex_to<numeric>(x/Pi).real().is_zero()) // cosh(I*x) -> cos(x)
794 if (is_exactly_a<function>(x)) {
796 // cosh(acosh(x)) -> x
797 if (is_ex_the_function(x, acosh))
799 // cosh(asinh(x)) -> (1+x^2)^(1/2)
800 if (is_ex_the_function(x, asinh))
801 return power(_ex1()+power(t,_ex2()),_ex1_2());
802 // cosh(atanh(x)) -> (1-x^2)^(-1/2)
803 if (is_ex_the_function(x, atanh))
804 return power(_ex1()-power(t,_ex2()),_ex_1_2());
807 return cosh(x).hold();
810 static ex cosh_deriv(const ex & x, unsigned deriv_param)
812 GINAC_ASSERT(deriv_param==0);
814 // d/dx cosh(x) -> sinh(x)
818 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
819 evalf_func(cosh_evalf).
820 derivative_func(cosh_deriv).
821 latex_name("\\cosh"));
824 // hyperbolic tangent (trigonometric function)
827 static ex tanh_evalf(const ex & x)
831 END_TYPECHECK(tanh(x))
833 return tanh(ex_to<numeric>(x)); // -> numeric tanh(numeric)
836 static ex tanh_eval(const ex & x)
838 if (x.info(info_flags::numeric)) {
839 if (x.is_zero()) // tanh(0) -> 0
841 if (!x.info(info_flags::crational)) // tanh(float) -> float
842 return tanh_evalf(x);
845 if ((x/Pi).info(info_flags::numeric) &&
846 ex_to<numeric>(x/Pi).real().is_zero()) // tanh(I*x) -> I*tan(x);
849 if (is_exactly_a<function>(x)) {
851 // tanh(atanh(x)) -> x
852 if (is_ex_the_function(x, atanh))
854 // tanh(asinh(x)) -> x*(1+x^2)^(-1/2)
855 if (is_ex_the_function(x, asinh))
856 return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
857 // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x
858 if (is_ex_the_function(x, acosh))
859 return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2())*power(t,_ex_1());
862 return tanh(x).hold();
865 static ex tanh_deriv(const ex & x, unsigned deriv_param)
867 GINAC_ASSERT(deriv_param==0);
869 // d/dx tanh(x) -> 1-tanh(x)^2
870 return _ex1()-power(tanh(x),_ex2());
873 static ex tanh_series(const ex &x,
874 const relational &rel,
878 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
880 // Taylor series where there is no pole falls back to tanh_deriv.
881 // On a pole simply expand sinh(x)/cosh(x).
882 const ex x_pt = x.subs(rel);
883 if (!(2*I*x_pt/Pi).info(info_flags::odd))
884 throw do_taylor(); // caught by function::series()
885 // if we got here we have to care for a simple pole
886 return (sinh(x)/cosh(x)).series(rel, order+2, options);
889 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
890 evalf_func(tanh_evalf).
891 derivative_func(tanh_deriv).
892 series_func(tanh_series).
893 latex_name("\\tanh"));
896 // inverse hyperbolic sine (trigonometric function)
899 static ex asinh_evalf(const ex & x)
903 END_TYPECHECK(asinh(x))
905 return asinh(ex_to<numeric>(x)); // -> numeric asinh(numeric)
908 static ex asinh_eval(const ex & x)
910 if (x.info(info_flags::numeric)) {
914 // asinh(float) -> float
915 if (!x.info(info_flags::crational))
916 return asinh_evalf(x);
919 return asinh(x).hold();
922 static ex asinh_deriv(const ex & x, unsigned deriv_param)
924 GINAC_ASSERT(deriv_param==0);
926 // d/dx asinh(x) -> 1/sqrt(1+x^2)
927 return power(_ex1()+power(x,_ex2()),_ex_1_2());
930 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
931 evalf_func(asinh_evalf).
932 derivative_func(asinh_deriv));
935 // inverse hyperbolic cosine (trigonometric function)
938 static ex acosh_evalf(const ex & x)
942 END_TYPECHECK(acosh(x))
944 return acosh(ex_to<numeric>(x)); // -> numeric acosh(numeric)
947 static ex acosh_eval(const ex & x)
949 if (x.info(info_flags::numeric)) {
950 // acosh(0) -> Pi*I/2
952 return Pi*I*numeric(1,2);
954 if (x.is_equal(_ex1()))
957 if (x.is_equal(_ex_1()))
959 // acosh(float) -> float
960 if (!x.info(info_flags::crational))
961 return acosh_evalf(x);
964 return acosh(x).hold();
967 static ex acosh_deriv(const ex & x, unsigned deriv_param)
969 GINAC_ASSERT(deriv_param==0);
971 // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
972 return power(x+_ex_1(),_ex_1_2())*power(x+_ex1(),_ex_1_2());
975 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
976 evalf_func(acosh_evalf).
977 derivative_func(acosh_deriv));
980 // inverse hyperbolic tangent (trigonometric function)
983 static ex atanh_evalf(const ex & x)
987 END_TYPECHECK(atanh(x))
989 return atanh(ex_to<numeric>(x)); // -> numeric atanh(numeric)
992 static ex atanh_eval(const ex & x)
994 if (x.info(info_flags::numeric)) {
998 // atanh({+|-}1) -> throw
999 if (x.is_equal(_ex1()) || x.is_equal(_ex_1()))
1000 throw (pole_error("atanh_eval(): logarithmic pole",0));
1001 // atanh(float) -> float
1002 if (!x.info(info_flags::crational))
1003 return atanh_evalf(x);
1006 return atanh(x).hold();
1009 static ex atanh_deriv(const ex & x, unsigned deriv_param)
1011 GINAC_ASSERT(deriv_param==0);
1013 // d/dx atanh(x) -> 1/(1-x^2)
1014 return power(_ex1()-power(x,_ex2()),_ex_1());
1017 static ex atanh_series(const ex &arg,
1018 const relational &rel,
1022 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
1024 // Taylor series where there is no pole or cut falls back to atanh_deriv.
1025 // There are two branch cuts, one runnig from 1 up the real axis and one
1026 // one running from -1 down the real axis. The points 1 and -1 are poles
1027 // On the branch cuts and the poles series expand
1028 // (log(1+x)-log(1-x))/2
1030 const ex arg_pt = arg.subs(rel);
1031 if (!(arg_pt).info(info_flags::real))
1032 throw do_taylor(); // Im(x) != 0
1033 if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1())
1034 throw do_taylor(); // Im(x) == 0, but abs(x)<1
1035 // care for the poles, using the defining formula for atanh()...
1036 if (arg_pt.is_equal(_ex1()) || arg_pt.is_equal(_ex_1()))
1037 return ((log(_ex1()+arg)-log(_ex1()-arg))*_ex1_2()).series(rel, order, options);
1038 // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
1039 if (!(options & series_options::suppress_branchcut)) {
1041 // This is the branch cut: assemble the primitive series manually and
1042 // then add the corresponding complex step function.
1043 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
1044 const ex point = rel.rhs();
1046 const ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point);
1047 ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2();
1049 Order0correction += log((arg_pt+_ex_1())/(arg_pt+_ex1()))*_ex1_2();
1051 Order0correction += log((arg_pt+_ex1())/(arg_pt+_ex_1()))*_ex_1_2();
1053 seq.push_back(expair(Order0correction, _ex0()));
1054 seq.push_back(expair(Order(_ex1()), order));
1055 return series(replarg - pseries(rel, seq), rel, order);
1060 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
1061 evalf_func(atanh_evalf).
1062 derivative_func(atanh_deriv).
1063 series_func(atanh_series));
1066 } // namespace GiNaC