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