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