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