merging 1.2 branch into main trunk
[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-2003 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
22  */
23
24 #include <vector>
25 #include <stdexcept>
26
27 #include "inifcns.h"
28 #include "ex.h"
29 #include "constant.h"
30 #include "numeric.h"
31 #include "power.h"
32 #include "operators.h"
33 #include "relational.h"
34 #include "symbol.h"
35 #include "pseries.h"
36 #include "utils.h"
37
38 namespace GiNaC {
39
40 //////////
41 // exponential function
42 //////////
43
44 static ex exp_evalf(const ex & x)
45 {
46         if (is_exactly_a<numeric>(x))
47                 return exp(ex_to<numeric>(x));
48         
49         return exp(x).hold();
50 }
51
52 static ex exp_eval(const ex & x)
53 {
54         // exp(0) -> 1
55         if (x.is_zero()) {
56                 return _ex1;
57         }
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                 const numeric z = mod(ex_to<numeric>(TwoExOverPiI),_num4);
62                 if (z.is_equal(_num0))
63                         return _ex1;
64                 if (z.is_equal(_num1))
65                         return ex(I);
66                 if (z.is_equal(_num2))
67                         return _ex_1;
68                 if (z.is_equal(_num3))
69                         return ex(-I);
70         }
71         // exp(log(x)) -> x
72         if (is_ex_the_function(x, log))
73                 return x.op(0);
74         
75         // exp(float)
76         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
77                 return exp(ex_to<numeric>(x));
78         
79         return exp(x).hold();
80 }
81
82 static ex exp_deriv(const ex & x, unsigned deriv_param)
83 {
84         GINAC_ASSERT(deriv_param==0);
85
86         // d/dx exp(x) -> exp(x)
87         return exp(x);
88 }
89
90 REGISTER_FUNCTION(exp, eval_func(exp_eval).
91                        evalf_func(exp_evalf).
92                        derivative_func(exp_deriv).
93                        latex_name("\\exp"));
94
95 //////////
96 // natural logarithm
97 //////////
98
99 static ex log_evalf(const ex & x)
100 {
101         if (is_exactly_a<numeric>(x))
102                 return log(ex_to<numeric>(x));
103         
104         return log(x).hold();
105 }
106
107 static ex log_eval(const ex & x)
108 {
109         if (x.info(info_flags::numeric)) {
110                 if (x.is_zero())         // 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
115                         return _ex0;
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);
120                 // log(float)
121                 if (!x.info(info_flags::crational))
122                         return log(ex_to<numeric>(x));
123         }
124         // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
125         if (is_ex_the_function(x, exp)) {
126                 const ex &t = x.op(0);
127                 if (t.info(info_flags::numeric)) {
128                         const numeric &nt = ex_to<numeric>(t);
129                         if (nt.is_real())
130                                 return t;
131                 }
132         }
133         
134         return log(x).hold();
135 }
136
137 static ex log_deriv(const ex & x, unsigned deriv_param)
138 {
139         GINAC_ASSERT(deriv_param==0);
140         
141         // d/dx log(x) -> 1/x
142         return power(x, _ex_1);
143 }
144
145 static ex log_series(const ex &arg,
146                      const relational &rel,
147                      int order,
148                      unsigned options)
149 {
150         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
151         ex arg_pt;
152         bool must_expand_arg = false;
153         // maybe substitution of rel into arg fails because of a pole
154         try {
155                 arg_pt = arg.subs(rel);
156         } catch (pole_error) {
157                 must_expand_arg = true;
158         }
159         // or we are at the branch point anyways
160         if (arg_pt.is_zero())
161                 must_expand_arg = true;
162         
163         if (must_expand_arg) {
164                 // method:
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 + x^(n+1) +...+ Order(x^(n+m))  ->  x^n * (1 + x +...+ 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                 pseries argser;          // series expansion of log's argument
173                 unsigned extra_ord = 0;  // extra expansion order
174                 do {
175                         // oops, the argument expanded to a pure Order(x^something)...
176                         argser = ex_to<pseries>(arg.series(rel, order+extra_ord, options));
177                         ++extra_ord;
178                 } while (!argser.is_terminating() && argser.nops()==1);
179
180                 const symbol &s = ex_to<symbol>(rel.lhs());
181                 const ex &point = rel.rhs();
182                 const int n = argser.ldegree(s);
183                 epvector seq;
184                 // construct what we carelessly called the n*log(x) term above
185                 const ex coeff = argser.coeff(s, n);
186                 // expand the log, but only if coeff is real and > 0, since otherwise
187                 // it would make the branch cut run into the wrong direction
188                 if (coeff.info(info_flags::positive))
189                         seq.push_back(expair(n*log(s-point)+log(coeff), _ex0));
190                 else
191                         seq.push_back(expair(log(coeff*pow(s-point, n)), _ex0));
192
193                 if (!argser.is_terminating() || argser.nops()!=1) {
194                         // in this case n more (or less) terms are needed
195                         // (sadly, to generate them, we have to start from the beginning)
196                         const ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
197                         return pseries(rel, seq).add_series(ex_to<pseries>(log(newarg).series(rel, order, options)));
198                 } else  // it was a monomial
199                         return pseries(rel, seq);
200         }
201         if (!(options & series_options::suppress_branchcut) &&
202              arg_pt.info(info_flags::negative)) {
203                 // method:
204                 // This is the branch cut: assemble the primitive series manually and
205                 // then add the corresponding complex step function.
206                 const symbol &s = ex_to<symbol>(rel.lhs());
207                 const ex &point = rel.rhs();
208                 const symbol foo;
209                 const ex replarg = series(log(arg), s==foo, order).subs(foo==point);
210                 epvector seq;
211                 seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0));
212                 seq.push_back(expair(Order(_ex1), order));
213                 return series(replarg - I*Pi + pseries(rel, seq), rel, order);
214         }
215         throw do_taylor();  // caught by function::series()
216 }
217
218 REGISTER_FUNCTION(log, eval_func(log_eval).
219                        evalf_func(log_evalf).
220                        derivative_func(log_deriv).
221                        series_func(log_series).
222                        latex_name("\\ln"));
223
224 //////////
225 // sine (trigonometric function)
226 //////////
227
228 static ex sin_evalf(const ex & x)
229 {
230         if (is_exactly_a<numeric>(x))
231                 return sin(ex_to<numeric>(x));
232         
233         return sin(x).hold();
234 }
235
236 static ex sin_eval(const ex & x)
237 {
238         // sin(n/d*Pi) -> { all known non-nested radicals }
239         const ex SixtyExOverPi = _ex60*x/Pi;
240         ex sign = _ex1;
241         if (SixtyExOverPi.info(info_flags::integer)) {
242                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120);
243                 if (z>=_num60) {
244                         // wrap to interval [0, Pi)
245                         z -= _num60;
246                         sign = _ex_1;
247                 }
248                 if (z>_num30) {
249                         // wrap to interval [0, Pi/2)
250                         z = _num60-z;
251                 }
252                 if (z.is_equal(_num0))  // sin(0)       -> 0
253                         return _ex0;
254                 if (z.is_equal(_num5))  // sin(Pi/12)   -> sqrt(6)/4*(1-sqrt(3)/3)
255                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
256                 if (z.is_equal(_num6))  // sin(Pi/10)   -> sqrt(5)/4-1/4
257                         return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
258                 if (z.is_equal(_num10)) // sin(Pi/6)    -> 1/2
259                         return sign*_ex1_2;
260                 if (z.is_equal(_num15)) // sin(Pi/4)    -> sqrt(2)/2
261                         return sign*_ex1_2*sqrt(_ex2);
262                 if (z.is_equal(_num18)) // sin(3/10*Pi) -> sqrt(5)/4+1/4
263                         return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
264                 if (z.is_equal(_num20)) // sin(Pi/3)    -> sqrt(3)/2
265                         return sign*_ex1_2*sqrt(_ex3);
266                 if (z.is_equal(_num25)) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
267                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
268                 if (z.is_equal(_num30)) // sin(Pi/2)    -> 1
269                         return sign;
270         }
271         
272         if (is_exactly_a<function>(x)) {
273                 const ex &t = x.op(0);
274                 // sin(asin(x)) -> x
275                 if (is_ex_the_function(x, asin))
276                         return t;
277                 // sin(acos(x)) -> sqrt(1-x^2)
278                 if (is_ex_the_function(x, acos))
279                         return sqrt(_ex1-power(t,_ex2));
280                 // sin(atan(x)) -> x/sqrt(1+x^2)
281                 if (is_ex_the_function(x, atan))
282                         return t*power(_ex1+power(t,_ex2),_ex_1_2);
283         }
284         
285         // sin(float) -> float
286         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
287                 return sin(ex_to<numeric>(x));
288         
289         return sin(x).hold();
290 }
291
292 static ex sin_deriv(const ex & x, unsigned deriv_param)
293 {
294         GINAC_ASSERT(deriv_param==0);
295         
296         // d/dx sin(x) -> cos(x)
297         return cos(x);
298 }
299
300 REGISTER_FUNCTION(sin, eval_func(sin_eval).
301                        evalf_func(sin_evalf).
302                        derivative_func(sin_deriv).
303                        latex_name("\\sin"));
304
305 //////////
306 // cosine (trigonometric function)
307 //////////
308
309 static ex cos_evalf(const ex & x)
310 {
311         if (is_exactly_a<numeric>(x))
312                 return cos(ex_to<numeric>(x));
313         
314         return cos(x).hold();
315 }
316
317 static ex cos_eval(const ex & x)
318 {
319         // cos(n/d*Pi) -> { all known non-nested radicals }
320         const ex SixtyExOverPi = _ex60*x/Pi;
321         ex sign = _ex1;
322         if (SixtyExOverPi.info(info_flags::integer)) {
323                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120);
324                 if (z>=_num60) {
325                         // wrap to interval [0, Pi)
326                         z = _num120-z;
327                 }
328                 if (z>=_num30) {
329                         // wrap to interval [0, Pi/2)
330                         z = _num60-z;
331                         sign = _ex_1;
332                 }
333                 if (z.is_equal(_num0))  // cos(0)       -> 1
334                         return sign;
335                 if (z.is_equal(_num5))  // cos(Pi/12)   -> sqrt(6)/4*(1+sqrt(3)/3)
336                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
337                 if (z.is_equal(_num10)) // cos(Pi/6)    -> sqrt(3)/2
338                         return sign*_ex1_2*sqrt(_ex3);
339                 if (z.is_equal(_num12)) // cos(Pi/5)    -> sqrt(5)/4+1/4
340                         return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
341                 if (z.is_equal(_num15)) // cos(Pi/4)    -> sqrt(2)/2
342                         return sign*_ex1_2*sqrt(_ex2);
343                 if (z.is_equal(_num20)) // cos(Pi/3)    -> 1/2
344                         return sign*_ex1_2;
345                 if (z.is_equal(_num24)) // cos(2/5*Pi)  -> sqrt(5)/4-1/4x
346                         return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
347                 if (z.is_equal(_num25)) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
348                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
349                 if (z.is_equal(_num30)) // cos(Pi/2)    -> 0
350                         return _ex0;
351         }
352         
353         if (is_exactly_a<function>(x)) {
354                 const ex &t = x.op(0);
355                 // cos(acos(x)) -> x
356                 if (is_ex_the_function(x, acos))
357                         return t;
358                 // cos(asin(x)) -> sqrt(1-x^2)
359                 if (is_ex_the_function(x, asin))
360                         return sqrt(_ex1-power(t,_ex2));
361                 // cos(atan(x)) -> 1/sqrt(1+x^2)
362                 if (is_ex_the_function(x, atan))
363                         return power(_ex1+power(t,_ex2),_ex_1_2);
364         }
365         
366         // cos(float) -> float
367         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
368                 return cos(ex_to<numeric>(x));
369         
370         return cos(x).hold();
371 }
372
373 static ex cos_deriv(const ex & x, unsigned deriv_param)
374 {
375         GINAC_ASSERT(deriv_param==0);
376
377         // d/dx cos(x) -> -sin(x)
378         return -sin(x);
379 }
380
381 REGISTER_FUNCTION(cos, eval_func(cos_eval).
382                        evalf_func(cos_evalf).
383                        derivative_func(cos_deriv).
384                        latex_name("\\cos"));
385
386 //////////
387 // tangent (trigonometric function)
388 //////////
389
390 static ex tan_evalf(const ex & x)
391 {
392         if (is_exactly_a<numeric>(x))
393                 return tan(ex_to<numeric>(x));
394         
395         return tan(x).hold();
396 }
397
398 static ex tan_eval(const ex & x)
399 {
400         // tan(n/d*Pi) -> { all known non-nested radicals }
401         const ex SixtyExOverPi = _ex60*x/Pi;
402         ex sign = _ex1;
403         if (SixtyExOverPi.info(info_flags::integer)) {
404                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num60);
405                 if (z>=_num60) {
406                         // wrap to interval [0, Pi)
407                         z -= _num60;
408                 }
409                 if (z>=_num30) {
410                         // wrap to interval [0, Pi/2)
411                         z = _num60-z;
412                         sign = _ex_1;
413                 }
414                 if (z.is_equal(_num0))  // tan(0)       -> 0
415                         return _ex0;
416                 if (z.is_equal(_num5))  // tan(Pi/12)   -> 2-sqrt(3)
417                         return sign*(_ex2-sqrt(_ex3));
418                 if (z.is_equal(_num10)) // tan(Pi/6)    -> sqrt(3)/3
419                         return sign*_ex1_3*sqrt(_ex3);
420                 if (z.is_equal(_num15)) // tan(Pi/4)    -> 1
421                         return sign;
422                 if (z.is_equal(_num20)) // tan(Pi/3)    -> sqrt(3)
423                         return sign*sqrt(_ex3);
424                 if (z.is_equal(_num25)) // tan(5/12*Pi) -> 2+sqrt(3)
425                         return sign*(sqrt(_ex3)+_ex2);
426                 if (z.is_equal(_num30)) // tan(Pi/2)    -> infinity
427                         throw (pole_error("tan_eval(): simple pole",1));
428         }
429         
430         if (is_exactly_a<function>(x)) {
431                 const ex &t = x.op(0);
432                 // tan(atan(x)) -> x
433                 if (is_ex_the_function(x, atan))
434                         return t;
435                 // tan(asin(x)) -> x/sqrt(1+x^2)
436                 if (is_ex_the_function(x, asin))
437                         return t*power(_ex1-power(t,_ex2),_ex_1_2);
438                 // tan(acos(x)) -> sqrt(1-x^2)/x
439                 if (is_ex_the_function(x, acos))
440                         return power(t,_ex_1)*sqrt(_ex1-power(t,_ex2));
441         }
442         
443         // tan(float) -> float
444         if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
445                 return tan(ex_to<numeric>(x));
446         }
447         
448         return tan(x).hold();
449 }
450
451 static ex tan_deriv(const ex & x, unsigned deriv_param)
452 {
453         GINAC_ASSERT(deriv_param==0);
454         
455         // d/dx tan(x) -> 1+tan(x)^2;
456         return (_ex1+power(tan(x),_ex2));
457 }
458
459 static ex tan_series(const ex &x,
460                      const relational &rel,
461                      int order,
462                      unsigned options)
463 {
464         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
465         // method:
466         // Taylor series where there is no pole falls back to tan_deriv.
467         // On a pole simply expand sin(x)/cos(x).
468         const ex x_pt = x.subs(rel);
469         if (!(2*x_pt/Pi).info(info_flags::odd))
470                 throw do_taylor();  // caught by function::series()
471         // if we got here we have to care for a simple pole
472         return (sin(x)/cos(x)).series(rel, order+2, options);
473 }
474
475 REGISTER_FUNCTION(tan, eval_func(tan_eval).
476                        evalf_func(tan_evalf).
477                        derivative_func(tan_deriv).
478                        series_func(tan_series).
479                        latex_name("\\tan"));
480
481 //////////
482 // inverse sine (arc sine)
483 //////////
484
485 static ex asin_evalf(const ex & x)
486 {
487         if (is_exactly_a<numeric>(x))
488                 return asin(ex_to<numeric>(x));
489         
490         return asin(x).hold();
491 }
492
493 static ex asin_eval(const ex & x)
494 {
495         if (x.info(info_flags::numeric)) {
496                 // asin(0) -> 0
497                 if (x.is_zero())
498                         return x;
499                 // asin(1/2) -> Pi/6
500                 if (x.is_equal(_ex1_2))
501                         return numeric(1,6)*Pi;
502                 // asin(1) -> Pi/2
503                 if (x.is_equal(_ex1))
504                         return _num1_2*Pi;
505                 // asin(-1/2) -> -Pi/6
506                 if (x.is_equal(_ex_1_2))
507                         return numeric(-1,6)*Pi;
508                 // asin(-1) -> -Pi/2
509                 if (x.is_equal(_ex_1))
510                         return _num_1_2*Pi;
511                 // asin(float) -> float
512                 if (!x.info(info_flags::crational))
513                         return asin(ex_to<numeric>(x));
514         }
515         
516         return asin(x).hold();
517 }
518
519 static ex asin_deriv(const ex & x, unsigned deriv_param)
520 {
521         GINAC_ASSERT(deriv_param==0);
522         
523         // d/dx asin(x) -> 1/sqrt(1-x^2)
524         return power(1-power(x,_ex2),_ex_1_2);
525 }
526
527 REGISTER_FUNCTION(asin, eval_func(asin_eval).
528                         evalf_func(asin_evalf).
529                         derivative_func(asin_deriv).
530                         latex_name("\\arcsin"));
531
532 //////////
533 // inverse cosine (arc cosine)
534 //////////
535
536 static ex acos_evalf(const ex & x)
537 {
538         if (is_exactly_a<numeric>(x))
539                 return acos(ex_to<numeric>(x));
540         
541         return acos(x).hold();
542 }
543
544 static ex acos_eval(const ex & x)
545 {
546         if (x.info(info_flags::numeric)) {
547                 // acos(1) -> 0
548                 if (x.is_equal(_ex1))
549                         return _ex0;
550                 // acos(1/2) -> Pi/3
551                 if (x.is_equal(_ex1_2))
552                         return _ex1_3*Pi;
553                 // acos(0) -> Pi/2
554                 if (x.is_zero())
555                         return _ex1_2*Pi;
556                 // acos(-1/2) -> 2/3*Pi
557                 if (x.is_equal(_ex_1_2))
558                         return numeric(2,3)*Pi;
559                 // acos(-1) -> Pi
560                 if (x.is_equal(_ex_1))
561                         return Pi;
562                 // acos(float) -> float
563                 if (!x.info(info_flags::crational))
564                         return acos(ex_to<numeric>(x));
565         }
566         
567         return acos(x).hold();
568 }
569
570 static ex acos_deriv(const ex & x, unsigned deriv_param)
571 {
572         GINAC_ASSERT(deriv_param==0);
573         
574         // d/dx acos(x) -> -1/sqrt(1-x^2)
575         return -power(1-power(x,_ex2),_ex_1_2);
576 }
577
578 REGISTER_FUNCTION(acos, eval_func(acos_eval).
579                         evalf_func(acos_evalf).
580                         derivative_func(acos_deriv).
581                         latex_name("\\arccos"));
582
583 //////////
584 // inverse tangent (arc tangent)
585 //////////
586
587 static ex atan_evalf(const ex & x)
588 {
589         if (is_exactly_a<numeric>(x))
590                 return atan(ex_to<numeric>(x));
591         
592         return atan(x).hold();
593 }
594
595 static ex atan_eval(const ex & x)
596 {
597         if (x.info(info_flags::numeric)) {
598                 // atan(0) -> 0
599                 if (x.is_zero())
600                         return _ex0;
601                 // atan(1) -> Pi/4
602                 if (x.is_equal(_ex1))
603                         return _ex1_4*Pi;
604                 // atan(-1) -> -Pi/4
605                 if (x.is_equal(_ex_1))
606                         return _ex_1_4*Pi;
607                 if (x.is_equal(I) || x.is_equal(-I))
608                         throw (pole_error("atan_eval(): logarithmic pole",0));
609                 // atan(float) -> float
610                 if (!x.info(info_flags::crational))
611                         return atan(ex_to<numeric>(x));
612         }
613         
614         return atan(x).hold();
615 }
616
617 static ex atan_deriv(const ex & x, unsigned deriv_param)
618 {
619         GINAC_ASSERT(deriv_param==0);
620
621         // d/dx atan(x) -> 1/(1+x^2)
622         return power(_ex1+power(x,_ex2), _ex_1);
623 }
624
625 static ex atan_series(const ex &arg,
626                       const relational &rel,
627                       int order,
628                       unsigned options)
629 {
630         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
631         // method:
632         // Taylor series where there is no pole or cut falls back to atan_deriv.
633         // There are two branch cuts, one runnig from I up the imaginary axis and
634         // one running from -I down the imaginary axis.  The points I and -I are
635         // poles.
636         // On the branch cuts and the poles series expand
637         //     (log(1+I*x)-log(1-I*x))/(2*I)
638         // instead.
639         const ex arg_pt = arg.subs(rel);
640         if (!(I*arg_pt).info(info_flags::real))
641                 throw do_taylor();     // Re(x) != 0
642         if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1)
643                 throw do_taylor();     // Re(x) == 0, but abs(x)<1
644         // care for the poles, using the defining formula for atan()...
645         if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
646                 return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
647         if (!(options & series_options::suppress_branchcut)) {
648                 // method:
649                 // This is the branch cut: assemble the primitive series manually and
650                 // then add the corresponding complex step function.
651                 const symbol &s = ex_to<symbol>(rel.lhs());
652                 const ex &point = rel.rhs();
653                 const symbol foo;
654                 const ex replarg = series(atan(arg), s==foo, order).subs(foo==point);
655                 ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2;
656                 if ((I*arg_pt)<_ex0)
657                         Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2;
658                 else
659                         Order0correction += log((I*arg_pt+_ex1)/(I*arg_pt+_ex_1))*I*_ex1_2;
660                 epvector seq;
661                 seq.push_back(expair(Order0correction, _ex0));
662                 seq.push_back(expair(Order(_ex1), order));
663                 return series(replarg - pseries(rel, seq), rel, order);
664         }
665         throw do_taylor();
666 }
667
668 REGISTER_FUNCTION(atan, eval_func(atan_eval).
669                         evalf_func(atan_evalf).
670                         derivative_func(atan_deriv).
671                         series_func(atan_series).
672                         latex_name("\\arctan"));
673
674 //////////
675 // inverse tangent (atan2(y,x))
676 //////////
677
678 static ex atan2_evalf(const ex &y, const ex &x)
679 {
680         if (is_exactly_a<numeric>(y) && is_exactly_a<numeric>(x))
681                 return atan2(ex_to<numeric>(y), ex_to<numeric>(x));
682         
683         return atan2(y, x).hold();
684 }
685
686 static ex atan2_eval(const ex & y, const ex & x)
687 {
688         if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
689                 x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
690                 return atan2_evalf(y,x);
691         }
692         
693         return atan2(y,x).hold();
694 }    
695
696 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
697 {
698         GINAC_ASSERT(deriv_param<2);
699         
700         if (deriv_param==0) {
701                 // d/dy atan(y,x)
702                 return x*power(power(x,_ex2)+power(y,_ex2),_ex_1);
703         }
704         // d/dx atan(y,x)
705         return -y*power(power(x,_ex2)+power(y,_ex2),_ex_1);
706 }
707
708 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
709                          evalf_func(atan2_evalf).
710                          derivative_func(atan2_deriv));
711
712 //////////
713 // hyperbolic sine (trigonometric function)
714 //////////
715
716 static ex sinh_evalf(const ex & x)
717 {
718         if (is_exactly_a<numeric>(x))
719                 return sinh(ex_to<numeric>(x));
720         
721         return sinh(x).hold();
722 }
723
724 static ex sinh_eval(const ex & x)
725 {
726         if (x.info(info_flags::numeric)) {
727                 if (x.is_zero())  // sinh(0) -> 0
728                         return _ex0;        
729                 if (!x.info(info_flags::crational))  // sinh(float) -> float
730                         return sinh(ex_to<numeric>(x));
731         }
732         
733         if ((x/Pi).info(info_flags::numeric) &&
734                 ex_to<numeric>(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
735                 return I*sin(x/I);
736         
737         if (is_exactly_a<function>(x)) {
738                 const ex &t = x.op(0);
739                 // sinh(asinh(x)) -> x
740                 if (is_ex_the_function(x, asinh))
741                         return t;
742                 // sinh(acosh(x)) -> sqrt(x-1) * sqrt(x+1)
743                 if (is_ex_the_function(x, acosh))
744                         return sqrt(t-_ex1)*sqrt(t+_ex1);
745                 // sinh(atanh(x)) -> x/sqrt(1-x^2)
746                 if (is_ex_the_function(x, atanh))
747                         return t*power(_ex1-power(t,_ex2),_ex_1_2);
748         }
749         
750         return sinh(x).hold();
751 }
752
753 static ex sinh_deriv(const ex & x, unsigned deriv_param)
754 {
755         GINAC_ASSERT(deriv_param==0);
756         
757         // d/dx sinh(x) -> cosh(x)
758         return cosh(x);
759 }
760
761 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
762                         evalf_func(sinh_evalf).
763                         derivative_func(sinh_deriv).
764                         latex_name("\\sinh"));
765
766 //////////
767 // hyperbolic cosine (trigonometric function)
768 //////////
769
770 static ex cosh_evalf(const ex & x)
771 {
772         if (is_exactly_a<numeric>(x))
773                 return cosh(ex_to<numeric>(x));
774         
775         return cosh(x).hold();
776 }
777
778 static ex cosh_eval(const ex & x)
779 {
780         if (x.info(info_flags::numeric)) {
781                 if (x.is_zero())  // cosh(0) -> 1
782                         return _ex1;
783                 if (!x.info(info_flags::crational))  // cosh(float) -> float
784                         return cosh(ex_to<numeric>(x));
785         }
786         
787         if ((x/Pi).info(info_flags::numeric) &&
788                 ex_to<numeric>(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
789                 return cos(x/I);
790         
791         if (is_exactly_a<function>(x)) {
792                 const ex &t = x.op(0);
793                 // cosh(acosh(x)) -> x
794                 if (is_ex_the_function(x, acosh))
795                         return t;
796                 // cosh(asinh(x)) -> sqrt(1+x^2)
797                 if (is_ex_the_function(x, asinh))
798                         return sqrt(_ex1+power(t,_ex2));
799                 // cosh(atanh(x)) -> 1/sqrt(1-x^2)
800                 if (is_ex_the_function(x, atanh))
801                         return power(_ex1-power(t,_ex2),_ex_1_2);
802         }
803         
804         return cosh(x).hold();
805 }
806
807 static ex cosh_deriv(const ex & x, unsigned deriv_param)
808 {
809         GINAC_ASSERT(deriv_param==0);
810         
811         // d/dx cosh(x) -> sinh(x)
812         return sinh(x);
813 }
814
815 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
816                         evalf_func(cosh_evalf).
817                         derivative_func(cosh_deriv).
818                         latex_name("\\cosh"));
819
820 //////////
821 // hyperbolic tangent (trigonometric function)
822 //////////
823
824 static ex tanh_evalf(const ex & x)
825 {
826         if (is_exactly_a<numeric>(x))
827                 return tanh(ex_to<numeric>(x));
828         
829         return tanh(x).hold();
830 }
831
832 static ex tanh_eval(const ex & x)
833 {
834         if (x.info(info_flags::numeric)) {
835                 if (x.is_zero())  // tanh(0) -> 0
836                         return _ex0;
837                 if (!x.info(info_flags::crational))  // tanh(float) -> float
838                         return tanh(ex_to<numeric>(x));
839         }
840         
841         if ((x/Pi).info(info_flags::numeric) &&
842                 ex_to<numeric>(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
843                 return I*tan(x/I);
844         
845         if (is_exactly_a<function>(x)) {
846                 const ex &t = x.op(0);
847                 // tanh(atanh(x)) -> x
848                 if (is_ex_the_function(x, atanh))
849                         return t;
850                 // tanh(asinh(x)) -> x/sqrt(1+x^2)
851                 if (is_ex_the_function(x, asinh))
852                         return t*power(_ex1+power(t,_ex2),_ex_1_2);
853                 // tanh(acosh(x)) -> sqrt(x-1)*sqrt(x+1)/x
854                 if (is_ex_the_function(x, acosh))
855                         return sqrt(t-_ex1)*sqrt(t+_ex1)*power(t,_ex_1);
856         }
857         
858         return tanh(x).hold();
859 }
860
861 static ex tanh_deriv(const ex & x, unsigned deriv_param)
862 {
863         GINAC_ASSERT(deriv_param==0);
864         
865         // d/dx tanh(x) -> 1-tanh(x)^2
866         return _ex1-power(tanh(x),_ex2);
867 }
868
869 static ex tanh_series(const ex &x,
870                       const relational &rel,
871                       int order,
872                       unsigned options)
873 {
874         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
875         // method:
876         // Taylor series where there is no pole falls back to tanh_deriv.
877         // On a pole simply expand sinh(x)/cosh(x).
878         const ex x_pt = x.subs(rel);
879         if (!(2*I*x_pt/Pi).info(info_flags::odd))
880                 throw do_taylor();  // caught by function::series()
881         // if we got here we have to care for a simple pole
882         return (sinh(x)/cosh(x)).series(rel, order+2, options);
883 }
884
885 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
886                         evalf_func(tanh_evalf).
887                         derivative_func(tanh_deriv).
888                         series_func(tanh_series).
889                         latex_name("\\tanh"));
890
891 //////////
892 // inverse hyperbolic sine (trigonometric function)
893 //////////
894
895 static ex asinh_evalf(const ex & x)
896 {
897         if (is_exactly_a<numeric>(x))
898                 return asinh(ex_to<numeric>(x));
899         
900         return asinh(x).hold();
901 }
902
903 static ex asinh_eval(const ex & x)
904 {
905         if (x.info(info_flags::numeric)) {
906                 // asinh(0) -> 0
907                 if (x.is_zero())
908                         return _ex0;
909                 // asinh(float) -> float
910                 if (!x.info(info_flags::crational))
911                         return asinh(ex_to<numeric>(x));
912         }
913         
914         return asinh(x).hold();
915 }
916
917 static ex asinh_deriv(const ex & x, unsigned deriv_param)
918 {
919         GINAC_ASSERT(deriv_param==0);
920         
921         // d/dx asinh(x) -> 1/sqrt(1+x^2)
922         return power(_ex1+power(x,_ex2),_ex_1_2);
923 }
924
925 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
926                          evalf_func(asinh_evalf).
927                          derivative_func(asinh_deriv));
928
929 //////////
930 // inverse hyperbolic cosine (trigonometric function)
931 //////////
932
933 static ex acosh_evalf(const ex & x)
934 {
935         if (is_exactly_a<numeric>(x))
936                 return acosh(ex_to<numeric>(x));
937         
938         return acosh(x).hold();
939 }
940
941 static ex acosh_eval(const ex & x)
942 {
943         if (x.info(info_flags::numeric)) {
944                 // acosh(0) -> Pi*I/2
945                 if (x.is_zero())
946                         return Pi*I*numeric(1,2);
947                 // acosh(1) -> 0
948                 if (x.is_equal(_ex1))
949                         return _ex0;
950                 // acosh(-1) -> Pi*I
951                 if (x.is_equal(_ex_1))
952                         return Pi*I;
953                 // acosh(float) -> float
954                 if (!x.info(info_flags::crational))
955                         return acosh(ex_to<numeric>(x));
956         }
957         
958         return acosh(x).hold();
959 }
960
961 static ex acosh_deriv(const ex & x, unsigned deriv_param)
962 {
963         GINAC_ASSERT(deriv_param==0);
964         
965         // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
966         return power(x+_ex_1,_ex_1_2)*power(x+_ex1,_ex_1_2);
967 }
968
969 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
970                          evalf_func(acosh_evalf).
971                          derivative_func(acosh_deriv));
972
973 //////////
974 // inverse hyperbolic tangent (trigonometric function)
975 //////////
976
977 static ex atanh_evalf(const ex & x)
978 {
979         if (is_exactly_a<numeric>(x))
980                 return atanh(ex_to<numeric>(x));
981         
982         return atanh(x).hold();
983 }
984
985 static ex atanh_eval(const ex & x)
986 {
987         if (x.info(info_flags::numeric)) {
988                 // atanh(0) -> 0
989                 if (x.is_zero())
990                         return _ex0;
991                 // atanh({+|-}1) -> throw
992                 if (x.is_equal(_ex1) || x.is_equal(_ex_1))
993                         throw (pole_error("atanh_eval(): logarithmic pole",0));
994                 // atanh(float) -> float
995                 if (!x.info(info_flags::crational))
996                         return atanh(ex_to<numeric>(x));
997         }
998         
999         return atanh(x).hold();
1000 }
1001
1002 static ex atanh_deriv(const ex & x, unsigned deriv_param)
1003 {
1004         GINAC_ASSERT(deriv_param==0);
1005         
1006         // d/dx atanh(x) -> 1/(1-x^2)
1007         return power(_ex1-power(x,_ex2),_ex_1);
1008 }
1009
1010 static ex atanh_series(const ex &arg,
1011                        const relational &rel,
1012                        int order,
1013                        unsigned options)
1014 {
1015         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
1016         // method:
1017         // Taylor series where there is no pole or cut falls back to atanh_deriv.
1018         // There are two branch cuts, one runnig from 1 up the real axis and one
1019         // one running from -1 down the real axis.  The points 1 and -1 are poles
1020         // On the branch cuts and the poles series expand
1021         //     (log(1+x)-log(1-x))/2
1022         // instead.
1023         const ex arg_pt = arg.subs(rel);
1024         if (!(arg_pt).info(info_flags::real))
1025                 throw do_taylor();     // Im(x) != 0
1026         if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1)
1027                 throw do_taylor();     // Im(x) == 0, but abs(x)<1
1028         // care for the poles, using the defining formula for atanh()...
1029         if (arg_pt.is_equal(_ex1) || arg_pt.is_equal(_ex_1))
1030                 return ((log(_ex1+arg)-log(_ex1-arg))*_ex1_2).series(rel, order, options);
1031         // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
1032         if (!(options & series_options::suppress_branchcut)) {
1033                 // method:
1034                 // This is the branch cut: assemble the primitive series manually and
1035                 // then add the corresponding complex step function.
1036                 const symbol &s = ex_to<symbol>(rel.lhs());
1037                 const ex &point = rel.rhs();
1038                 const symbol foo;
1039                 const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point);
1040                 ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2;
1041                 if (arg_pt<_ex0)
1042                         Order0correction += log((arg_pt+_ex_1)/(arg_pt+_ex1))*_ex1_2;
1043                 else
1044                         Order0correction += log((arg_pt+_ex1)/(arg_pt+_ex_1))*_ex_1_2;
1045                 epvector seq;
1046                 seq.push_back(expair(Order0correction, _ex0));
1047                 seq.push_back(expair(Order(_ex1), order));
1048                 return series(replarg - pseries(rel, seq), rel, order);
1049         }
1050         throw do_taylor();
1051 }
1052
1053 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
1054                          evalf_func(atanh_evalf).
1055                          derivative_func(atanh_deriv).
1056                          series_func(atanh_series));
1057
1058
1059 } // namespace GiNaC