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