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