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 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));
98 static ex log_evalf(const ex & x)
102 END_TYPECHECK(log(x))
104 return log(ex_to_numeric(x)); // -> numeric log(numeric)
107 static ex log_eval(const ex & x)
109 if (x.info(info_flags::numeric)) {
110 if (x.is_equal(_ex0())) // log(0) -> infinity
111 throw(pole_error("log_eval(): log(0)",0));
112 if (x.info(info_flags::real) && x.info(info_flags::negative))
113 return (log(-x)+I*Pi);
114 if (x.is_equal(_ex1())) // log(1) -> 0
116 if (x.is_equal(I)) // log(I) -> Pi*I/2
117 return (Pi*I*_num1_2());
118 if (x.is_equal(-I)) // log(-I) -> -Pi*I/2
119 return (Pi*I*_num_1_2());
121 if (!x.info(info_flags::crational))
124 // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
125 if (is_ex_the_function(x, exp)) {
127 if (t.info(info_flags::numeric)) {
128 numeric nt = ex_to_numeric(t);
134 return log(x).hold();
137 static ex log_deriv(const ex & x, unsigned deriv_param)
139 GINAC_ASSERT(deriv_param==0);
141 // d/dx log(x) -> 1/x
142 return power(x, _ex_1());
145 static ex log_series(const ex &arg,
146 const relational &rel,
150 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
152 bool must_expand_arg = false;
153 // maybe substitution of rel into arg fails because of a pole
155 arg_pt = arg.subs(rel);
156 } catch (pole_error) {
157 must_expand_arg = true;
159 // or we are at the branch point anyways
160 if (arg_pt.is_zero())
161 must_expand_arg = true;
163 if (must_expand_arg) {
165 // This is the branch point: Series expand the argument first, then
166 // trivially factorize it to isolate that part which has constant
167 // leading coefficient in this fashion:
168 // x^n + Order(x^(n+m)) -> x^n * (1 + Order(x^m)).
169 // Return a plain n*log(x) for the x^n part and series expand the
170 // other part. Add them together and reexpand again in order to have
171 // one unnested pseries object. All this also works for negative n.
172 const pseries argser = ex_to_pseries(arg.series(rel, order, options));
173 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
174 const ex point = rel.rhs();
175 const int n = argser.ldegree(*s);
177 // construct what we carelessly called the n*log(x) term above
178 ex coeff = argser.coeff(*s, n);
179 // expand the log, but only if coeff is real and > 0, since otherwise
180 // it would make the branch cut run into the wrong direction
181 if (coeff.info(info_flags::positive))
182 seq.push_back(expair(n*log(*s-point)+log(coeff), _ex0()));
184 seq.push_back(expair(log(coeff*pow(*s-point, n)), _ex0()));
185 if (!argser.is_terminating() || argser.nops()!=1) {
186 // in this case n more terms are needed
187 // (sadly, to generate them, we have to start from the beginning)
188 ex newarg = ex_to_pseries((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
189 return pseries(rel, seq).add_series(ex_to_pseries(log(newarg).series(rel, order, options)));
190 } else // it was a monomial
191 return pseries(rel, seq);
193 if (!(options & series_options::suppress_branchcut) &&
194 arg_pt.info(info_flags::negative)) {
196 // This is the branch cut: assemble the primitive series manually and
197 // then add the corresponding complex step function.
198 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
199 const ex point = rel.rhs();
201 ex replarg = series(log(arg), *s==foo, order).subs(foo==point);
203 seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
204 seq.push_back(expair(Order(_ex1()), order));
205 return series(replarg - I*Pi + pseries(rel, seq), rel, order);
207 throw do_taylor(); // caught by function::series()
210 REGISTER_FUNCTION(log, eval_func(log_eval).
211 evalf_func(log_evalf).
212 derivative_func(log_deriv).
213 series_func(log_series));
216 // sine (trigonometric function)
219 static ex sin_evalf(const ex & x)
223 END_TYPECHECK(sin(x))
225 return sin(ex_to_numeric(x)); // -> numeric sin(numeric)
228 static ex sin_eval(const ex & x)
230 // sin(n/d*Pi) -> { all known non-nested radicals }
231 ex SixtyExOverPi = _ex60()*x/Pi;
233 if (SixtyExOverPi.info(info_flags::integer)) {
234 numeric z = mod(ex_to_numeric(SixtyExOverPi),_num120());
236 // wrap to interval [0, Pi)
241 // wrap to interval [0, Pi/2)
244 if (z.is_equal(_num0())) // sin(0) -> 0
246 if (z.is_equal(_num5())) // sin(Pi/12) -> sqrt(6)/4*(1-sqrt(3)/3)
247 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
248 if (z.is_equal(_num6())) // sin(Pi/10) -> sqrt(5)/4-1/4
249 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
250 if (z.is_equal(_num10())) // sin(Pi/6) -> 1/2
251 return sign*_ex1_2();
252 if (z.is_equal(_num15())) // sin(Pi/4) -> sqrt(2)/2
253 return sign*_ex1_2()*power(_ex2(),_ex1_2());
254 if (z.is_equal(_num18())) // sin(3/10*Pi) -> sqrt(5)/4+1/4
255 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
256 if (z.is_equal(_num20())) // sin(Pi/3) -> sqrt(3)/2
257 return sign*_ex1_2()*power(_ex3(),_ex1_2());
258 if (z.is_equal(_num25())) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
259 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
260 if (z.is_equal(_num30())) // sin(Pi/2) -> 1
264 if (is_ex_exactly_of_type(x, function)) {
267 if (is_ex_the_function(x, asin))
269 // sin(acos(x)) -> sqrt(1-x^2)
270 if (is_ex_the_function(x, acos))
271 return power(_ex1()-power(t,_ex2()),_ex1_2());
272 // sin(atan(x)) -> x*(1+x^2)^(-1/2)
273 if (is_ex_the_function(x, atan))
274 return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
277 // sin(float) -> float
278 if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
281 return sin(x).hold();
284 static ex sin_deriv(const ex & x, unsigned deriv_param)
286 GINAC_ASSERT(deriv_param==0);
288 // d/dx sin(x) -> cos(x)
292 REGISTER_FUNCTION(sin, eval_func(sin_eval).
293 evalf_func(sin_evalf).
294 derivative_func(sin_deriv));
297 // cosine (trigonometric function)
300 static ex cos_evalf(const ex & x)
304 END_TYPECHECK(cos(x))
306 return cos(ex_to_numeric(x)); // -> numeric cos(numeric)
309 static ex cos_eval(const ex & x)
311 // cos(n/d*Pi) -> { all known non-nested radicals }
312 ex SixtyExOverPi = _ex60()*x/Pi;
314 if (SixtyExOverPi.info(info_flags::integer)) {
315 numeric z = mod(ex_to_numeric(SixtyExOverPi),_num120());
317 // wrap to interval [0, Pi)
321 // wrap to interval [0, Pi/2)
325 if (z.is_equal(_num0())) // cos(0) -> 1
327 if (z.is_equal(_num5())) // cos(Pi/12) -> sqrt(6)/4*(1+sqrt(3)/3)
328 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
329 if (z.is_equal(_num10())) // cos(Pi/6) -> sqrt(3)/2
330 return sign*_ex1_2()*power(_ex3(),_ex1_2());
331 if (z.is_equal(_num12())) // cos(Pi/5) -> sqrt(5)/4+1/4
332 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
333 if (z.is_equal(_num15())) // cos(Pi/4) -> sqrt(2)/2
334 return sign*_ex1_2()*power(_ex2(),_ex1_2());
335 if (z.is_equal(_num20())) // cos(Pi/3) -> 1/2
336 return sign*_ex1_2();
337 if (z.is_equal(_num24())) // cos(2/5*Pi) -> sqrt(5)/4-1/4x
338 return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
339 if (z.is_equal(_num25())) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
340 return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
341 if (z.is_equal(_num30())) // cos(Pi/2) -> 0
345 if (is_ex_exactly_of_type(x, function)) {
348 if (is_ex_the_function(x, acos))
350 // cos(asin(x)) -> (1-x^2)^(1/2)
351 if (is_ex_the_function(x, asin))
352 return power(_ex1()-power(t,_ex2()),_ex1_2());
353 // cos(atan(x)) -> (1+x^2)^(-1/2)
354 if (is_ex_the_function(x, atan))
355 return power(_ex1()+power(t,_ex2()),_ex_1_2());
358 // cos(float) -> float
359 if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
362 return cos(x).hold();
365 static ex cos_deriv(const ex & x, unsigned deriv_param)
367 GINAC_ASSERT(deriv_param==0);
369 // d/dx cos(x) -> -sin(x)
370 return _ex_1()*sin(x);
373 REGISTER_FUNCTION(cos, eval_func(cos_eval).
374 evalf_func(cos_evalf).
375 derivative_func(cos_deriv));
378 // tangent (trigonometric function)
381 static ex tan_evalf(const ex & x)
385 END_TYPECHECK(tan(x)) // -> numeric tan(numeric)
387 return tan(ex_to_numeric(x));
390 static ex tan_eval(const ex & x)
392 // tan(n/d*Pi) -> { all known non-nested radicals }
393 ex SixtyExOverPi = _ex60()*x/Pi;
395 if (SixtyExOverPi.info(info_flags::integer)) {
396 numeric z = mod(ex_to_numeric(SixtyExOverPi),_num60());
398 // wrap to interval [0, Pi)
402 // wrap to interval [0, Pi/2)
406 if (z.is_equal(_num0())) // tan(0) -> 0
408 if (z.is_equal(_num5())) // tan(Pi/12) -> 2-sqrt(3)
409 return sign*(_ex2()-power(_ex3(),_ex1_2()));
410 if (z.is_equal(_num10())) // tan(Pi/6) -> sqrt(3)/3
411 return sign*_ex1_3()*power(_ex3(),_ex1_2());
412 if (z.is_equal(_num15())) // tan(Pi/4) -> 1
414 if (z.is_equal(_num20())) // tan(Pi/3) -> sqrt(3)
415 return sign*power(_ex3(),_ex1_2());
416 if (z.is_equal(_num25())) // tan(5/12*Pi) -> 2+sqrt(3)
417 return sign*(power(_ex3(),_ex1_2())+_ex2());
418 if (z.is_equal(_num30())) // tan(Pi/2) -> infinity
419 throw (pole_error("tan_eval(): simple pole",1));
422 if (is_ex_exactly_of_type(x, function)) {
425 if (is_ex_the_function(x, atan))
427 // tan(asin(x)) -> x*(1+x^2)^(-1/2)
428 if (is_ex_the_function(x, asin))
429 return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
430 // tan(acos(x)) -> (1-x^2)^(1/2)/x
431 if (is_ex_the_function(x, acos))
432 return power(t,_ex_1())*power(_ex1()-power(t,_ex2()),_ex1_2());
435 // tan(float) -> float
436 if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
440 return tan(x).hold();
443 static ex tan_deriv(const ex & x, unsigned deriv_param)
445 GINAC_ASSERT(deriv_param==0);
447 // d/dx tan(x) -> 1+tan(x)^2;
448 return (_ex1()+power(tan(x),_ex2()));
451 static ex tan_series(const ex &x,
452 const relational &rel,
456 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
458 // Taylor series where there is no pole falls back to tan_deriv.
459 // On a pole simply expand sin(x)/cos(x).
460 const ex x_pt = x.subs(rel);
461 if (!(2*x_pt/Pi).info(info_flags::odd))
462 throw do_taylor(); // caught by function::series()
463 // if we got here we have to care for a simple pole
464 return (sin(x)/cos(x)).series(rel, order+2, options);
467 REGISTER_FUNCTION(tan, eval_func(tan_eval).
468 evalf_func(tan_evalf).
469 derivative_func(tan_deriv).
470 series_func(tan_series));
473 // inverse sine (arc sine)
476 static ex asin_evalf(const ex & x)
480 END_TYPECHECK(asin(x))
482 return asin(ex_to_numeric(x)); // -> numeric asin(numeric)
485 static ex asin_eval(const ex & x)
487 if (x.info(info_flags::numeric)) {
492 if (x.is_equal(_ex1_2()))
493 return numeric(1,6)*Pi;
495 if (x.is_equal(_ex1()))
497 // asin(-1/2) -> -Pi/6
498 if (x.is_equal(_ex_1_2()))
499 return numeric(-1,6)*Pi;
501 if (x.is_equal(_ex_1()))
502 return _num_1_2()*Pi;
503 // asin(float) -> float
504 if (!x.info(info_flags::crational))
505 return asin_evalf(x);
508 return asin(x).hold();
511 static ex asin_deriv(const ex & x, unsigned deriv_param)
513 GINAC_ASSERT(deriv_param==0);
515 // d/dx asin(x) -> 1/sqrt(1-x^2)
516 return power(1-power(x,_ex2()),_ex_1_2());
519 REGISTER_FUNCTION(asin, eval_func(asin_eval).
520 evalf_func(asin_evalf).
521 derivative_func(asin_deriv));
524 // inverse cosine (arc cosine)
527 static ex acos_evalf(const ex & x)
531 END_TYPECHECK(acos(x))
533 return acos(ex_to_numeric(x)); // -> numeric acos(numeric)
536 static ex acos_eval(const ex & x)
538 if (x.info(info_flags::numeric)) {
540 if (x.is_equal(_ex1()))
543 if (x.is_equal(_ex1_2()))
548 // acos(-1/2) -> 2/3*Pi
549 if (x.is_equal(_ex_1_2()))
550 return numeric(2,3)*Pi;
552 if (x.is_equal(_ex_1()))
554 // acos(float) -> float
555 if (!x.info(info_flags::crational))
556 return acos_evalf(x);
559 return acos(x).hold();
562 static ex acos_deriv(const ex & x, unsigned deriv_param)
564 GINAC_ASSERT(deriv_param==0);
566 // d/dx acos(x) -> -1/sqrt(1-x^2)
567 return _ex_1()*power(1-power(x,_ex2()),_ex_1_2());
570 REGISTER_FUNCTION(acos, eval_func(acos_eval).
571 evalf_func(acos_evalf).
572 derivative_func(acos_deriv));
575 // inverse tangent (arc tangent)
578 static ex atan_evalf(const ex & x)
582 END_TYPECHECK(atan(x))
584 return atan(ex_to_numeric(x)); // -> numeric atan(numeric)
587 static ex atan_eval(const ex & x)
589 if (x.info(info_flags::numeric)) {
591 if (x.is_equal(_ex0()))
594 if (x.is_equal(_ex1()))
597 if (x.is_equal(_ex_1()))
599 if (x.is_equal(I) || x.is_equal(-I))
600 throw (pole_error("atan_eval(): logarithmic pole",0));
601 // atan(float) -> float
602 if (!x.info(info_flags::crational))
603 return atan_evalf(x);
606 return atan(x).hold();
609 static ex atan_deriv(const ex & x, unsigned deriv_param)
611 GINAC_ASSERT(deriv_param==0);
613 // d/dx atan(x) -> 1/(1+x^2)
614 return power(_ex1()+power(x,_ex2()), _ex_1());
617 static ex atan_series(const ex &arg,
618 const relational &rel,
622 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
624 // Taylor series where there is no pole or cut falls back to atan_deriv.
625 // There are two branch cuts, one runnig from I up the imaginary axis and
626 // one running from -I down the imaginary axis. The points I and -I are
628 // On the branch cuts and the poles series expand
629 // (log(1+I*x)-log(1-I*x))/(2*I)
631 const ex arg_pt = arg.subs(rel);
632 if (!(I*arg_pt).info(info_flags::real))
633 throw do_taylor(); // Re(x) != 0
634 if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1())
635 throw do_taylor(); // Re(x) == 0, but abs(x)<1
636 // care for the poles, using the defining formula for atan()...
637 if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
638 return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
639 if (!(options & series_options::suppress_branchcut)) {
641 // This is the branch cut: assemble the primitive series manually and
642 // then add the corresponding complex step function.
643 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
644 const ex point = rel.rhs();
646 ex replarg = series(atan(arg), *s==foo, order).subs(foo==point);
647 ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2();
648 if ((I*arg_pt)<_ex0())
649 Order0correction += log((I*arg_pt+_ex_1())/(I*arg_pt+_ex1()))*I*_ex_1_2();
651 Order0correction += log((I*arg_pt+_ex1())/(I*arg_pt+_ex_1()))*I*_ex1_2();
653 seq.push_back(expair(Order0correction, _ex0()));
654 seq.push_back(expair(Order(_ex1()), order));
655 return series(replarg - pseries(rel, seq), rel, order);
660 REGISTER_FUNCTION(atan, eval_func(atan_eval).
661 evalf_func(atan_evalf).
662 derivative_func(atan_deriv).
663 series_func(atan_series));
666 // inverse tangent (atan2(y,x))
669 static ex atan2_evalf(const ex & y, const ex & x)
674 END_TYPECHECK(atan2(y,x))
676 return atan(ex_to_numeric(y),ex_to_numeric(x)); // -> numeric atan(numeric)
679 static ex atan2_eval(const ex & y, const ex & x)
681 if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
682 x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
683 return atan2_evalf(y,x);
686 return atan2(y,x).hold();
689 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
691 GINAC_ASSERT(deriv_param<2);
693 if (deriv_param==0) {
695 return x*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
698 return -y*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
701 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
702 evalf_func(atan2_evalf).
703 derivative_func(atan2_deriv));
706 // hyperbolic sine (trigonometric function)
709 static ex sinh_evalf(const ex & x)
713 END_TYPECHECK(sinh(x))
715 return sinh(ex_to_numeric(x)); // -> numeric sinh(numeric)
718 static ex sinh_eval(const ex & x)
720 if (x.info(info_flags::numeric)) {
721 if (x.is_zero()) // sinh(0) -> 0
723 if (!x.info(info_flags::crational)) // sinh(float) -> float
724 return sinh_evalf(x);
727 if ((x/Pi).info(info_flags::numeric) &&
728 ex_to_numeric(x/Pi).real().is_zero()) // sinh(I*x) -> I*sin(x)
731 if (is_ex_exactly_of_type(x, function)) {
733 // sinh(asinh(x)) -> x
734 if (is_ex_the_function(x, asinh))
736 // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2)
737 if (is_ex_the_function(x, acosh))
738 return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2());
739 // sinh(atanh(x)) -> x*(1-x^2)^(-1/2)
740 if (is_ex_the_function(x, atanh))
741 return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
744 return sinh(x).hold();
747 static ex sinh_deriv(const ex & x, unsigned deriv_param)
749 GINAC_ASSERT(deriv_param==0);
751 // d/dx sinh(x) -> cosh(x)
755 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
756 evalf_func(sinh_evalf).
757 derivative_func(sinh_deriv));
760 // hyperbolic cosine (trigonometric function)
763 static ex cosh_evalf(const ex & x)
767 END_TYPECHECK(cosh(x))
769 return cosh(ex_to_numeric(x)); // -> numeric cosh(numeric)
772 static ex cosh_eval(const ex & x)
774 if (x.info(info_flags::numeric)) {
775 if (x.is_zero()) // cosh(0) -> 1
777 if (!x.info(info_flags::crational)) // cosh(float) -> float
778 return cosh_evalf(x);
781 if ((x/Pi).info(info_flags::numeric) &&
782 ex_to_numeric(x/Pi).real().is_zero()) // cosh(I*x) -> cos(x)
785 if (is_ex_exactly_of_type(x, function)) {
787 // cosh(acosh(x)) -> x
788 if (is_ex_the_function(x, acosh))
790 // cosh(asinh(x)) -> (1+x^2)^(1/2)
791 if (is_ex_the_function(x, asinh))
792 return power(_ex1()+power(t,_ex2()),_ex1_2());
793 // cosh(atanh(x)) -> (1-x^2)^(-1/2)
794 if (is_ex_the_function(x, atanh))
795 return power(_ex1()-power(t,_ex2()),_ex_1_2());
798 return cosh(x).hold();
801 static ex cosh_deriv(const ex & x, unsigned deriv_param)
803 GINAC_ASSERT(deriv_param==0);
805 // d/dx cosh(x) -> sinh(x)
809 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
810 evalf_func(cosh_evalf).
811 derivative_func(cosh_deriv));
815 // hyperbolic tangent (trigonometric function)
818 static ex tanh_evalf(const ex & x)
822 END_TYPECHECK(tanh(x))
824 return tanh(ex_to_numeric(x)); // -> numeric tanh(numeric)
827 static ex tanh_eval(const ex & x)
829 if (x.info(info_flags::numeric)) {
830 if (x.is_zero()) // tanh(0) -> 0
832 if (!x.info(info_flags::crational)) // tanh(float) -> float
833 return tanh_evalf(x);
836 if ((x/Pi).info(info_flags::numeric) &&
837 ex_to_numeric(x/Pi).real().is_zero()) // tanh(I*x) -> I*tan(x);
840 if (is_ex_exactly_of_type(x, function)) {
842 // tanh(atanh(x)) -> x
843 if (is_ex_the_function(x, atanh))
845 // tanh(asinh(x)) -> x*(1+x^2)^(-1/2)
846 if (is_ex_the_function(x, asinh))
847 return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
848 // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x
849 if (is_ex_the_function(x, acosh))
850 return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2())*power(t,_ex_1());
853 return tanh(x).hold();
856 static ex tanh_deriv(const ex & x, unsigned deriv_param)
858 GINAC_ASSERT(deriv_param==0);
860 // d/dx tanh(x) -> 1-tanh(x)^2
861 return _ex1()-power(tanh(x),_ex2());
864 static ex tanh_series(const ex &x,
865 const relational &rel,
869 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
871 // Taylor series where there is no pole falls back to tanh_deriv.
872 // On a pole simply expand sinh(x)/cosh(x).
873 const ex x_pt = x.subs(rel);
874 if (!(2*I*x_pt/Pi).info(info_flags::odd))
875 throw do_taylor(); // caught by function::series()
876 // if we got here we have to care for a simple pole
877 return (sinh(x)/cosh(x)).series(rel, order+2, options);
880 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
881 evalf_func(tanh_evalf).
882 derivative_func(tanh_deriv).
883 series_func(tanh_series));
886 // inverse hyperbolic sine (trigonometric function)
889 static ex asinh_evalf(const ex & x)
893 END_TYPECHECK(asinh(x))
895 return asinh(ex_to_numeric(x)); // -> numeric asinh(numeric)
898 static ex asinh_eval(const ex & x)
900 if (x.info(info_flags::numeric)) {
904 // asinh(float) -> float
905 if (!x.info(info_flags::crational))
906 return asinh_evalf(x);
909 return asinh(x).hold();
912 static ex asinh_deriv(const ex & x, unsigned deriv_param)
914 GINAC_ASSERT(deriv_param==0);
916 // d/dx asinh(x) -> 1/sqrt(1+x^2)
917 return power(_ex1()+power(x,_ex2()),_ex_1_2());
920 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
921 evalf_func(asinh_evalf).
922 derivative_func(asinh_deriv));
925 // inverse hyperbolic cosine (trigonometric function)
928 static ex acosh_evalf(const ex & x)
932 END_TYPECHECK(acosh(x))
934 return acosh(ex_to_numeric(x)); // -> numeric acosh(numeric)
937 static ex acosh_eval(const ex & x)
939 if (x.info(info_flags::numeric)) {
940 // acosh(0) -> Pi*I/2
942 return Pi*I*numeric(1,2);
944 if (x.is_equal(_ex1()))
947 if (x.is_equal(_ex_1()))
949 // acosh(float) -> float
950 if (!x.info(info_flags::crational))
951 return acosh_evalf(x);
954 return acosh(x).hold();
957 static ex acosh_deriv(const ex & x, unsigned deriv_param)
959 GINAC_ASSERT(deriv_param==0);
961 // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
962 return power(x+_ex_1(),_ex_1_2())*power(x+_ex1(),_ex_1_2());
965 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
966 evalf_func(acosh_evalf).
967 derivative_func(acosh_deriv));
970 // inverse hyperbolic tangent (trigonometric function)
973 static ex atanh_evalf(const ex & x)
977 END_TYPECHECK(atanh(x))
979 return atanh(ex_to_numeric(x)); // -> numeric atanh(numeric)
982 static ex atanh_eval(const ex & x)
984 if (x.info(info_flags::numeric)) {
988 // atanh({+|-}1) -> throw
989 if (x.is_equal(_ex1()) || x.is_equal(_ex_1()))
990 throw (pole_error("atanh_eval(): logarithmic pole",0));
991 // atanh(float) -> float
992 if (!x.info(info_flags::crational))
993 return atanh_evalf(x);
996 return atanh(x).hold();
999 static ex atanh_deriv(const ex & x, unsigned deriv_param)
1001 GINAC_ASSERT(deriv_param==0);
1003 // d/dx atanh(x) -> 1/(1-x^2)
1004 return power(_ex1()-power(x,_ex2()),_ex_1());
1007 static ex atanh_series(const ex &arg,
1008 const relational &rel,
1012 GINAC_ASSERT(is_ex_exactly_of_type(rel.lhs(),symbol));
1014 // Taylor series where there is no pole or cut falls back to atanh_deriv.
1015 // There are two branch cuts, one runnig from 1 up the real axis and one
1016 // one running from -1 down the real axis. The points 1 and -1 are poles
1017 // On the branch cuts and the poles series expand
1018 // (log(1+x)-log(1-x))/2
1020 const ex arg_pt = arg.subs(rel);
1021 if (!(arg_pt).info(info_flags::real))
1022 throw do_taylor(); // Im(x) != 0
1023 if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1())
1024 throw do_taylor(); // Im(x) == 0, but abs(x)<1
1025 // care for the poles, using the defining formula for atanh()...
1026 if (arg_pt.is_equal(_ex1()) || arg_pt.is_equal(_ex_1()))
1027 return ((log(_ex1()+arg)-log(_ex1()-arg))*_ex1_2()).series(rel, order, options);
1028 // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
1029 if (!(options & series_options::suppress_branchcut)) {
1031 // This is the branch cut: assemble the primitive series manually and
1032 // then add the corresponding complex step function.
1033 const symbol *s = static_cast<symbol *>(rel.lhs().bp);
1034 const ex point = rel.rhs();
1036 ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point);
1037 ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2();
1039 Order0correction += log((arg_pt+_ex_1())/(arg_pt+_ex1()))*_ex1_2();
1041 Order0correction += log((arg_pt+_ex1())/(arg_pt+_ex_1()))*_ex_1_2();
1043 seq.push_back(expair(Order0correction, _ex0()));
1044 seq.push_back(expair(Order(_ex1()), order));
1045 return series(replarg - pseries(rel, seq), rel, order);
1050 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
1051 evalf_func(atanh_evalf).
1052 derivative_func(atanh_deriv).
1053 series_func(atanh_series));
1056 } // namespace GiNaC