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