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