- switched to automake build environment
[ginac.git] / ginac / inifcns_trans.cpp
1 /** @file inifcns_trans.cpp
2  *
3  *  Implementation of transcendental (and trigonometric and hyperbolic)
4  *  functions. */
5
6 #include <vector>
7 #include <stdexcept>
8
9 #include "ginac.h"
10
11 //////////
12 // exponential function
13 //////////
14
15 ex exp_evalf(ex const & x)
16 {
17     BEGIN_TYPECHECK
18         TYPECHECK(x,numeric)
19     END_TYPECHECK(exp(x))
20     
21     return exp(ex_to_numeric(x)); // -> numeric exp(numeric)
22 }
23
24 ex exp_eval(ex const & x)
25 {
26     // exp(0) -> 1
27     if (x.is_zero()) {
28         return exONE();
29     }
30     // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
31     ex TwoExOverPiI=(2*x)/(Pi*I);
32     if (TwoExOverPiI.info(info_flags::integer)) {
33         numeric z=mod(ex_to_numeric(TwoExOverPiI),numeric(4));
34         if (z.is_equal(numZERO()))
35             return exONE();
36         if (z.is_equal(numONE()))
37             return ex(I);
38         if (z.is_equal(numTWO()))
39             return exMINUSONE();
40         if (z.is_equal(numTHREE()))
41             return ex(-I);
42     }
43     // exp(log(x)) -> x
44     if (is_ex_the_function(x, log))
45         return x.op(0);
46     
47     // exp(float)
48     if (x.info(info_flags::numeric) && !x.info(info_flags::rational))
49         return exp_evalf(x);
50     
51     return exp(x).hold();
52 }    
53
54 ex exp_diff(ex const & x, unsigned diff_param)
55 {
56     ASSERT(diff_param==0);
57
58     return exp(x);
59 }
60
61 REGISTER_FUNCTION(exp, exp_eval, exp_evalf, exp_diff, NULL);
62
63 //////////
64 // natural logarithm
65 //////////
66
67 ex log_evalf(ex const & x)
68 {
69     BEGIN_TYPECHECK
70         TYPECHECK(x,numeric)
71     END_TYPECHECK(log(x))
72     
73     return log(ex_to_numeric(x)); // -> numeric log(numeric)
74 }
75
76 ex log_eval(ex const & x)
77 {
78     if (x.info(info_flags::numeric)) {
79         // log(1) -> 0
80         if (x.is_equal(exONE()))
81             return exZERO();
82         // log(-1) -> I*Pi
83         if (x.is_equal(exMINUSONE()))
84             return (I*Pi);
85         // log(I) -> Pi*I/2
86         if (x.is_equal(I))
87             return (I*Pi*numeric(1,2));
88         // log(-I) -> -Pi*I/2
89         if (x.is_equal(-I))
90             return (I*Pi*numeric(-1,2));
91         // log(0) -> throw singularity
92         if (x.is_equal(exZERO()))
93             throw(std::domain_error("log_eval(): log(0)"));
94         // log(float)
95         if (!x.info(info_flags::rational))
96             return log_evalf(x);
97     }
98     
99     return log(x).hold();
100 }    
101
102 ex log_diff(ex const & x, unsigned diff_param)
103 {
104     ASSERT(diff_param==0);
105
106     return power(x, -1);
107 }
108
109 REGISTER_FUNCTION(log, log_eval, log_evalf, log_diff, NULL);
110
111 //////////
112 // sine (trigonometric function)
113 //////////
114
115 ex sin_evalf(ex const & x)
116 {
117     BEGIN_TYPECHECK
118        TYPECHECK(x,numeric)
119     END_TYPECHECK(sin(x))
120     
121     return sin(ex_to_numeric(x)); // -> numeric sin(numeric)
122 }
123
124 ex sin_eval(ex const & x)
125 {
126     // sin(n*Pi) -> 0
127     ex xOverPi=x/Pi;
128     if (xOverPi.info(info_flags::integer))
129         return exZERO();
130     
131     // sin((2n+1)*Pi/2) -> {+|-}1
132     ex xOverPiMinusHalf=xOverPi-exHALF();
133     if (xOverPiMinusHalf.info(info_flags::even))
134         return exONE();
135     else if (xOverPiMinusHalf.info(info_flags::odd))
136         return exMINUSONE();
137     
138     if (is_ex_exactly_of_type(x, function)) {
139         ex t=x.op(0);
140         // sin(asin(x)) -> x
141         if (is_ex_the_function(x, asin))
142             return t;
143         // sin(acos(x)) -> (1-x^2)^(1/2)
144         if (is_ex_the_function(x, acos))
145             return power(exONE()-power(t,exTWO()),exHALF());
146         // sin(atan(x)) -> x*(1+x^2)^(-1/2)
147         if (is_ex_the_function(x, atan))
148             return t*power(exONE()+power(t,exTWO()),exMINUSHALF());
149     }
150     
151     // sin(float) -> float
152     if (x.info(info_flags::numeric) && !x.info(info_flags::rational))
153         return sin_evalf(x);
154     
155     return sin(x).hold();
156 }
157
158 ex sin_diff(ex const & x, unsigned diff_param)
159 {
160     ASSERT(diff_param==0);
161     
162     return cos(x);
163 }
164
165 REGISTER_FUNCTION(sin, sin_eval, sin_evalf, sin_diff, NULL);
166
167 //////////
168 // cosine (trigonometric function)
169 //////////
170
171 ex cos_evalf(ex const & x)
172 {
173     BEGIN_TYPECHECK
174         TYPECHECK(x,numeric)
175     END_TYPECHECK(cos(x))
176     
177     return cos(ex_to_numeric(x)); // -> numeric cos(numeric)
178 }
179
180 ex cos_eval(ex const & x)
181 {
182     // cos(n*Pi) -> {+|-}1
183     ex xOverPi=x/Pi;
184     if (xOverPi.info(info_flags::even))
185         return exONE();
186     else if (xOverPi.info(info_flags::odd))
187         return exMINUSONE();
188     
189     // cos((2n+1)*Pi/2) -> 0
190     ex xOverPiMinusHalf=xOverPi-exHALF();
191     if (xOverPiMinusHalf.info(info_flags::integer))
192         return exZERO();
193     
194     if (is_ex_exactly_of_type(x, function)) {
195         ex t=x.op(0);
196         // cos(acos(x)) -> x
197         if (is_ex_the_function(x, acos))
198             return t;
199         // cos(asin(x)) -> (1-x^2)^(1/2)
200         if (is_ex_the_function(x, asin))
201             return power(exONE()-power(t,exTWO()),exHALF());
202         // cos(atan(x)) -> (1+x^2)^(-1/2)
203         if (is_ex_the_function(x, atan))
204             return power(exONE()+power(t,exTWO()),exMINUSHALF());
205     }
206     
207     // cos(float) -> float
208     if (x.info(info_flags::numeric) && !x.info(info_flags::rational))
209         return cos_evalf(x);
210     
211     return cos(x).hold();
212 }
213
214 ex cos_diff(ex const & x, unsigned diff_param)
215 {
216     ASSERT(diff_param==0);
217
218     return numMINUSONE()*sin(x);
219 }
220
221 REGISTER_FUNCTION(cos, cos_eval, cos_evalf, cos_diff, NULL);
222
223 //////////
224 // tangent (trigonometric function)
225 //////////
226
227 ex tan_evalf(ex const & x)
228 {
229     BEGIN_TYPECHECK
230        TYPECHECK(x,numeric)
231     END_TYPECHECK(tan(x)) // -> numeric tan(numeric)
232     
233     return tan(ex_to_numeric(x));
234 }
235
236 ex tan_eval(ex const & x)
237 {
238     // tan(n*Pi/3) -> {0|3^(1/2)|-(3^(1/2))}
239     ex ThreeExOverPi=numTHREE()*x/Pi;
240     if (ThreeExOverPi.info(info_flags::integer)) {
241         numeric z=mod(ex_to_numeric(ThreeExOverPi),numeric(3));
242         if (z.is_equal(numZERO()))
243             return exZERO();
244         if (z.is_equal(numONE()))
245             return power(exTHREE(),exHALF());
246         if (z.is_equal(numTWO()))
247             return -power(exTHREE(),exHALF());
248     }
249     
250     // tan((2n+1)*Pi/2) -> throw
251     ex ExOverPiMinusHalf=x/Pi-exHALF();
252     if (ExOverPiMinusHalf.info(info_flags::integer))
253         throw (std::domain_error("tan_eval(): infinity"));
254     
255     if (is_ex_exactly_of_type(x, function)) {
256         ex t=x.op(0);
257         // tan(atan(x)) -> x
258         if (is_ex_the_function(x, atan))
259             return t;
260         // tan(asin(x)) -> x*(1+x^2)^(-1/2)
261         if (is_ex_the_function(x, asin))
262             return t*power(exONE()-power(t,exTWO()),exMINUSHALF());
263         // tan(acos(x)) -> (1-x^2)^(1/2)/x
264         if (is_ex_the_function(x, acos))
265             return power(t,exMINUSONE())*power(exONE()-power(t,exTWO()),exHALF());
266     }
267     
268     // tan(float) -> float
269     if (x.info(info_flags::numeric) && !x.info(info_flags::rational)) {
270         return tan_evalf(x);
271     }
272     
273     return tan(x).hold();
274 }
275
276 ex tan_diff(ex const & x, unsigned diff_param)
277 {
278     ASSERT(diff_param==0);
279     
280     return (1+power(tan(x),exTWO()));
281 }
282
283 REGISTER_FUNCTION(tan, tan_eval, tan_evalf, tan_diff, NULL);
284
285 //////////
286 // inverse sine (arc sine)
287 //////////
288
289 ex asin_evalf(ex const & x)
290 {
291     BEGIN_TYPECHECK
292        TYPECHECK(x,numeric)
293     END_TYPECHECK(asin(x))
294     
295     return asin(ex_to_numeric(x)); // -> numeric asin(numeric)
296 }
297
298 ex asin_eval(ex const & x)
299 {
300     if (x.info(info_flags::numeric)) {
301         // asin(0) -> 0
302         if (x.is_zero())
303             return x;
304         // asin(1/2) -> Pi/6
305         if (x.is_equal(exHALF()))
306             return numeric(1,6)*Pi;
307         // asin(1) -> Pi/2
308         if (x.is_equal(exONE()))
309             return numeric(1,2)*Pi;
310         // asin(-1/2) -> -Pi/6
311         if (x.is_equal(exMINUSHALF()))
312             return numeric(-1,6)*Pi;
313         // asin(-1) -> -Pi/2
314         if (x.is_equal(exMINUSONE()))
315             return numeric(-1,2)*Pi;
316         // asin(float) -> float
317         if (!x.info(info_flags::rational))
318             return asin_evalf(x);
319     }
320     
321     return asin(x).hold();
322 }
323
324 ex asin_diff(ex const & x, unsigned diff_param)
325 {
326     ASSERT(diff_param==0);
327     
328     return power(1-power(x,exTWO()),exMINUSHALF());
329 }
330
331 REGISTER_FUNCTION(asin, asin_eval, asin_evalf, asin_diff, NULL);
332
333 //////////
334 // inverse cosine (arc cosine)
335 //////////
336
337 ex acos_evalf(ex const & x)
338 {
339     BEGIN_TYPECHECK
340        TYPECHECK(x,numeric)
341     END_TYPECHECK(acos(x))
342     
343     return acos(ex_to_numeric(x)); // -> numeric acos(numeric)
344 }
345
346 ex acos_eval(ex const & x)
347 {
348     if (x.info(info_flags::numeric)) {
349         // acos(1) -> 0
350         if (x.is_equal(exONE()))
351             return exZERO();
352         // acos(1/2) -> Pi/3
353         if (x.is_equal(exHALF()))
354             return numeric(1,3)*Pi;
355         // acos(0) -> Pi/2
356         if (x.is_zero())
357             return numeric(1,2)*Pi;
358         // acos(-1/2) -> 2/3*Pi
359         if (x.is_equal(exMINUSHALF()))
360             return numeric(2,3)*Pi;
361         // acos(-1) -> Pi
362         if (x.is_equal(exMINUSONE()))
363             return Pi;
364         // acos(float) -> float
365         if (!x.info(info_flags::rational))
366             return acos_evalf(x);
367     }
368     
369     return acos(x).hold();
370 }
371
372 ex acos_diff(ex const & x, unsigned diff_param)
373 {
374     ASSERT(diff_param==0);
375     
376     return numMINUSONE()*power(1-power(x,exTWO()),exMINUSHALF());
377 }
378
379 REGISTER_FUNCTION(acos, acos_eval, acos_evalf, acos_diff, NULL);
380
381 //////////
382 // inverse tangent (arc tangent)
383 //////////
384
385 ex atan_evalf(ex const & x)
386 {
387     BEGIN_TYPECHECK
388         TYPECHECK(x,numeric)
389     END_TYPECHECK(atan(x))
390     
391     return atan(ex_to_numeric(x)); // -> numeric atan(numeric)
392 }
393
394 ex atan_eval(ex const & x)
395 {
396     if (x.info(info_flags::numeric)) {
397         // atan(0) -> 0
398         if (x.is_equal(exZERO()))
399             return exZERO();
400         // atan(float) -> float
401         if (!x.info(info_flags::rational))
402             return atan_evalf(x);
403     }
404     
405     return atan(x).hold();
406 }    
407
408 ex atan_diff(ex const & x, unsigned diff_param)
409 {
410     ASSERT(diff_param==0);
411
412     return power(1+x*x, -1);
413 }
414
415 REGISTER_FUNCTION(atan, atan_eval, atan_evalf, atan_diff, NULL);
416
417 //////////
418 // inverse tangent (atan2(y,x))
419 //////////
420
421 ex atan2_evalf(ex const & y, ex const & x)
422 {
423     BEGIN_TYPECHECK
424         TYPECHECK(y,numeric)
425         TYPECHECK(x,numeric)
426     END_TYPECHECK(atan2(y,x))
427     
428     return atan(ex_to_numeric(y),ex_to_numeric(x)); // -> numeric atan(numeric)
429 }
430
431 ex atan2_eval(ex const & y, ex const & x)
432 {
433     if (y.info(info_flags::numeric) && !y.info(info_flags::rational) &&
434         x.info(info_flags::numeric) && !x.info(info_flags::rational)) {
435         return atan2_evalf(y,x);
436     }
437     
438     return atan2(y,x).hold();
439 }    
440
441 ex atan2_diff(ex const & y, ex const & x, unsigned diff_param)
442 {
443     ASSERT(diff_param<2);
444
445     if (diff_param==0) {
446         // d/dy atan(y,x)
447         return power(x*(1+y*y/(x*x)),-1);
448     }
449     // d/dx atan(y,x)
450     return -y*power(x*x+y*y,-1);
451 }
452
453 REGISTER_FUNCTION(atan2, atan2_eval, atan2_evalf, atan2_diff, NULL);
454
455 //////////
456 // hyperbolic sine (trigonometric function)
457 //////////
458
459 ex sinh_evalf(ex const & x)
460 {
461     BEGIN_TYPECHECK
462        TYPECHECK(x,numeric)
463     END_TYPECHECK(sinh(x))
464     
465     return sinh(ex_to_numeric(x)); // -> numeric sinh(numeric)
466 }
467
468 ex sinh_eval(ex const & x)
469 {
470     if (x.info(info_flags::numeric)) {
471         // sinh(0) -> 0
472         if (x.is_zero())
473             return exZERO();
474         // sinh(float) -> float
475         if (!x.info(info_flags::rational))
476             return sinh_evalf(x);
477     }
478     
479     if (is_ex_exactly_of_type(x, function)) {
480         ex t=x.op(0);
481         // sinh(asinh(x)) -> x
482         if (is_ex_the_function(x, asinh))
483             return t;
484         // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2)
485         if (is_ex_the_function(x, acosh))
486             return power(t-exONE(),exHALF())*power(t+exONE(),exHALF());
487         // sinh(atanh(x)) -> x*(1-x^2)^(-1/2)
488         if (is_ex_the_function(x, atanh))
489             return t*power(exONE()-power(t,exTWO()),exMINUSHALF());
490     }
491     
492     return sinh(x).hold();
493 }
494
495 ex sinh_diff(ex const & x, unsigned diff_param)
496 {
497     ASSERT(diff_param==0);
498     
499     return cosh(x);
500 }
501
502 REGISTER_FUNCTION(sinh, sinh_eval, sinh_evalf, sinh_diff, NULL);
503
504 //////////
505 // hyperbolic cosine (trigonometric function)
506 //////////
507
508 ex cosh_evalf(ex const & x)
509 {
510     BEGIN_TYPECHECK
511        TYPECHECK(x,numeric)
512     END_TYPECHECK(cosh(x))
513     
514     return cosh(ex_to_numeric(x)); // -> numeric cosh(numeric)
515 }
516
517 ex cosh_eval(ex const & x)
518 {
519     if (x.info(info_flags::numeric)) {
520         // cosh(0) -> 1
521         if (x.is_zero())
522             return exONE();
523         // cosh(float) -> float
524         if (!x.info(info_flags::rational))
525             return cosh_evalf(x);
526     }
527     
528     if (is_ex_exactly_of_type(x, function)) {
529         ex t=x.op(0);
530         // cosh(acosh(x)) -> x
531         if (is_ex_the_function(x, acosh))
532             return t;
533         // cosh(asinh(x)) -> (1+x^2)^(1/2)
534         if (is_ex_the_function(x, asinh))
535             return power(exONE()+power(t,exTWO()),exHALF());
536         // cosh(atanh(x)) -> (1-x^2)^(-1/2)
537         if (is_ex_the_function(x, atanh))
538             return power(exONE()-power(t,exTWO()),exMINUSHALF());
539     }
540     
541     return cosh(x).hold();
542 }
543
544 ex cosh_diff(ex const & x, unsigned diff_param)
545 {
546     ASSERT(diff_param==0);
547     
548     return sinh(x);
549 }
550
551 REGISTER_FUNCTION(cosh, cosh_eval, cosh_evalf, cosh_diff, NULL);
552
553 //////////
554 // hyperbolic tangent (trigonometric function)
555 //////////
556
557 ex tanh_evalf(ex const & x)
558 {
559     BEGIN_TYPECHECK
560        TYPECHECK(x,numeric)
561     END_TYPECHECK(tanh(x))
562     
563     return tanh(ex_to_numeric(x)); // -> numeric tanh(numeric)
564 }
565
566 ex tanh_eval(ex const & x)
567 {
568     if (x.info(info_flags::numeric)) {
569         // tanh(0) -> 0
570         if (x.is_zero())
571             return exZERO();
572         // tanh(float) -> float
573         if (!x.info(info_flags::rational))
574             return tanh_evalf(x);
575     }
576     
577     if (is_ex_exactly_of_type(x, function)) {
578         ex t=x.op(0);
579         // tanh(atanh(x)) -> x
580         if (is_ex_the_function(x, atanh))
581             return t;
582         // tanh(asinh(x)) -> x*(1+x^2)^(-1/2)
583         if (is_ex_the_function(x, asinh))
584             return t*power(exONE()+power(t,exTWO()),exMINUSHALF());
585         // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x
586         if (is_ex_the_function(x, acosh))
587             return power(t-exONE(),exHALF())*power(t+exONE(),exHALF())*power(t,exMINUSONE());
588     }
589     
590     return tanh(x).hold();
591 }
592
593 ex tanh_diff(ex const & x, unsigned diff_param)
594 {
595     ASSERT(diff_param==0);
596     
597     return exONE()-power(tanh(x),exTWO());
598 }
599
600 REGISTER_FUNCTION(tanh, tanh_eval, tanh_evalf, tanh_diff, NULL);
601
602 //////////
603 // inverse hyperbolic sine (trigonometric function)
604 //////////
605
606 ex asinh_evalf(ex const & x)
607 {
608     BEGIN_TYPECHECK
609        TYPECHECK(x,numeric)
610     END_TYPECHECK(asinh(x))
611     
612     return asinh(ex_to_numeric(x)); // -> numeric asinh(numeric)
613 }
614
615 ex asinh_eval(ex const & x)
616 {
617     if (x.info(info_flags::numeric)) {
618         // asinh(0) -> 0
619         if (x.is_zero())
620             return exZERO();
621         // asinh(float) -> float
622         if (!x.info(info_flags::rational))
623             return asinh_evalf(x);
624     }
625     
626     return asinh(x).hold();
627 }
628
629 ex asinh_diff(ex const & x, unsigned diff_param)
630 {
631     ASSERT(diff_param==0);
632     
633     return power(1+power(x,exTWO()),exMINUSHALF());
634 }
635
636 REGISTER_FUNCTION(asinh, asinh_eval, asinh_evalf, asinh_diff, NULL);
637
638 //////////
639 // inverse hyperbolic cosine (trigonometric function)
640 //////////
641
642 ex acosh_evalf(ex const & x)
643 {
644     BEGIN_TYPECHECK
645        TYPECHECK(x,numeric)
646     END_TYPECHECK(acosh(x))
647     
648     return acosh(ex_to_numeric(x)); // -> numeric acosh(numeric)
649 }
650
651 ex acosh_eval(ex const & x)
652 {
653     if (x.info(info_flags::numeric)) {
654         // acosh(0) -> Pi*I/2
655         if (x.is_zero())
656             return Pi*I*numeric(1,2);
657         // acosh(1) -> 0
658         if (x.is_equal(exONE()))
659             return exZERO();
660         // acosh(-1) -> Pi*I
661         if (x.is_equal(exMINUSONE()))
662             return Pi*I;
663         // acosh(float) -> float
664         if (!x.info(info_flags::rational))
665             return acosh_evalf(x);
666     }
667     
668     return acosh(x).hold();
669 }
670
671 ex acosh_diff(ex const & x, unsigned diff_param)
672 {
673     ASSERT(diff_param==0);
674     
675     return power(x-1,exMINUSHALF())*power(x+1,exMINUSHALF());
676 }
677
678 REGISTER_FUNCTION(acosh, acosh_eval, acosh_evalf, acosh_diff, NULL);
679
680 //////////
681 // inverse hyperbolic tangent (trigonometric function)
682 //////////
683
684 ex atanh_evalf(ex const & x)
685 {
686     BEGIN_TYPECHECK
687        TYPECHECK(x,numeric)
688     END_TYPECHECK(atanh(x))
689     
690     return atanh(ex_to_numeric(x)); // -> numeric atanh(numeric)
691 }
692
693 ex atanh_eval(ex const & x)
694 {
695     if (x.info(info_flags::numeric)) {
696         // atanh(0) -> 0
697         if (x.is_zero())
698             return exZERO();
699         // atanh({+|-}1) -> throw
700         if (x.is_equal(exONE()) || x.is_equal(exONE()))
701             throw (std::domain_error("atanh_eval(): infinity"));
702         // atanh(float) -> float
703         if (!x.info(info_flags::rational))
704             return atanh_evalf(x);
705     }
706     
707     return atanh(x).hold();
708 }
709
710 ex atanh_diff(ex const & x, unsigned diff_param)
711 {
712     ASSERT(diff_param==0);
713     
714     return power(exONE()-power(x,exTWO()),exMINUSONE());
715 }
716
717 REGISTER_FUNCTION(atanh, atanh_eval, atanh_evalf, atanh_diff, NULL);