[PATCH 1/3] Automatic evaluation of (e^t)^s = e^(ts).
[ginac.git] / ginac / inifcns_trans.cpp
1 /** @file inifcns_trans.cpp
2  *
3  *  Implementation of transcendental (and trigonometric and hyperbolic)
4  *  functions. */
5
6 /*
7  *  GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
8  *
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.
13  *
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.
18  *
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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22  */
23
24 #include "inifcns.h"
25 #include "ex.h"
26 #include "constant.h"
27 #include "add.h"
28 #include "mul.h"
29 #include "numeric.h"
30 #include "power.h"
31 #include "operators.h"
32 #include "relational.h"
33 #include "symbol.h"
34 #include "pseries.h"
35 #include "utils.h"
36
37 #include <stdexcept>
38 #include <vector>
39
40 namespace GiNaC {
41
42 //////////
43 // exponential function
44 //////////
45
46 static ex exp_evalf(const ex & x)
47 {
48         if (is_exactly_a<numeric>(x))
49                 return exp(ex_to<numeric>(x));
50         
51         return exp(x).hold();
52 }
53
54 static ex exp_eval(const ex & x)
55 {
56         // exp(0) -> 1
57         if (x.is_zero()) {
58                 return _ex1;
59         }
60
61         // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
62         const ex TwoExOverPiI=(_ex2*x)/(Pi*I);
63         if (TwoExOverPiI.info(info_flags::integer)) {
64                 const numeric z = mod(ex_to<numeric>(TwoExOverPiI),*_num4_p);
65                 if (z.is_equal(*_num0_p))
66                         return _ex1;
67                 if (z.is_equal(*_num1_p))
68                         return ex(I);
69                 if (z.is_equal(*_num2_p))
70                         return _ex_1;
71                 if (z.is_equal(*_num3_p))
72                         return ex(-I);
73         }
74
75         // exp(log(x)) -> x
76         if (is_ex_the_function(x, log))
77                 return x.op(0);
78         
79         // exp(float) -> float
80         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
81                 return exp(ex_to<numeric>(x));
82         
83         return exp(x).hold();
84 }
85
86 static ex exp_expand(const ex & arg, unsigned options)
87 {
88         ex exp_arg;
89         if (options & expand_options::expand_function_args)
90                 exp_arg = arg.expand(options);
91         else
92                 exp_arg=arg;
93
94         if ((options & expand_options::expand_transcendental)
95                 && is_exactly_a<add>(exp_arg)) {
96                 exvector prodseq;
97                 prodseq.reserve(exp_arg.nops());
98                 for (const_iterator i = exp_arg.begin(); i != exp_arg.end(); ++i)
99                         prodseq.push_back(exp(*i));
100
101                 return dynallocate<mul>(prodseq).setflag(status_flags::expanded);
102         }
103
104         return exp(exp_arg).hold();
105 }
106
107 static ex exp_deriv(const ex & x, unsigned deriv_param)
108 {
109         GINAC_ASSERT(deriv_param==0);
110
111         // d/dx exp(x) -> exp(x)
112         return exp(x);
113 }
114
115 static ex exp_real_part(const ex & x)
116 {
117         return exp(GiNaC::real_part(x))*cos(GiNaC::imag_part(x));
118 }
119
120 static ex exp_imag_part(const ex & x)
121 {
122         return exp(GiNaC::real_part(x))*sin(GiNaC::imag_part(x));
123 }
124
125 static ex exp_conjugate(const ex & x)
126 {
127         // conjugate(exp(x))==exp(conjugate(x))
128         return exp(x.conjugate());
129 }
130
131 static ex exp_power(const ex & x, const ex & a)
132 {
133         /*
134          * The power law (e^x)^a=e^(x*a) is used in two cases:
135          * a) a is an integer and x may be complex;
136          * b) both x and a are reals.
137          * Negative a is excluded to keep automatic simplifications like exp(x)/exp(x)=1.
138          */
139         if (a.info(info_flags::nonnegative)
140             && (a.info(info_flags::integer) || (x.info(info_flags::real) && a.info(info_flags::real))))
141                 return exp(x*a);
142         else if (a.info(info_flags::negative)
143                  && (a.info(info_flags::integer) || (x.info(info_flags::real) && a.info(info_flags::real))))
144                 return power(exp(-x*a), _ex_1).hold();
145
146         return power(exp(x), a).hold();
147 }
148
149 REGISTER_FUNCTION(exp, eval_func(exp_eval).
150                        evalf_func(exp_evalf).
151                        expand_func(exp_expand).
152                        derivative_func(exp_deriv).
153                        real_part_func(exp_real_part).
154                        imag_part_func(exp_imag_part).
155                        conjugate_func(exp_conjugate).
156                        power_func(exp_power).
157                        latex_name("\\exp"));
158
159 //////////
160 // natural logarithm
161 //////////
162
163 static ex log_evalf(const ex & x)
164 {
165         if (is_exactly_a<numeric>(x))
166                 return log(ex_to<numeric>(x));
167         
168         return log(x).hold();
169 }
170
171 static ex log_eval(const ex & x)
172 {
173         if (x.info(info_flags::numeric)) {
174                 if (x.is_zero())         // log(0) -> infinity
175                         throw(pole_error("log_eval(): log(0)",0));
176                 if (x.info(info_flags::rational) && x.info(info_flags::negative))
177                         return (log(-x)+I*Pi);
178                 if (x.is_equal(_ex1))  // log(1) -> 0
179                         return _ex0;
180                 if (x.is_equal(I))       // log(I) -> Pi*I/2
181                         return (Pi*I*_ex1_2);
182                 if (x.is_equal(-I))      // log(-I) -> -Pi*I/2
183                         return (Pi*I*_ex_1_2);
184
185                 // log(float) -> float
186                 if (!x.info(info_flags::crational))
187                         return log(ex_to<numeric>(x));
188         }
189
190         // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
191         if (is_ex_the_function(x, exp)) {
192                 const ex &t = x.op(0);
193                 if (t.info(info_flags::real))
194                         return t;
195         }
196
197         return log(x).hold();
198 }
199
200 static ex log_deriv(const ex & x, unsigned deriv_param)
201 {
202         GINAC_ASSERT(deriv_param==0);
203         
204         // d/dx log(x) -> 1/x
205         return power(x, _ex_1);
206 }
207
208 static ex log_series(const ex &arg,
209                      const relational &rel,
210                      int order,
211                      unsigned options)
212 {
213         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
214         ex arg_pt;
215         bool must_expand_arg = false;
216         // maybe substitution of rel into arg fails because of a pole
217         try {
218                 arg_pt = arg.subs(rel, subs_options::no_pattern);
219         } catch (pole_error &) {
220                 must_expand_arg = true;
221         }
222         // or we are at the branch point anyways
223         if (arg_pt.is_zero())
224                 must_expand_arg = true;
225         
226         if (arg.diff(ex_to<symbol>(rel.lhs())).is_zero()) {
227                 throw do_taylor();
228         }
229
230         if (must_expand_arg) {
231                 // method:
232                 // This is the branch point: Series expand the argument first, then
233                 // trivially factorize it to isolate that part which has constant
234                 // leading coefficient in this fashion:
235                 //   x^n + x^(n+1) +...+ Order(x^(n+m))  ->  x^n * (1 + x +...+ Order(x^m)).
236                 // Return a plain n*log(x) for the x^n part and series expand the
237                 // other part.  Add them together and reexpand again in order to have
238                 // one unnested pseries object.  All this also works for negative n.
239                 pseries argser;          // series expansion of log's argument
240                 unsigned extra_ord = 0;  // extra expansion order
241                 do {
242                         // oops, the argument expanded to a pure Order(x^something)...
243                         argser = ex_to<pseries>(arg.series(rel, order+extra_ord, options));
244                         ++extra_ord;
245                 } while (!argser.is_terminating() && argser.nops()==1);
246
247                 const symbol &s = ex_to<symbol>(rel.lhs());
248                 const ex &point = rel.rhs();
249                 const int n = argser.ldegree(s);
250                 epvector seq;
251                 // construct what we carelessly called the n*log(x) term above
252                 const ex coeff = argser.coeff(s, n);
253                 // expand the log, but only if coeff is real and > 0, since otherwise
254                 // it would make the branch cut run into the wrong direction
255                 if (coeff.info(info_flags::positive))
256                         seq.push_back(expair(n*log(s-point)+log(coeff), _ex0));
257                 else
258                         seq.push_back(expair(log(coeff*pow(s-point, n)), _ex0));
259
260                 if (!argser.is_terminating() || argser.nops()!=1) {
261                         // in this case n more (or less) terms are needed
262                         // (sadly, to generate them, we have to start from the beginning)
263                         if (n == 0 && coeff == 1) {
264                                 ex rest = pseries(rel, epvector{expair(-1, _ex0), expair(Order(_ex1), order)}).add_series(argser);
265                                 ex acc = dynallocate<pseries>(rel, epvector());
266                                 for (int i = order-1; i>0; --i) {
267                                         epvector cterm { expair(i%2 ? _ex1/i : _ex_1/i, _ex0) };
268                                         acc = pseries(rel, std::move(cterm)).add_series(ex_to<pseries>(acc));
269                                         acc = (ex_to<pseries>(rest)).mul_series(ex_to<pseries>(acc));
270                                 }
271                                 return acc;
272                         }
273                         const ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
274                         return pseries(rel, std::move(seq)).add_series(ex_to<pseries>(log(newarg).series(rel, order, options)));
275                 } else  // it was a monomial
276                         return pseries(rel, std::move(seq));
277         }
278         if (!(options & series_options::suppress_branchcut) &&
279              arg_pt.info(info_flags::negative)) {
280                 // method:
281                 // This is the branch cut: assemble the primitive series manually and
282                 // then add the corresponding complex step function.
283                 const symbol &s = ex_to<symbol>(rel.lhs());
284                 const ex &point = rel.rhs();
285                 const symbol foo;
286                 const ex replarg = series(log(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
287                 epvector seq;
288                 if (order > 0) {
289                         seq.reserve(2);
290                         seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0));
291                 }
292                 seq.push_back(expair(Order(_ex1), order));
293                 return series(replarg - I*Pi + pseries(rel, std::move(seq)), rel, order);
294         }
295         throw do_taylor();  // caught by function::series()
296 }
297
298 static ex log_real_part(const ex & x)
299 {
300         if (x.info(info_flags::nonnegative))
301                 return log(x).hold();
302         return log(abs(x));
303 }
304
305 static ex log_imag_part(const ex & x)
306 {
307         if (x.info(info_flags::nonnegative))
308                 return 0;
309         return atan2(GiNaC::imag_part(x), GiNaC::real_part(x));
310 }
311
312 static ex log_expand(const ex & arg, unsigned options)
313 {
314         if ((options & expand_options::expand_transcendental)
315                 && is_exactly_a<mul>(arg) && !arg.info(info_flags::indefinite)) {
316                 exvector sumseq;
317                 exvector prodseq;
318                 sumseq.reserve(arg.nops());
319                 prodseq.reserve(arg.nops());
320                 bool possign=true;
321
322                 // searching for positive/negative factors
323                 for (const_iterator i = arg.begin(); i != arg.end(); ++i) {
324                         ex e;
325                         if (options & expand_options::expand_function_args)
326                                 e=i->expand(options);
327                         else
328                                 e=*i;
329                         if (e.info(info_flags::positive))
330                                 sumseq.push_back(log(e));
331                         else if (e.info(info_flags::negative)) {
332                                 sumseq.push_back(log(-e));
333                                 possign = !possign;
334                         } else
335                                 prodseq.push_back(e);
336                 }
337
338                 if (sumseq.size() > 0) {
339                         ex newarg;
340                         if (options & expand_options::expand_function_args)
341                                 newarg=((possign?_ex1:_ex_1)*mul(prodseq)).expand(options);
342                         else {
343                                 newarg=(possign?_ex1:_ex_1)*mul(prodseq);
344                                 ex_to<basic>(newarg).setflag(status_flags::purely_indefinite);
345                         }
346                         return add(sumseq)+log(newarg);
347                 } else {
348                         if (!(options & expand_options::expand_function_args))
349                                 ex_to<basic>(arg).setflag(status_flags::purely_indefinite);
350                 }
351         }
352
353         if (options & expand_options::expand_function_args)
354                 return log(arg.expand(options)).hold();
355         else
356                 return log(arg).hold();
357 }
358
359 static ex log_conjugate(const ex & x)
360 {
361         // conjugate(log(x))==log(conjugate(x)) unless on the branch cut which
362         // runs along the negative real axis.
363         if (x.info(info_flags::positive)) {
364                 return log(x);
365         }
366         if (is_exactly_a<numeric>(x) &&
367             !x.imag_part().is_zero()) {
368                 return log(x.conjugate());
369         }
370         return conjugate_function(log(x)).hold();
371 }
372
373 REGISTER_FUNCTION(log, eval_func(log_eval).
374                        evalf_func(log_evalf).
375                        expand_func(log_expand).
376                        derivative_func(log_deriv).
377                        series_func(log_series).
378                        real_part_func(log_real_part).
379                        imag_part_func(log_imag_part).
380                        conjugate_func(log_conjugate).
381                        latex_name("\\ln"));
382
383 //////////
384 // sine (trigonometric function)
385 //////////
386
387 static ex sin_evalf(const ex & x)
388 {
389         if (is_exactly_a<numeric>(x))
390                 return sin(ex_to<numeric>(x));
391         
392         return sin(x).hold();
393 }
394
395 static ex sin_eval(const ex & x)
396 {
397         // sin(n/d*Pi) -> { all known non-nested radicals }
398         const ex SixtyExOverPi = _ex60*x/Pi;
399         ex sign = _ex1;
400         if (SixtyExOverPi.info(info_flags::integer)) {
401                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p);
402                 if (z>=*_num60_p) {
403                         // wrap to interval [0, Pi)
404                         z -= *_num60_p;
405                         sign = _ex_1;
406                 }
407                 if (z>*_num30_p) {
408                         // wrap to interval [0, Pi/2)
409                         z = *_num60_p-z;
410                 }
411                 if (z.is_equal(*_num0_p))  // sin(0)       -> 0
412                         return _ex0;
413                 if (z.is_equal(*_num5_p))  // sin(Pi/12)   -> sqrt(6)/4*(1-sqrt(3)/3)
414                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
415                 if (z.is_equal(*_num6_p))  // sin(Pi/10)   -> sqrt(5)/4-1/4
416                         return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
417                 if (z.is_equal(*_num10_p)) // sin(Pi/6)    -> 1/2
418                         return sign*_ex1_2;
419                 if (z.is_equal(*_num15_p)) // sin(Pi/4)    -> sqrt(2)/2
420                         return sign*_ex1_2*sqrt(_ex2);
421                 if (z.is_equal(*_num18_p)) // sin(3/10*Pi) -> sqrt(5)/4+1/4
422                         return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
423                 if (z.is_equal(*_num20_p)) // sin(Pi/3)    -> sqrt(3)/2
424                         return sign*_ex1_2*sqrt(_ex3);
425                 if (z.is_equal(*_num25_p)) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
426                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
427                 if (z.is_equal(*_num30_p)) // sin(Pi/2)    -> 1
428                         return sign;
429         }
430
431         if (is_exactly_a<function>(x)) {
432                 const ex &t = x.op(0);
433
434                 // sin(asin(x)) -> x
435                 if (is_ex_the_function(x, asin))
436                         return t;
437
438                 // sin(acos(x)) -> sqrt(1-x^2)
439                 if (is_ex_the_function(x, acos))
440                         return sqrt(_ex1-power(t,_ex2));
441
442                 // sin(atan(x)) -> x/sqrt(1+x^2)
443                 if (is_ex_the_function(x, atan))
444                         return t*power(_ex1+power(t,_ex2),_ex_1_2);
445         }
446         
447         // sin(float) -> float
448         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
449                 return sin(ex_to<numeric>(x));
450
451         // sin() is odd
452         if (x.info(info_flags::negative))
453                 return -sin(-x);
454         
455         return sin(x).hold();
456 }
457
458 static ex sin_deriv(const ex & x, unsigned deriv_param)
459 {
460         GINAC_ASSERT(deriv_param==0);
461         
462         // d/dx sin(x) -> cos(x)
463         return cos(x);
464 }
465
466 static ex sin_real_part(const ex & x)
467 {
468         return cosh(GiNaC::imag_part(x))*sin(GiNaC::real_part(x));
469 }
470
471 static ex sin_imag_part(const ex & x)
472 {
473         return sinh(GiNaC::imag_part(x))*cos(GiNaC::real_part(x));
474 }
475
476 static ex sin_conjugate(const ex & x)
477 {
478         // conjugate(sin(x))==sin(conjugate(x))
479         return sin(x.conjugate());
480 }
481
482 REGISTER_FUNCTION(sin, eval_func(sin_eval).
483                        evalf_func(sin_evalf).
484                        derivative_func(sin_deriv).
485                        real_part_func(sin_real_part).
486                        imag_part_func(sin_imag_part).
487                        conjugate_func(sin_conjugate).
488                        latex_name("\\sin"));
489
490 //////////
491 // cosine (trigonometric function)
492 //////////
493
494 static ex cos_evalf(const ex & x)
495 {
496         if (is_exactly_a<numeric>(x))
497                 return cos(ex_to<numeric>(x));
498         
499         return cos(x).hold();
500 }
501
502 static ex cos_eval(const ex & x)
503 {
504         // cos(n/d*Pi) -> { all known non-nested radicals }
505         const ex SixtyExOverPi = _ex60*x/Pi;
506         ex sign = _ex1;
507         if (SixtyExOverPi.info(info_flags::integer)) {
508                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p);
509                 if (z>=*_num60_p) {
510                         // wrap to interval [0, Pi)
511                         z = *_num120_p-z;
512                 }
513                 if (z>=*_num30_p) {
514                         // wrap to interval [0, Pi/2)
515                         z = *_num60_p-z;
516                         sign = _ex_1;
517                 }
518                 if (z.is_equal(*_num0_p))  // cos(0)       -> 1
519                         return sign;
520                 if (z.is_equal(*_num5_p))  // cos(Pi/12)   -> sqrt(6)/4*(1+sqrt(3)/3)
521                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
522                 if (z.is_equal(*_num10_p)) // cos(Pi/6)    -> sqrt(3)/2
523                         return sign*_ex1_2*sqrt(_ex3);
524                 if (z.is_equal(*_num12_p)) // cos(Pi/5)    -> sqrt(5)/4+1/4
525                         return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
526                 if (z.is_equal(*_num15_p)) // cos(Pi/4)    -> sqrt(2)/2
527                         return sign*_ex1_2*sqrt(_ex2);
528                 if (z.is_equal(*_num20_p)) // cos(Pi/3)    -> 1/2
529                         return sign*_ex1_2;
530                 if (z.is_equal(*_num24_p)) // cos(2/5*Pi)  -> sqrt(5)/4-1/4x
531                         return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
532                 if (z.is_equal(*_num25_p)) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
533                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
534                 if (z.is_equal(*_num30_p)) // cos(Pi/2)    -> 0
535                         return _ex0;
536         }
537
538         if (is_exactly_a<function>(x)) {
539                 const ex &t = x.op(0);
540
541                 // cos(acos(x)) -> x
542                 if (is_ex_the_function(x, acos))
543                         return t;
544
545                 // cos(asin(x)) -> sqrt(1-x^2)
546                 if (is_ex_the_function(x, asin))
547                         return sqrt(_ex1-power(t,_ex2));
548
549                 // cos(atan(x)) -> 1/sqrt(1+x^2)
550                 if (is_ex_the_function(x, atan))
551                         return power(_ex1+power(t,_ex2),_ex_1_2);
552         }
553         
554         // cos(float) -> float
555         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
556                 return cos(ex_to<numeric>(x));
557         
558         // cos() is even
559         if (x.info(info_flags::negative))
560                 return cos(-x);
561         
562         return cos(x).hold();
563 }
564
565 static ex cos_deriv(const ex & x, unsigned deriv_param)
566 {
567         GINAC_ASSERT(deriv_param==0);
568
569         // d/dx cos(x) -> -sin(x)
570         return -sin(x);
571 }
572
573 static ex cos_real_part(const ex & x)
574 {
575         return cosh(GiNaC::imag_part(x))*cos(GiNaC::real_part(x));
576 }
577
578 static ex cos_imag_part(const ex & x)
579 {
580         return -sinh(GiNaC::imag_part(x))*sin(GiNaC::real_part(x));
581 }
582
583 static ex cos_conjugate(const ex & x)
584 {
585         // conjugate(cos(x))==cos(conjugate(x))
586         return cos(x.conjugate());
587 }
588
589 REGISTER_FUNCTION(cos, eval_func(cos_eval).
590                        evalf_func(cos_evalf).
591                        derivative_func(cos_deriv).
592                        real_part_func(cos_real_part).
593                        imag_part_func(cos_imag_part).
594                        conjugate_func(cos_conjugate).
595                        latex_name("\\cos"));
596
597 //////////
598 // tangent (trigonometric function)
599 //////////
600
601 static ex tan_evalf(const ex & x)
602 {
603         if (is_exactly_a<numeric>(x))
604                 return tan(ex_to<numeric>(x));
605         
606         return tan(x).hold();
607 }
608
609 static ex tan_eval(const ex & x)
610 {
611         // tan(n/d*Pi) -> { all known non-nested radicals }
612         const ex SixtyExOverPi = _ex60*x/Pi;
613         ex sign = _ex1;
614         if (SixtyExOverPi.info(info_flags::integer)) {
615                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num60_p);
616                 if (z>=*_num60_p) {
617                         // wrap to interval [0, Pi)
618                         z -= *_num60_p;
619                 }
620                 if (z>=*_num30_p) {
621                         // wrap to interval [0, Pi/2)
622                         z = *_num60_p-z;
623                         sign = _ex_1;
624                 }
625                 if (z.is_equal(*_num0_p))  // tan(0)       -> 0
626                         return _ex0;
627                 if (z.is_equal(*_num5_p))  // tan(Pi/12)   -> 2-sqrt(3)
628                         return sign*(_ex2-sqrt(_ex3));
629                 if (z.is_equal(*_num10_p)) // tan(Pi/6)    -> sqrt(3)/3
630                         return sign*_ex1_3*sqrt(_ex3);
631                 if (z.is_equal(*_num15_p)) // tan(Pi/4)    -> 1
632                         return sign;
633                 if (z.is_equal(*_num20_p)) // tan(Pi/3)    -> sqrt(3)
634                         return sign*sqrt(_ex3);
635                 if (z.is_equal(*_num25_p)) // tan(5/12*Pi) -> 2+sqrt(3)
636                         return sign*(sqrt(_ex3)+_ex2);
637                 if (z.is_equal(*_num30_p)) // tan(Pi/2)    -> infinity
638                         throw (pole_error("tan_eval(): simple pole",1));
639         }
640
641         if (is_exactly_a<function>(x)) {
642                 const ex &t = x.op(0);
643
644                 // tan(atan(x)) -> x
645                 if (is_ex_the_function(x, atan))
646                         return t;
647
648                 // tan(asin(x)) -> x/sqrt(1+x^2)
649                 if (is_ex_the_function(x, asin))
650                         return t*power(_ex1-power(t,_ex2),_ex_1_2);
651
652                 // tan(acos(x)) -> sqrt(1-x^2)/x
653                 if (is_ex_the_function(x, acos))
654                         return power(t,_ex_1)*sqrt(_ex1-power(t,_ex2));
655         }
656         
657         // tan(float) -> float
658         if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
659                 return tan(ex_to<numeric>(x));
660         }
661         
662         // tan() is odd
663         if (x.info(info_flags::negative))
664                 return -tan(-x);
665         
666         return tan(x).hold();
667 }
668
669 static ex tan_deriv(const ex & x, unsigned deriv_param)
670 {
671         GINAC_ASSERT(deriv_param==0);
672         
673         // d/dx tan(x) -> 1+tan(x)^2;
674         return (_ex1+power(tan(x),_ex2));
675 }
676
677 static ex tan_real_part(const ex & x)
678 {
679         ex a = GiNaC::real_part(x);
680         ex b = GiNaC::imag_part(x);
681         return tan(a)/(1+power(tan(a),2)*power(tan(b),2));
682 }
683
684 static ex tan_imag_part(const ex & x)
685 {
686         ex a = GiNaC::real_part(x);
687         ex b = GiNaC::imag_part(x);
688         return tanh(b)/(1+power(tan(a),2)*power(tan(b),2));
689 }
690
691 static ex tan_series(const ex &x,
692                      const relational &rel,
693                      int order,
694                      unsigned options)
695 {
696         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
697         // method:
698         // Taylor series where there is no pole falls back to tan_deriv.
699         // On a pole simply expand sin(x)/cos(x).
700         const ex x_pt = x.subs(rel, subs_options::no_pattern);
701         if (!(2*x_pt/Pi).info(info_flags::odd))
702                 throw do_taylor();  // caught by function::series()
703         // if we got here we have to care for a simple pole
704         return (sin(x)/cos(x)).series(rel, order, options);
705 }
706
707 static ex tan_conjugate(const ex & x)
708 {
709         // conjugate(tan(x))==tan(conjugate(x))
710         return tan(x.conjugate());
711 }
712
713 REGISTER_FUNCTION(tan, eval_func(tan_eval).
714                        evalf_func(tan_evalf).
715                        derivative_func(tan_deriv).
716                        series_func(tan_series).
717                        real_part_func(tan_real_part).
718                        imag_part_func(tan_imag_part).
719                        conjugate_func(tan_conjugate).
720                        latex_name("\\tan"));
721
722 //////////
723 // inverse sine (arc sine)
724 //////////
725
726 static ex asin_evalf(const ex & x)
727 {
728         if (is_exactly_a<numeric>(x))
729                 return asin(ex_to<numeric>(x));
730         
731         return asin(x).hold();
732 }
733
734 static ex asin_eval(const ex & x)
735 {
736         if (x.info(info_flags::numeric)) {
737
738                 // asin(0) -> 0
739                 if (x.is_zero())
740                         return x;
741
742                 // asin(1/2) -> Pi/6
743                 if (x.is_equal(_ex1_2))
744                         return numeric(1,6)*Pi;
745
746                 // asin(1) -> Pi/2
747                 if (x.is_equal(_ex1))
748                         return _ex1_2*Pi;
749
750                 // asin(-1/2) -> -Pi/6
751                 if (x.is_equal(_ex_1_2))
752                         return numeric(-1,6)*Pi;
753
754                 // asin(-1) -> -Pi/2
755                 if (x.is_equal(_ex_1))
756                         return _ex_1_2*Pi;
757
758                 // asin(float) -> float
759                 if (!x.info(info_flags::crational))
760                         return asin(ex_to<numeric>(x));
761
762                 // asin() is odd
763                 if (x.info(info_flags::negative))
764                         return -asin(-x);
765         }
766         
767         return asin(x).hold();
768 }
769
770 static ex asin_deriv(const ex & x, unsigned deriv_param)
771 {
772         GINAC_ASSERT(deriv_param==0);
773         
774         // d/dx asin(x) -> 1/sqrt(1-x^2)
775         return power(1-power(x,_ex2),_ex_1_2);
776 }
777
778 static ex asin_conjugate(const ex & x)
779 {
780         // conjugate(asin(x))==asin(conjugate(x)) unless on the branch cuts which
781         // run along the real axis outside the interval [-1, +1].
782         if (is_exactly_a<numeric>(x) &&
783             (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) {
784                 return asin(x.conjugate());
785         }
786         return conjugate_function(asin(x)).hold();
787 }
788
789 REGISTER_FUNCTION(asin, eval_func(asin_eval).
790                         evalf_func(asin_evalf).
791                         derivative_func(asin_deriv).
792                         conjugate_func(asin_conjugate).
793                         latex_name("\\arcsin"));
794
795 //////////
796 // inverse cosine (arc cosine)
797 //////////
798
799 static ex acos_evalf(const ex & x)
800 {
801         if (is_exactly_a<numeric>(x))
802                 return acos(ex_to<numeric>(x));
803         
804         return acos(x).hold();
805 }
806
807 static ex acos_eval(const ex & x)
808 {
809         if (x.info(info_flags::numeric)) {
810
811                 // acos(1) -> 0
812                 if (x.is_equal(_ex1))
813                         return _ex0;
814
815                 // acos(1/2) -> Pi/3
816                 if (x.is_equal(_ex1_2))
817                         return _ex1_3*Pi;
818
819                 // acos(0) -> Pi/2
820                 if (x.is_zero())
821                         return _ex1_2*Pi;
822
823                 // acos(-1/2) -> 2/3*Pi
824                 if (x.is_equal(_ex_1_2))
825                         return numeric(2,3)*Pi;
826
827                 // acos(-1) -> Pi
828                 if (x.is_equal(_ex_1))
829                         return Pi;
830
831                 // acos(float) -> float
832                 if (!x.info(info_flags::crational))
833                         return acos(ex_to<numeric>(x));
834
835                 // acos(-x) -> Pi-acos(x)
836                 if (x.info(info_flags::negative))
837                         return Pi-acos(-x);
838         }
839         
840         return acos(x).hold();
841 }
842
843 static ex acos_deriv(const ex & x, unsigned deriv_param)
844 {
845         GINAC_ASSERT(deriv_param==0);
846         
847         // d/dx acos(x) -> -1/sqrt(1-x^2)
848         return -power(1-power(x,_ex2),_ex_1_2);
849 }
850
851 static ex acos_conjugate(const ex & x)
852 {
853         // conjugate(acos(x))==acos(conjugate(x)) unless on the branch cuts which
854         // run along the real axis outside the interval [-1, +1].
855         if (is_exactly_a<numeric>(x) &&
856             (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) {
857                 return acos(x.conjugate());
858         }
859         return conjugate_function(acos(x)).hold();
860 }
861
862 REGISTER_FUNCTION(acos, eval_func(acos_eval).
863                         evalf_func(acos_evalf).
864                         derivative_func(acos_deriv).
865                         conjugate_func(acos_conjugate).
866                         latex_name("\\arccos"));
867
868 //////////
869 // inverse tangent (arc tangent)
870 //////////
871
872 static ex atan_evalf(const ex & x)
873 {
874         if (is_exactly_a<numeric>(x))
875                 return atan(ex_to<numeric>(x));
876         
877         return atan(x).hold();
878 }
879
880 static ex atan_eval(const ex & x)
881 {
882         if (x.info(info_flags::numeric)) {
883
884                 // atan(0) -> 0
885                 if (x.is_zero())
886                         return _ex0;
887
888                 // atan(1) -> Pi/4
889                 if (x.is_equal(_ex1))
890                         return _ex1_4*Pi;
891
892                 // atan(-1) -> -Pi/4
893                 if (x.is_equal(_ex_1))
894                         return _ex_1_4*Pi;
895
896                 if (x.is_equal(I) || x.is_equal(-I))
897                         throw (pole_error("atan_eval(): logarithmic pole",0));
898
899                 // atan(float) -> float
900                 if (!x.info(info_flags::crational))
901                         return atan(ex_to<numeric>(x));
902
903                 // atan() is odd
904                 if (x.info(info_flags::negative))
905                         return -atan(-x);
906         }
907         
908         return atan(x).hold();
909 }
910
911 static ex atan_deriv(const ex & x, unsigned deriv_param)
912 {
913         GINAC_ASSERT(deriv_param==0);
914
915         // d/dx atan(x) -> 1/(1+x^2)
916         return power(_ex1+power(x,_ex2), _ex_1);
917 }
918
919 static ex atan_series(const ex &arg,
920                       const relational &rel,
921                       int order,
922                       unsigned options)
923 {
924         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
925         // method:
926         // Taylor series where there is no pole or cut falls back to atan_deriv.
927         // There are two branch cuts, one runnig from I up the imaginary axis and
928         // one running from -I down the imaginary axis.  The points I and -I are
929         // poles.
930         // On the branch cuts and the poles series expand
931         //     (log(1+I*x)-log(1-I*x))/(2*I)
932         // instead.
933         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
934         if (!(I*arg_pt).info(info_flags::real))
935                 throw do_taylor();     // Re(x) != 0
936         if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1)
937                 throw do_taylor();     // Re(x) == 0, but abs(x)<1
938         // care for the poles, using the defining formula for atan()...
939         if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
940                 return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
941         if (!(options & series_options::suppress_branchcut)) {
942                 // method:
943                 // This is the branch cut: assemble the primitive series manually and
944                 // then add the corresponding complex step function.
945                 const symbol &s = ex_to<symbol>(rel.lhs());
946                 const ex &point = rel.rhs();
947                 const symbol foo;
948                 const ex replarg = series(atan(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
949                 ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2;
950                 if ((I*arg_pt)<_ex0)
951                         Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2;
952                 else
953                         Order0correction += log((I*arg_pt+_ex1)/(I*arg_pt+_ex_1))*I*_ex1_2;
954                 epvector seq;
955                 if (order > 0) {
956                         seq.reserve(2);
957                         seq.push_back(expair(Order0correction, _ex0));
958                 }
959                 seq.push_back(expair(Order(_ex1), order));
960                 return series(replarg - pseries(rel, std::move(seq)), rel, order);
961         }
962         throw do_taylor();
963 }
964
965 static ex atan_conjugate(const ex & x)
966 {
967         // conjugate(atan(x))==atan(conjugate(x)) unless on the branch cuts which
968         // run along the imaginary axis outside the interval [-I, +I].
969         if (x.info(info_flags::real))
970                 return atan(x);
971         if (is_exactly_a<numeric>(x)) {
972                 const numeric x_re = ex_to<numeric>(x.real_part());
973                 const numeric x_im = ex_to<numeric>(x.imag_part());
974                 if (!x_re.is_zero() ||
975                     (x_im > *_num_1_p && x_im < *_num1_p))
976                         return atan(x.conjugate());
977         }
978         return conjugate_function(atan(x)).hold();
979 }
980
981 REGISTER_FUNCTION(atan, eval_func(atan_eval).
982                         evalf_func(atan_evalf).
983                         derivative_func(atan_deriv).
984                         series_func(atan_series).
985                         conjugate_func(atan_conjugate).
986                         latex_name("\\arctan"));
987
988 //////////
989 // inverse tangent (atan2(y,x))
990 //////////
991
992 static ex atan2_evalf(const ex &y, const ex &x)
993 {
994         if (is_exactly_a<numeric>(y) && is_exactly_a<numeric>(x))
995                 return atan(ex_to<numeric>(y), ex_to<numeric>(x));
996         
997         return atan2(y, x).hold();
998 }
999
1000 static ex atan2_eval(const ex & y, const ex & x)
1001 {
1002         if (y.is_zero()) {
1003
1004                 // atan2(0, 0) -> 0
1005                 if (x.is_zero())
1006                         return _ex0;
1007
1008                 // atan2(0, x), x real and positive -> 0
1009                 if (x.info(info_flags::positive))
1010                         return _ex0;
1011
1012                 // atan2(0, x), x real and negative -> Pi
1013                 if (x.info(info_flags::negative))
1014                         return Pi;
1015         }
1016
1017         if (x.is_zero()) {
1018
1019                 // atan2(y, 0), y real and positive -> Pi/2
1020                 if (y.info(info_flags::positive))
1021                         return _ex1_2*Pi;
1022
1023                 // atan2(y, 0), y real and negative -> -Pi/2
1024                 if (y.info(info_flags::negative))
1025                         return _ex_1_2*Pi;
1026         }
1027
1028         if (y.is_equal(x)) {
1029
1030                 // atan2(y, y), y real and positive -> Pi/4
1031                 if (y.info(info_flags::positive))
1032                         return _ex1_4*Pi;
1033
1034                 // atan2(y, y), y real and negative -> -3/4*Pi
1035                 if (y.info(info_flags::negative))
1036                         return numeric(-3, 4)*Pi;
1037         }
1038
1039         if (y.is_equal(-x)) {
1040
1041                 // atan2(y, -y), y real and positive -> 3*Pi/4
1042                 if (y.info(info_flags::positive))
1043                         return numeric(3, 4)*Pi;
1044
1045                 // atan2(y, -y), y real and negative -> -Pi/4
1046                 if (y.info(info_flags::negative))
1047                         return _ex_1_4*Pi;
1048         }
1049
1050         // atan2(float, float) -> float
1051         if (is_a<numeric>(y) && !y.info(info_flags::crational) &&
1052             is_a<numeric>(x) && !x.info(info_flags::crational))
1053                 return atan(ex_to<numeric>(y), ex_to<numeric>(x));
1054
1055         // atan2(real, real) -> atan(y/x) +/- Pi
1056         if (y.info(info_flags::real) && x.info(info_flags::real)) {
1057                 if (x.info(info_flags::positive))
1058                         return atan(y/x);
1059
1060                 if (x.info(info_flags::negative)) {
1061                         if (y.info(info_flags::positive))
1062                                 return atan(y/x)+Pi;
1063                         if (y.info(info_flags::negative))
1064                                 return atan(y/x)-Pi;
1065                 }
1066         }
1067
1068         return atan2(y, x).hold();
1069 }    
1070
1071 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
1072 {
1073         GINAC_ASSERT(deriv_param<2);
1074         
1075         if (deriv_param==0) {
1076                 // d/dy atan2(y,x)
1077                 return x*power(power(x,_ex2)+power(y,_ex2),_ex_1);
1078         }
1079         // d/dx atan2(y,x)
1080         return -y*power(power(x,_ex2)+power(y,_ex2),_ex_1);
1081 }
1082
1083 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
1084                          evalf_func(atan2_evalf).
1085                          derivative_func(atan2_deriv));
1086
1087 //////////
1088 // hyperbolic sine (trigonometric function)
1089 //////////
1090
1091 static ex sinh_evalf(const ex & x)
1092 {
1093         if (is_exactly_a<numeric>(x))
1094                 return sinh(ex_to<numeric>(x));
1095         
1096         return sinh(x).hold();
1097 }
1098
1099 static ex sinh_eval(const ex & x)
1100 {
1101         if (x.info(info_flags::numeric)) {
1102
1103                 // sinh(0) -> 0
1104                 if (x.is_zero())
1105                         return _ex0;        
1106
1107                 // sinh(float) -> float
1108                 if (!x.info(info_flags::crational))
1109                         return sinh(ex_to<numeric>(x));
1110
1111                 // sinh() is odd
1112                 if (x.info(info_flags::negative))
1113                         return -sinh(-x);
1114         }
1115         
1116         if ((x/Pi).info(info_flags::numeric) &&
1117                 ex_to<numeric>(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
1118                 return I*sin(x/I);
1119         
1120         if (is_exactly_a<function>(x)) {
1121                 const ex &t = x.op(0);
1122
1123                 // sinh(asinh(x)) -> x
1124                 if (is_ex_the_function(x, asinh))
1125                         return t;
1126
1127                 // sinh(acosh(x)) -> sqrt(x-1) * sqrt(x+1)
1128                 if (is_ex_the_function(x, acosh))
1129                         return sqrt(t-_ex1)*sqrt(t+_ex1);
1130
1131                 // sinh(atanh(x)) -> x/sqrt(1-x^2)
1132                 if (is_ex_the_function(x, atanh))
1133                         return t*power(_ex1-power(t,_ex2),_ex_1_2);
1134         }
1135         
1136         return sinh(x).hold();
1137 }
1138
1139 static ex sinh_deriv(const ex & x, unsigned deriv_param)
1140 {
1141         GINAC_ASSERT(deriv_param==0);
1142         
1143         // d/dx sinh(x) -> cosh(x)
1144         return cosh(x);
1145 }
1146
1147 static ex sinh_real_part(const ex & x)
1148 {
1149         return sinh(GiNaC::real_part(x))*cos(GiNaC::imag_part(x));
1150 }
1151
1152 static ex sinh_imag_part(const ex & x)
1153 {
1154         return cosh(GiNaC::real_part(x))*sin(GiNaC::imag_part(x));
1155 }
1156
1157 static ex sinh_conjugate(const ex & x)
1158 {
1159         // conjugate(sinh(x))==sinh(conjugate(x))
1160         return sinh(x.conjugate());
1161 }
1162
1163 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
1164                         evalf_func(sinh_evalf).
1165                         derivative_func(sinh_deriv).
1166                         real_part_func(sinh_real_part).
1167                         imag_part_func(sinh_imag_part).
1168                         conjugate_func(sinh_conjugate).
1169                         latex_name("\\sinh"));
1170
1171 //////////
1172 // hyperbolic cosine (trigonometric function)
1173 //////////
1174
1175 static ex cosh_evalf(const ex & x)
1176 {
1177         if (is_exactly_a<numeric>(x))
1178                 return cosh(ex_to<numeric>(x));
1179         
1180         return cosh(x).hold();
1181 }
1182
1183 static ex cosh_eval(const ex & x)
1184 {
1185         if (x.info(info_flags::numeric)) {
1186
1187                 // cosh(0) -> 1
1188                 if (x.is_zero())
1189                         return _ex1;
1190
1191                 // cosh(float) -> float
1192                 if (!x.info(info_flags::crational))
1193                         return cosh(ex_to<numeric>(x));
1194
1195                 // cosh() is even
1196                 if (x.info(info_flags::negative))
1197                         return cosh(-x);
1198         }
1199         
1200         if ((x/Pi).info(info_flags::numeric) &&
1201                 ex_to<numeric>(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
1202                 return cos(x/I);
1203         
1204         if (is_exactly_a<function>(x)) {
1205                 const ex &t = x.op(0);
1206
1207                 // cosh(acosh(x)) -> x
1208                 if (is_ex_the_function(x, acosh))
1209                         return t;
1210
1211                 // cosh(asinh(x)) -> sqrt(1+x^2)
1212                 if (is_ex_the_function(x, asinh))
1213                         return sqrt(_ex1+power(t,_ex2));
1214
1215                 // cosh(atanh(x)) -> 1/sqrt(1-x^2)
1216                 if (is_ex_the_function(x, atanh))
1217                         return power(_ex1-power(t,_ex2),_ex_1_2);
1218         }
1219         
1220         return cosh(x).hold();
1221 }
1222
1223 static ex cosh_deriv(const ex & x, unsigned deriv_param)
1224 {
1225         GINAC_ASSERT(deriv_param==0);
1226         
1227         // d/dx cosh(x) -> sinh(x)
1228         return sinh(x);
1229 }
1230
1231 static ex cosh_real_part(const ex & x)
1232 {
1233         return cosh(GiNaC::real_part(x))*cos(GiNaC::imag_part(x));
1234 }
1235
1236 static ex cosh_imag_part(const ex & x)
1237 {
1238         return sinh(GiNaC::real_part(x))*sin(GiNaC::imag_part(x));
1239 }
1240
1241 static ex cosh_conjugate(const ex & x)
1242 {
1243         // conjugate(cosh(x))==cosh(conjugate(x))
1244         return cosh(x.conjugate());
1245 }
1246
1247 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
1248                         evalf_func(cosh_evalf).
1249                         derivative_func(cosh_deriv).
1250                         real_part_func(cosh_real_part).
1251                         imag_part_func(cosh_imag_part).
1252                         conjugate_func(cosh_conjugate).
1253                         latex_name("\\cosh"));
1254
1255 //////////
1256 // hyperbolic tangent (trigonometric function)
1257 //////////
1258
1259 static ex tanh_evalf(const ex & x)
1260 {
1261         if (is_exactly_a<numeric>(x))
1262                 return tanh(ex_to<numeric>(x));
1263         
1264         return tanh(x).hold();
1265 }
1266
1267 static ex tanh_eval(const ex & x)
1268 {
1269         if (x.info(info_flags::numeric)) {
1270
1271                 // tanh(0) -> 0
1272                 if (x.is_zero())
1273                         return _ex0;
1274
1275                 // tanh(float) -> float
1276                 if (!x.info(info_flags::crational))
1277                         return tanh(ex_to<numeric>(x));
1278
1279                 // tanh() is odd
1280                 if (x.info(info_flags::negative))
1281                         return -tanh(-x);
1282         }
1283         
1284         if ((x/Pi).info(info_flags::numeric) &&
1285                 ex_to<numeric>(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
1286                 return I*tan(x/I);
1287         
1288         if (is_exactly_a<function>(x)) {
1289                 const ex &t = x.op(0);
1290
1291                 // tanh(atanh(x)) -> x
1292                 if (is_ex_the_function(x, atanh))
1293                         return t;
1294
1295                 // tanh(asinh(x)) -> x/sqrt(1+x^2)
1296                 if (is_ex_the_function(x, asinh))
1297                         return t*power(_ex1+power(t,_ex2),_ex_1_2);
1298
1299                 // tanh(acosh(x)) -> sqrt(x-1)*sqrt(x+1)/x
1300                 if (is_ex_the_function(x, acosh))
1301                         return sqrt(t-_ex1)*sqrt(t+_ex1)*power(t,_ex_1);
1302         }
1303         
1304         return tanh(x).hold();
1305 }
1306
1307 static ex tanh_deriv(const ex & x, unsigned deriv_param)
1308 {
1309         GINAC_ASSERT(deriv_param==0);
1310         
1311         // d/dx tanh(x) -> 1-tanh(x)^2
1312         return _ex1-power(tanh(x),_ex2);
1313 }
1314
1315 static ex tanh_series(const ex &x,
1316                       const relational &rel,
1317                       int order,
1318                       unsigned options)
1319 {
1320         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
1321         // method:
1322         // Taylor series where there is no pole falls back to tanh_deriv.
1323         // On a pole simply expand sinh(x)/cosh(x).
1324         const ex x_pt = x.subs(rel, subs_options::no_pattern);
1325         if (!(2*I*x_pt/Pi).info(info_flags::odd))
1326                 throw do_taylor();  // caught by function::series()
1327         // if we got here we have to care for a simple pole
1328         return (sinh(x)/cosh(x)).series(rel, order, options);
1329 }
1330
1331 static ex tanh_real_part(const ex & x)
1332 {
1333         ex a = GiNaC::real_part(x);
1334         ex b = GiNaC::imag_part(x);
1335         return tanh(a)/(1+power(tanh(a),2)*power(tan(b),2));
1336 }
1337
1338 static ex tanh_imag_part(const ex & x)
1339 {
1340         ex a = GiNaC::real_part(x);
1341         ex b = GiNaC::imag_part(x);
1342         return tan(b)/(1+power(tanh(a),2)*power(tan(b),2));
1343 }
1344
1345 static ex tanh_conjugate(const ex & x)
1346 {
1347         // conjugate(tanh(x))==tanh(conjugate(x))
1348         return tanh(x.conjugate());
1349 }
1350
1351 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
1352                         evalf_func(tanh_evalf).
1353                         derivative_func(tanh_deriv).
1354                         series_func(tanh_series).
1355                         real_part_func(tanh_real_part).
1356                         imag_part_func(tanh_imag_part).
1357                         conjugate_func(tanh_conjugate).
1358                         latex_name("\\tanh"));
1359
1360 //////////
1361 // inverse hyperbolic sine (trigonometric function)
1362 //////////
1363
1364 static ex asinh_evalf(const ex & x)
1365 {
1366         if (is_exactly_a<numeric>(x))
1367                 return asinh(ex_to<numeric>(x));
1368         
1369         return asinh(x).hold();
1370 }
1371
1372 static ex asinh_eval(const ex & x)
1373 {
1374         if (x.info(info_flags::numeric)) {
1375
1376                 // asinh(0) -> 0
1377                 if (x.is_zero())
1378                         return _ex0;
1379
1380                 // asinh(float) -> float
1381                 if (!x.info(info_flags::crational))
1382                         return asinh(ex_to<numeric>(x));
1383
1384                 // asinh() is odd
1385                 if (x.info(info_flags::negative))
1386                         return -asinh(-x);
1387         }
1388         
1389         return asinh(x).hold();
1390 }
1391
1392 static ex asinh_deriv(const ex & x, unsigned deriv_param)
1393 {
1394         GINAC_ASSERT(deriv_param==0);
1395         
1396         // d/dx asinh(x) -> 1/sqrt(1+x^2)
1397         return power(_ex1+power(x,_ex2),_ex_1_2);
1398 }
1399
1400 static ex asinh_conjugate(const ex & x)
1401 {
1402         // conjugate(asinh(x))==asinh(conjugate(x)) unless on the branch cuts which
1403         // run along the imaginary axis outside the interval [-I, +I].
1404         if (x.info(info_flags::real))
1405                 return asinh(x);
1406         if (is_exactly_a<numeric>(x)) {
1407                 const numeric x_re = ex_to<numeric>(x.real_part());
1408                 const numeric x_im = ex_to<numeric>(x.imag_part());
1409                 if (!x_re.is_zero() ||
1410                     (x_im > *_num_1_p && x_im < *_num1_p))
1411                         return asinh(x.conjugate());
1412         }
1413         return conjugate_function(asinh(x)).hold();
1414 }
1415
1416 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
1417                          evalf_func(asinh_evalf).
1418                          derivative_func(asinh_deriv).
1419                          conjugate_func(asinh_conjugate));
1420
1421 //////////
1422 // inverse hyperbolic cosine (trigonometric function)
1423 //////////
1424
1425 static ex acosh_evalf(const ex & x)
1426 {
1427         if (is_exactly_a<numeric>(x))
1428                 return acosh(ex_to<numeric>(x));
1429         
1430         return acosh(x).hold();
1431 }
1432
1433 static ex acosh_eval(const ex & x)
1434 {
1435         if (x.info(info_flags::numeric)) {
1436
1437                 // acosh(0) -> Pi*I/2
1438                 if (x.is_zero())
1439                         return Pi*I*numeric(1,2);
1440
1441                 // acosh(1) -> 0
1442                 if (x.is_equal(_ex1))
1443                         return _ex0;
1444
1445                 // acosh(-1) -> Pi*I
1446                 if (x.is_equal(_ex_1))
1447                         return Pi*I;
1448
1449                 // acosh(float) -> float
1450                 if (!x.info(info_flags::crational))
1451                         return acosh(ex_to<numeric>(x));
1452
1453                 // acosh(-x) -> Pi*I-acosh(x)
1454                 if (x.info(info_flags::negative))
1455                         return Pi*I-acosh(-x);
1456         }
1457         
1458         return acosh(x).hold();
1459 }
1460
1461 static ex acosh_deriv(const ex & x, unsigned deriv_param)
1462 {
1463         GINAC_ASSERT(deriv_param==0);
1464         
1465         // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
1466         return power(x+_ex_1,_ex_1_2)*power(x+_ex1,_ex_1_2);
1467 }
1468
1469 static ex acosh_conjugate(const ex & x)
1470 {
1471         // conjugate(acosh(x))==acosh(conjugate(x)) unless on the branch cut
1472         // which runs along the real axis from +1 to -inf.
1473         if (is_exactly_a<numeric>(x) &&
1474             (!x.imag_part().is_zero() || x > *_num1_p)) {
1475                 return acosh(x.conjugate());
1476         }
1477         return conjugate_function(acosh(x)).hold();
1478 }
1479
1480 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
1481                          evalf_func(acosh_evalf).
1482                          derivative_func(acosh_deriv).
1483                          conjugate_func(acosh_conjugate));
1484
1485 //////////
1486 // inverse hyperbolic tangent (trigonometric function)
1487 //////////
1488
1489 static ex atanh_evalf(const ex & x)
1490 {
1491         if (is_exactly_a<numeric>(x))
1492                 return atanh(ex_to<numeric>(x));
1493         
1494         return atanh(x).hold();
1495 }
1496
1497 static ex atanh_eval(const ex & x)
1498 {
1499         if (x.info(info_flags::numeric)) {
1500
1501                 // atanh(0) -> 0
1502                 if (x.is_zero())
1503                         return _ex0;
1504
1505                 // atanh({+|-}1) -> throw
1506                 if (x.is_equal(_ex1) || x.is_equal(_ex_1))
1507                         throw (pole_error("atanh_eval(): logarithmic pole",0));
1508
1509                 // atanh(float) -> float
1510                 if (!x.info(info_flags::crational))
1511                         return atanh(ex_to<numeric>(x));
1512
1513                 // atanh() is odd
1514                 if (x.info(info_flags::negative))
1515                         return -atanh(-x);
1516         }
1517         
1518         return atanh(x).hold();
1519 }
1520
1521 static ex atanh_deriv(const ex & x, unsigned deriv_param)
1522 {
1523         GINAC_ASSERT(deriv_param==0);
1524         
1525         // d/dx atanh(x) -> 1/(1-x^2)
1526         return power(_ex1-power(x,_ex2),_ex_1);
1527 }
1528
1529 static ex atanh_series(const ex &arg,
1530                        const relational &rel,
1531                        int order,
1532                        unsigned options)
1533 {
1534         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
1535         // method:
1536         // Taylor series where there is no pole or cut falls back to atanh_deriv.
1537         // There are two branch cuts, one runnig from 1 up the real axis and one
1538         // one running from -1 down the real axis.  The points 1 and -1 are poles
1539         // On the branch cuts and the poles series expand
1540         //     (log(1+x)-log(1-x))/2
1541         // instead.
1542         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
1543         if (!(arg_pt).info(info_flags::real))
1544                 throw do_taylor();     // Im(x) != 0
1545         if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1)
1546                 throw do_taylor();     // Im(x) == 0, but abs(x)<1
1547         // care for the poles, using the defining formula for atanh()...
1548         if (arg_pt.is_equal(_ex1) || arg_pt.is_equal(_ex_1))
1549                 return ((log(_ex1+arg)-log(_ex1-arg))*_ex1_2).series(rel, order, options);
1550         // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
1551         if (!(options & series_options::suppress_branchcut)) {
1552                 // method:
1553                 // This is the branch cut: assemble the primitive series manually and
1554                 // then add the corresponding complex step function.
1555                 const symbol &s = ex_to<symbol>(rel.lhs());
1556                 const ex &point = rel.rhs();
1557                 const symbol foo;
1558                 const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
1559                 ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2;
1560                 if (arg_pt<_ex0)
1561                         Order0correction += log((arg_pt+_ex_1)/(arg_pt+_ex1))*_ex1_2;
1562                 else
1563                         Order0correction += log((arg_pt+_ex1)/(arg_pt+_ex_1))*_ex_1_2;
1564                 epvector seq;
1565                 if (order > 0) {
1566                         seq.reserve(2);
1567                         seq.push_back(expair(Order0correction, _ex0));
1568                 }
1569                 seq.push_back(expair(Order(_ex1), order));
1570                 return series(replarg - pseries(rel, std::move(seq)), rel, order);
1571         }
1572         throw do_taylor();
1573 }
1574
1575 static ex atanh_conjugate(const ex & x)
1576 {
1577         // conjugate(atanh(x))==atanh(conjugate(x)) unless on the branch cuts which
1578         // run along the real axis outside the interval [-1, +1].
1579         if (is_exactly_a<numeric>(x) &&
1580             (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) {
1581                 return atanh(x.conjugate());
1582         }
1583         return conjugate_function(atanh(x)).hold();
1584 }
1585
1586 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
1587                          evalf_func(atanh_evalf).
1588                          derivative_func(atanh_deriv).
1589                          series_func(atanh_series).
1590                          conjugate_func(atanh_conjugate));
1591
1592
1593 } // namespace GiNaC