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