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