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