- Cleanups: My evil plot of making ex::bp private may finally be carried
[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-2001 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 #include "relational.h"
33 #include "symbol.h"
34 #include "pseries.h"
35 #include "utils.h"
36
37 namespace GiNaC {
38
39 //////////
40 // exponential function
41 //////////
42
43 static ex exp_evalf(const ex & x)
44 {
45         if (is_exactly_a<numeric>(x))
46                 return exp(ex_to<numeric>(x));
47         
48         return exp(x).hold();
49 }
50
51 static ex exp_eval(const ex & x)
52 {
53         // exp(0) -> 1
54         if (x.is_zero()) {
55                 return _ex1();
56         }
57         // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
58         const ex TwoExOverPiI=(_ex2()*x)/(Pi*I);
59         if (TwoExOverPiI.info(info_flags::integer)) {
60                 numeric z=mod(ex_to<numeric>(TwoExOverPiI),_num4());
61                 if (z.is_equal(_num0()))
62                         return _ex1();
63                 if (z.is_equal(_num1()))
64                         return ex(I);
65                 if (z.is_equal(_num2()))
66                         return _ex_1();
67                 if (z.is_equal(_num3()))
68                         return ex(-I);
69         }
70         // exp(log(x)) -> x
71         if (is_ex_the_function(x, log))
72                 return x.op(0);
73         
74         // exp(float)
75         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
76                 return exp(ex_to<numeric>(x));
77         
78         return exp(x).hold();
79 }
80
81 static ex exp_deriv(const ex & x, unsigned deriv_param)
82 {
83         GINAC_ASSERT(deriv_param==0);
84
85         // d/dx exp(x) -> exp(x)
86         return exp(x);
87 }
88
89 REGISTER_FUNCTION(exp, eval_func(exp_eval).
90                        evalf_func(exp_evalf).
91                        derivative_func(exp_deriv).
92                        latex_name("\\exp"));
93
94 //////////
95 // natural logarithm
96 //////////
97
98 static ex log_evalf(const ex & x)
99 {
100         if (is_exactly_a<numeric>(x))
101                 return log(ex_to<numeric>(x));
102         
103         return log(x).hold();
104 }
105
106 static ex log_eval(const ex & x)
107 {
108         if (x.info(info_flags::numeric)) {
109                 if (x.is_zero())         // log(0) -> infinity
110                         throw(pole_error("log_eval(): log(0)",0));
111                 if (x.info(info_flags::real) && x.info(info_flags::negative))
112                         return (log(-x)+I*Pi);
113                 if (x.is_equal(_ex1()))  // log(1) -> 0
114                         return _ex0();
115                 if (x.is_equal(I))       // log(I) -> Pi*I/2
116                         return (Pi*I*_num1_2());
117                 if (x.is_equal(-I))      // log(-I) -> -Pi*I/2
118                         return (Pi*I*_num_1_2());
119                 // log(float)
120                 if (!x.info(info_flags::crational))
121                         return log(ex_to<numeric>(x));
122         }
123         // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
124         if (is_ex_the_function(x, exp)) {
125                 ex t = x.op(0);
126                 if (t.info(info_flags::numeric)) {
127                         numeric nt = ex_to<numeric>(t);
128                         if (nt.is_real())
129                                 return t;
130                 }
131         }
132         
133         return log(x).hold();
134 }
135
136 static ex log_deriv(const ex & x, unsigned deriv_param)
137 {
138         GINAC_ASSERT(deriv_param==0);
139         
140         // d/dx log(x) -> 1/x
141         return power(x, _ex_1());
142 }
143
144 static ex log_series(const ex &arg,
145                      const relational &rel,
146                      int order,
147                      unsigned options)
148 {
149         GINAC_ASSERT(is_exactly_a<symbol>(rel.lhs()));
150         ex arg_pt;
151         bool must_expand_arg = false;
152         // maybe substitution of rel into arg fails because of a pole
153         try {
154                 arg_pt = arg.subs(rel);
155         } catch (pole_error) {
156                 must_expand_arg = true;
157         }
158         // or we are at the branch point anyways
159         if (arg_pt.is_zero())
160                 must_expand_arg = true;
161         
162         if (must_expand_arg) {
163                 // method:
164                 // This is the branch point: Series expand the argument first, then
165                 // trivially factorize it to isolate that part which has constant
166                 // leading coefficient in this fashion:
167                 //   x^n + x^(n+1) +...+ Order(x^(n+m))  ->  x^n * (1 + x +...+ Order(x^m)).
168                 // Return a plain n*log(x) for the x^n part and series expand the
169                 // other part.  Add them together and reexpand again in order to have
170                 // one unnested pseries object.  All this also works for negative n.
171                 pseries argser;          // series expansion of log's argument
172                 unsigned extra_ord = 0;  // extra expansion order
173                 do {
174                         // oops, the argument expanded to a pure Order(x^something)...
175                         argser = ex_to<pseries>(arg.series(rel, order+extra_ord, options));
176                         ++extra_ord;
177                 } while (!argser.is_terminating() && argser.nops()==1);
178
179                 const symbol &s = ex_to<symbol>(rel.lhs());
180                 const ex point = rel.rhs();
181                 const int n = argser.ldegree(s);
182                 epvector seq;
183                 // construct what we carelessly called the n*log(x) term above
184                 const ex coeff = argser.coeff(s, n);
185                 // expand the log, but only if coeff is real and > 0, since otherwise
186                 // it would make the branch cut run into the wrong direction
187                 if (coeff.info(info_flags::positive))
188                         seq.push_back(expair(n*log(s-point)+log(coeff), _ex0()));
189                 else
190                         seq.push_back(expair(log(coeff*pow(s-point, n)), _ex0()));
191
192                 if (!argser.is_terminating() || argser.nops()!=1) {
193                         // in this case n more (or less) terms are needed
194                         // (sadly, to generate them, we have to start from the beginning)
195                         const ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
196                         return pseries(rel, seq).add_series(ex_to<pseries>(log(newarg).series(rel, order, options)));
197                 } else  // it was a monomial
198                         return pseries(rel, seq);
199         }
200         if (!(options & series_options::suppress_branchcut) &&
201              arg_pt.info(info_flags::negative)) {
202                 // method:
203                 // This is the branch cut: assemble the primitive series manually and
204                 // then add the corresponding complex step function.
205                 const symbol &s = ex_to<symbol>(rel.lhs());
206                 const ex point = rel.rhs();
207                 const symbol foo;
208                 const ex replarg = series(log(arg), s==foo, order).subs(foo==point);
209                 epvector seq;
210                 seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
211                 seq.push_back(expair(Order(_ex1()), order));
212                 return series(replarg - I*Pi + pseries(rel, seq), rel, order);
213         }
214         throw do_taylor();  // caught by function::series()
215 }
216
217 REGISTER_FUNCTION(log, eval_func(log_eval).
218                        evalf_func(log_evalf).
219                        derivative_func(log_deriv).
220                        series_func(log_series).
221                        latex_name("\\ln"));
222
223 //////////
224 // sine (trigonometric function)
225 //////////
226
227 static ex sin_evalf(const ex & x)
228 {
229         if (is_exactly_a<numeric>(x))
230                 return sin(ex_to<numeric>(x));
231         
232         return sin(x).hold();
233 }
234
235 static ex sin_eval(const ex & x)
236 {
237         // sin(n/d*Pi) -> { all known non-nested radicals }
238         const ex SixtyExOverPi = _ex60()*x/Pi;
239         ex sign = _ex1();
240         if (SixtyExOverPi.info(info_flags::integer)) {
241                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
242                 if (z>=_num60()) {
243                         // wrap to interval [0, Pi)
244                         z -= _num60();
245                         sign = _ex_1();
246                 }
247                 if (z>_num30()) {
248                         // wrap to interval [0, Pi/2)
249                         z = _num60()-z;
250                 }
251                 if (z.is_equal(_num0()))  // sin(0)       -> 0
252                         return _ex0();
253                 if (z.is_equal(_num5()))  // sin(Pi/12)   -> sqrt(6)/4*(1-sqrt(3)/3)
254                         return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
255                 if (z.is_equal(_num6()))  // sin(Pi/10)   -> sqrt(5)/4-1/4
256                         return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
257                 if (z.is_equal(_num10())) // sin(Pi/6)    -> 1/2
258                         return sign*_ex1_2();
259                 if (z.is_equal(_num15())) // sin(Pi/4)    -> sqrt(2)/2
260                         return sign*_ex1_2()*power(_ex2(),_ex1_2());
261                 if (z.is_equal(_num18())) // sin(3/10*Pi) -> sqrt(5)/4+1/4
262                         return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
263                 if (z.is_equal(_num20())) // sin(Pi/3)    -> sqrt(3)/2
264                         return sign*_ex1_2()*power(_ex3(),_ex1_2());
265                 if (z.is_equal(_num25())) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
266                         return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
267                 if (z.is_equal(_num30())) // sin(Pi/2)    -> 1
268                         return sign*_ex1();
269         }
270         
271         if (is_exactly_a<function>(x)) {
272                 ex t = x.op(0);
273                 // sin(asin(x)) -> x
274                 if (is_ex_the_function(x, asin))
275                         return t;
276                 // sin(acos(x)) -> sqrt(1-x^2)
277                 if (is_ex_the_function(x, acos))
278                         return power(_ex1()-power(t,_ex2()),_ex1_2());
279                 // sin(atan(x)) -> x*(1+x^2)^(-1/2)
280                 if (is_ex_the_function(x, atan))
281                         return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
282         }
283         
284         // sin(float) -> float
285         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
286                 return sin(ex_to<numeric>(x));
287         
288         return sin(x).hold();
289 }
290
291 static ex sin_deriv(const ex & x, unsigned deriv_param)
292 {
293         GINAC_ASSERT(deriv_param==0);
294         
295         // d/dx sin(x) -> cos(x)
296         return cos(x);
297 }
298
299 REGISTER_FUNCTION(sin, eval_func(sin_eval).
300                        evalf_func(sin_evalf).
301                        derivative_func(sin_deriv).
302                        latex_name("\\sin"));
303
304 //////////
305 // cosine (trigonometric function)
306 //////////
307
308 static ex cos_evalf(const ex & x)
309 {
310         if (is_exactly_a<numeric>(x))
311                 return cos(ex_to<numeric>(x));
312         
313         return cos(x).hold();
314 }
315
316 static ex cos_eval(const ex & x)
317 {
318         // cos(n/d*Pi) -> { all known non-nested radicals }
319         const ex SixtyExOverPi = _ex60()*x/Pi;
320         ex sign = _ex1();
321         if (SixtyExOverPi.info(info_flags::integer)) {
322                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
323                 if (z>=_num60()) {
324                         // wrap to interval [0, Pi)
325                         z = _num120()-z;
326                 }
327                 if (z>=_num30()) {
328                         // wrap to interval [0, Pi/2)
329                         z = _num60()-z;
330                         sign = _ex_1();
331                 }
332                 if (z.is_equal(_num0()))  // cos(0)       -> 1
333                         return sign*_ex1();
334                 if (z.is_equal(_num5()))  // cos(Pi/12)   -> sqrt(6)/4*(1+sqrt(3)/3)
335                         return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex1_3()*power(_ex3(),_ex1_2()));
336                 if (z.is_equal(_num10())) // cos(Pi/6)    -> sqrt(3)/2
337                         return sign*_ex1_2()*power(_ex3(),_ex1_2());
338                 if (z.is_equal(_num12())) // cos(Pi/5)    -> sqrt(5)/4+1/4
339                         return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex1_4());
340                 if (z.is_equal(_num15())) // cos(Pi/4)    -> sqrt(2)/2
341                         return sign*_ex1_2()*power(_ex2(),_ex1_2());
342                 if (z.is_equal(_num20())) // cos(Pi/3)    -> 1/2
343                         return sign*_ex1_2();
344                 if (z.is_equal(_num24())) // cos(2/5*Pi)  -> sqrt(5)/4-1/4x
345                         return sign*(_ex1_4()*power(_ex5(),_ex1_2())+_ex_1_4());
346                 if (z.is_equal(_num25())) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
347                         return sign*_ex1_4()*power(_ex6(),_ex1_2())*(_ex1()+_ex_1_3()*power(_ex3(),_ex1_2()));
348                 if (z.is_equal(_num30())) // cos(Pi/2)    -> 0
349                         return sign*_ex0();
350         }
351         
352         if (is_exactly_a<function>(x)) {
353                 ex t = x.op(0);
354                 // cos(acos(x)) -> x
355                 if (is_ex_the_function(x, acos))
356                         return t;
357                 // cos(asin(x)) -> (1-x^2)^(1/2)
358                 if (is_ex_the_function(x, asin))
359                         return power(_ex1()-power(t,_ex2()),_ex1_2());
360                 // cos(atan(x)) -> (1+x^2)^(-1/2)
361                 if (is_ex_the_function(x, atan))
362                         return power(_ex1()+power(t,_ex2()),_ex_1_2());
363         }
364         
365         // cos(float) -> float
366         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
367                 return cos(ex_to<numeric>(x));
368         
369         return cos(x).hold();
370 }
371
372 static ex cos_deriv(const ex & x, unsigned deriv_param)
373 {
374         GINAC_ASSERT(deriv_param==0);
375
376         // d/dx cos(x) -> -sin(x)
377         return _ex_1()*sin(x);
378 }
379
380 REGISTER_FUNCTION(cos, eval_func(cos_eval).
381                        evalf_func(cos_evalf).
382                        derivative_func(cos_deriv).
383                        latex_name("\\cos"));
384
385 //////////
386 // tangent (trigonometric function)
387 //////////
388
389 static ex tan_evalf(const ex & x)
390 {
391         if (is_exactly_a<numeric>(x))
392                 return tan(ex_to<numeric>(x));
393         
394         return tan(x).hold();
395 }
396
397 static ex tan_eval(const ex & x)
398 {
399         // tan(n/d*Pi) -> { all known non-nested radicals }
400         const ex SixtyExOverPi = _ex60()*x/Pi;
401         ex sign = _ex1();
402         if (SixtyExOverPi.info(info_flags::integer)) {
403                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num60());
404                 if (z>=_num60()) {
405                         // wrap to interval [0, Pi)
406                         z -= _num60();
407                 }
408                 if (z>=_num30()) {
409                         // wrap to interval [0, Pi/2)
410                         z = _num60()-z;
411                         sign = _ex_1();
412                 }
413                 if (z.is_equal(_num0()))  // tan(0)       -> 0
414                         return _ex0();
415                 if (z.is_equal(_num5()))  // tan(Pi/12)   -> 2-sqrt(3)
416                         return sign*(_ex2()-power(_ex3(),_ex1_2()));
417                 if (z.is_equal(_num10())) // tan(Pi/6)    -> sqrt(3)/3
418                         return sign*_ex1_3()*power(_ex3(),_ex1_2());
419                 if (z.is_equal(_num15())) // tan(Pi/4)    -> 1
420                         return sign*_ex1();
421                 if (z.is_equal(_num20())) // tan(Pi/3)    -> sqrt(3)
422                         return sign*power(_ex3(),_ex1_2());
423                 if (z.is_equal(_num25())) // tan(5/12*Pi) -> 2+sqrt(3)
424                         return sign*(power(_ex3(),_ex1_2())+_ex2());
425                 if (z.is_equal(_num30())) // tan(Pi/2)    -> infinity
426                         throw (pole_error("tan_eval(): simple pole",1));
427         }
428         
429         if (is_exactly_a<function>(x)) {
430                 ex t = x.op(0);
431                 // tan(atan(x)) -> x
432                 if (is_ex_the_function(x, atan))
433                         return t;
434                 // tan(asin(x)) -> x*(1+x^2)^(-1/2)
435                 if (is_ex_the_function(x, asin))
436                         return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
437                 // tan(acos(x)) -> (1-x^2)^(1/2)/x
438                 if (is_ex_the_function(x, acos))
439                         return power(t,_ex_1())*power(_ex1()-power(t,_ex2()),_ex1_2());
440         }
441         
442         // tan(float) -> float
443         if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
444                 return tan(ex_to<numeric>(x));
445         }
446         
447         return tan(x).hold();
448 }
449
450 static ex tan_deriv(const ex & x, unsigned deriv_param)
451 {
452         GINAC_ASSERT(deriv_param==0);
453         
454         // d/dx tan(x) -> 1+tan(x)^2;
455         return (_ex1()+power(tan(x),_ex2()));
456 }
457
458 static ex tan_series(const ex &x,
459                      const relational &rel,
460                      int order,
461                      unsigned options)
462 {
463         GINAC_ASSERT(is_exactly_a<symbol>(rel.lhs()));
464         // method:
465         // Taylor series where there is no pole falls back to tan_deriv.
466         // On a pole simply expand sin(x)/cos(x).
467         const ex x_pt = x.subs(rel);
468         if (!(2*x_pt/Pi).info(info_flags::odd))
469                 throw do_taylor();  // caught by function::series()
470         // if we got here we have to care for a simple pole
471         return (sin(x)/cos(x)).series(rel, order+2, options);
472 }
473
474 REGISTER_FUNCTION(tan, eval_func(tan_eval).
475                        evalf_func(tan_evalf).
476                        derivative_func(tan_deriv).
477                        series_func(tan_series).
478                        latex_name("\\tan"));
479
480 //////////
481 // inverse sine (arc sine)
482 //////////
483
484 static ex asin_evalf(const ex & x)
485 {
486         if (is_exactly_a<numeric>(x))
487                 return asin(ex_to<numeric>(x));
488         
489         return asin(x).hold();
490 }
491
492 static ex asin_eval(const ex & x)
493 {
494         if (x.info(info_flags::numeric)) {
495                 // asin(0) -> 0
496                 if (x.is_zero())
497                         return x;
498                 // asin(1/2) -> Pi/6
499                 if (x.is_equal(_ex1_2()))
500                         return numeric(1,6)*Pi;
501                 // asin(1) -> Pi/2
502                 if (x.is_equal(_ex1()))
503                         return _num1_2()*Pi;
504                 // asin(-1/2) -> -Pi/6
505                 if (x.is_equal(_ex_1_2()))
506                         return numeric(-1,6)*Pi;
507                 // asin(-1) -> -Pi/2
508                 if (x.is_equal(_ex_1()))
509                         return _num_1_2()*Pi;
510                 // asin(float) -> float
511                 if (!x.info(info_flags::crational))
512                         return asin(ex_to<numeric>(x));
513         }
514         
515         return asin(x).hold();
516 }
517
518 static ex asin_deriv(const ex & x, unsigned deriv_param)
519 {
520         GINAC_ASSERT(deriv_param==0);
521         
522         // d/dx asin(x) -> 1/sqrt(1-x^2)
523         return power(1-power(x,_ex2()),_ex_1_2());
524 }
525
526 REGISTER_FUNCTION(asin, eval_func(asin_eval).
527                         evalf_func(asin_evalf).
528                         derivative_func(asin_deriv).
529                         latex_name("\\arcsin"));
530
531 //////////
532 // inverse cosine (arc cosine)
533 //////////
534
535 static ex acos_evalf(const ex & x)
536 {
537         if (is_exactly_a<numeric>(x))
538                 return acos(ex_to<numeric>(x));
539         
540         return acos(x).hold();
541 }
542
543 static ex acos_eval(const ex & x)
544 {
545         if (x.info(info_flags::numeric)) {
546                 // acos(1) -> 0
547                 if (x.is_equal(_ex1()))
548                         return _ex0();
549                 // acos(1/2) -> Pi/3
550                 if (x.is_equal(_ex1_2()))
551                         return _ex1_3()*Pi;
552                 // acos(0) -> Pi/2
553                 if (x.is_zero())
554                         return _ex1_2()*Pi;
555                 // acos(-1/2) -> 2/3*Pi
556                 if (x.is_equal(_ex_1_2()))
557                         return numeric(2,3)*Pi;
558                 // acos(-1) -> Pi
559                 if (x.is_equal(_ex_1()))
560                         return Pi;
561                 // acos(float) -> float
562                 if (!x.info(info_flags::crational))
563                         return acos(ex_to<numeric>(x));
564         }
565         
566         return acos(x).hold();
567 }
568
569 static ex acos_deriv(const ex & x, unsigned deriv_param)
570 {
571         GINAC_ASSERT(deriv_param==0);
572         
573         // d/dx acos(x) -> -1/sqrt(1-x^2)
574         return _ex_1()*power(1-power(x,_ex2()),_ex_1_2());
575 }
576
577 REGISTER_FUNCTION(acos, eval_func(acos_eval).
578                         evalf_func(acos_evalf).
579                         derivative_func(acos_deriv).
580                         latex_name("\\arccos"));
581
582 //////////
583 // inverse tangent (arc tangent)
584 //////////
585
586 static ex atan_evalf(const ex & x)
587 {
588         if (is_exactly_a<numeric>(x))
589                 return atan(ex_to<numeric>(x));
590         
591         return atan(x).hold();
592 }
593
594 static ex atan_eval(const ex & x)
595 {
596         if (x.info(info_flags::numeric)) {
597                 // atan(0) -> 0
598                 if (x.is_zero())
599                         return _ex0();
600                 // atan(1) -> Pi/4
601                 if (x.is_equal(_ex1()))
602                         return _ex1_4()*Pi;
603                 // atan(-1) -> -Pi/4
604                 if (x.is_equal(_ex_1()))
605                         return _ex_1_4()*Pi;
606                 if (x.is_equal(I) || x.is_equal(-I))
607                         throw (pole_error("atan_eval(): logarithmic pole",0));
608                 // atan(float) -> float
609                 if (!x.info(info_flags::crational))
610                         return atan(ex_to<numeric>(x));
611         }
612         
613         return atan(x).hold();
614 }
615
616 static ex atan_deriv(const ex & x, unsigned deriv_param)
617 {
618         GINAC_ASSERT(deriv_param==0);
619
620         // d/dx atan(x) -> 1/(1+x^2)
621         return power(_ex1()+power(x,_ex2()), _ex_1());
622 }
623
624 static ex atan_series(const ex &arg,
625                       const relational &rel,
626                       int order,
627                       unsigned options)
628 {
629         GINAC_ASSERT(is_exactly_a<symbol>(rel.lhs()));
630         // method:
631         // Taylor series where there is no pole or cut falls back to atan_deriv.
632         // There are two branch cuts, one runnig from I up the imaginary axis and
633         // one running from -I down the imaginary axis.  The points I and -I are
634         // poles.
635         // On the branch cuts and the poles series expand
636         //     (log(1+I*x)-log(1-I*x))/(2*I)
637         // instead.
638         const ex arg_pt = arg.subs(rel);
639         if (!(I*arg_pt).info(info_flags::real))
640                 throw do_taylor();     // Re(x) != 0
641         if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1())
642                 throw do_taylor();     // Re(x) == 0, but abs(x)<1
643         // care for the poles, using the defining formula for atan()...
644         if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
645                 return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
646         if (!(options & series_options::suppress_branchcut)) {
647                 // method:
648                 // This is the branch cut: assemble the primitive series manually and
649                 // then add the corresponding complex step function.
650                 const symbol &s = ex_to<symbol>(rel.lhs());
651                 const ex point = rel.rhs();
652                 const symbol foo;
653                 const ex replarg = series(atan(arg), s==foo, order).subs(foo==point);
654                 ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2();
655                 if ((I*arg_pt)<_ex0())
656                         Order0correction += log((I*arg_pt+_ex_1())/(I*arg_pt+_ex1()))*I*_ex_1_2();
657                 else
658                         Order0correction += log((I*arg_pt+_ex1())/(I*arg_pt+_ex_1()))*I*_ex1_2();
659                 epvector seq;
660                 seq.push_back(expair(Order0correction, _ex0()));
661                 seq.push_back(expair(Order(_ex1()), order));
662                 return series(replarg - pseries(rel, seq), rel, order);
663         }
664         throw do_taylor();
665 }
666
667 REGISTER_FUNCTION(atan, eval_func(atan_eval).
668                         evalf_func(atan_evalf).
669                         derivative_func(atan_deriv).
670                         series_func(atan_series).
671                         latex_name("\\arctan"));
672
673 //////////
674 // inverse tangent (atan2(y,x))
675 //////////
676
677 static ex atan2_evalf(const ex &y, const ex &x)
678 {
679         if (is_exactly_a<numeric>(y) && is_exactly_a<numeric>(x))
680                 return atan2(ex_to<numeric>(y), ex_to<numeric>(x));
681         
682         return atan2(y, x).hold();
683 }
684
685 static ex atan2_eval(const ex & y, const ex & x)
686 {
687         if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
688                 x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
689                 return atan2_evalf(y,x);
690         }
691         
692         return atan2(y,x).hold();
693 }    
694
695 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
696 {
697         GINAC_ASSERT(deriv_param<2);
698         
699         if (deriv_param==0) {
700                 // d/dy atan(y,x)
701                 return x*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
702         }
703         // d/dx atan(y,x)
704         return -y*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
705 }
706
707 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
708                          evalf_func(atan2_evalf).
709                          derivative_func(atan2_deriv));
710
711 //////////
712 // hyperbolic sine (trigonometric function)
713 //////////
714
715 static ex sinh_evalf(const ex & x)
716 {
717         if (is_exactly_a<numeric>(x))
718                 return sinh(ex_to<numeric>(x));
719         
720         return sinh(x).hold();
721 }
722
723 static ex sinh_eval(const ex & x)
724 {
725         if (x.info(info_flags::numeric)) {
726                 if (x.is_zero())  // sinh(0) -> 0
727                         return _ex0();        
728                 if (!x.info(info_flags::crational))  // sinh(float) -> float
729                         return sinh(ex_to<numeric>(x));
730         }
731         
732         if ((x/Pi).info(info_flags::numeric) &&
733                 ex_to<numeric>(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
734                 return I*sin(x/I);
735         
736         if (is_exactly_a<function>(x)) {
737                 ex t = x.op(0);
738                 // sinh(asinh(x)) -> x
739                 if (is_ex_the_function(x, asinh))
740                         return t;
741                 // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2)
742                 if (is_ex_the_function(x, acosh))
743                         return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2());
744                 // sinh(atanh(x)) -> x*(1-x^2)^(-1/2)
745                 if (is_ex_the_function(x, atanh))
746                         return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
747         }
748         
749         return sinh(x).hold();
750 }
751
752 static ex sinh_deriv(const ex & x, unsigned deriv_param)
753 {
754         GINAC_ASSERT(deriv_param==0);
755         
756         // d/dx sinh(x) -> cosh(x)
757         return cosh(x);
758 }
759
760 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
761                         evalf_func(sinh_evalf).
762                         derivative_func(sinh_deriv).
763                         latex_name("\\sinh"));
764
765 //////////
766 // hyperbolic cosine (trigonometric function)
767 //////////
768
769 static ex cosh_evalf(const ex & x)
770 {
771         if (is_exactly_a<numeric>(x))
772                 return cosh(ex_to<numeric>(x));
773         
774         return cosh(x).hold();
775 }
776
777 static ex cosh_eval(const ex & x)
778 {
779         if (x.info(info_flags::numeric)) {
780                 if (x.is_zero())  // cosh(0) -> 1
781                         return _ex1();
782                 if (!x.info(info_flags::crational))  // cosh(float) -> float
783                         return cosh(ex_to<numeric>(x));
784         }
785         
786         if ((x/Pi).info(info_flags::numeric) &&
787                 ex_to<numeric>(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
788                 return cos(x/I);
789         
790         if (is_exactly_a<function>(x)) {
791                 ex t = x.op(0);
792                 // cosh(acosh(x)) -> x
793                 if (is_ex_the_function(x, acosh))
794                         return t;
795                 // cosh(asinh(x)) -> (1+x^2)^(1/2)
796                 if (is_ex_the_function(x, asinh))
797                         return power(_ex1()+power(t,_ex2()),_ex1_2());
798                 // cosh(atanh(x)) -> (1-x^2)^(-1/2)
799                 if (is_ex_the_function(x, atanh))
800                         return power(_ex1()-power(t,_ex2()),_ex_1_2());
801         }
802         
803         return cosh(x).hold();
804 }
805
806 static ex cosh_deriv(const ex & x, unsigned deriv_param)
807 {
808         GINAC_ASSERT(deriv_param==0);
809         
810         // d/dx cosh(x) -> sinh(x)
811         return sinh(x);
812 }
813
814 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
815                         evalf_func(cosh_evalf).
816                         derivative_func(cosh_deriv).
817                         latex_name("\\cosh"));
818
819 //////////
820 // hyperbolic tangent (trigonometric function)
821 //////////
822
823 static ex tanh_evalf(const ex & x)
824 {
825         if (is_exactly_a<numeric>(x))
826                 return tanh(ex_to<numeric>(x));
827         
828         return tanh(x).hold();
829 }
830
831 static ex tanh_eval(const ex & x)
832 {
833         if (x.info(info_flags::numeric)) {
834                 if (x.is_zero())  // tanh(0) -> 0
835                         return _ex0();
836                 if (!x.info(info_flags::crational))  // tanh(float) -> float
837                         return tanh(ex_to<numeric>(x));
838         }
839         
840         if ((x/Pi).info(info_flags::numeric) &&
841                 ex_to<numeric>(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
842                 return I*tan(x/I);
843         
844         if (is_exactly_a<function>(x)) {
845                 ex t = x.op(0);
846                 // tanh(atanh(x)) -> x
847                 if (is_ex_the_function(x, atanh))
848                         return t;
849                 // tanh(asinh(x)) -> x*(1+x^2)^(-1/2)
850                 if (is_ex_the_function(x, asinh))
851                         return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
852                 // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x
853                 if (is_ex_the_function(x, acosh))
854                         return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2())*power(t,_ex_1());
855         }
856         
857         return tanh(x).hold();
858 }
859
860 static ex tanh_deriv(const ex & x, unsigned deriv_param)
861 {
862         GINAC_ASSERT(deriv_param==0);
863         
864         // d/dx tanh(x) -> 1-tanh(x)^2
865         return _ex1()-power(tanh(x),_ex2());
866 }
867
868 static ex tanh_series(const ex &x,
869                       const relational &rel,
870                       int order,
871                       unsigned options)
872 {
873         GINAC_ASSERT(is_exactly_a<symbol>(rel.lhs()));
874         // method:
875         // Taylor series where there is no pole falls back to tanh_deriv.
876         // On a pole simply expand sinh(x)/cosh(x).
877         const ex x_pt = x.subs(rel);
878         if (!(2*I*x_pt/Pi).info(info_flags::odd))
879                 throw do_taylor();  // caught by function::series()
880         // if we got here we have to care for a simple pole
881         return (sinh(x)/cosh(x)).series(rel, order+2, options);
882 }
883
884 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
885                         evalf_func(tanh_evalf).
886                         derivative_func(tanh_deriv).
887                         series_func(tanh_series).
888                         latex_name("\\tanh"));
889
890 //////////
891 // inverse hyperbolic sine (trigonometric function)
892 //////////
893
894 static ex asinh_evalf(const ex & x)
895 {
896         if (is_exactly_a<numeric>(x))
897                 return asinh(ex_to<numeric>(x));
898         
899         return asinh(x).hold();
900 }
901
902 static ex asinh_eval(const ex & x)
903 {
904         if (x.info(info_flags::numeric)) {
905                 // asinh(0) -> 0
906                 if (x.is_zero())
907                         return _ex0();
908                 // asinh(float) -> float
909                 if (!x.info(info_flags::crational))
910                         return asinh(ex_to<numeric>(x));
911         }
912         
913         return asinh(x).hold();
914 }
915
916 static ex asinh_deriv(const ex & x, unsigned deriv_param)
917 {
918         GINAC_ASSERT(deriv_param==0);
919         
920         // d/dx asinh(x) -> 1/sqrt(1+x^2)
921         return power(_ex1()+power(x,_ex2()),_ex_1_2());
922 }
923
924 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
925                          evalf_func(asinh_evalf).
926                          derivative_func(asinh_deriv));
927
928 //////////
929 // inverse hyperbolic cosine (trigonometric function)
930 //////////
931
932 static ex acosh_evalf(const ex & x)
933 {
934         if (is_exactly_a<numeric>(x))
935                 return acosh(ex_to<numeric>(x));
936         
937         return acosh(x).hold();
938 }
939
940 static ex acosh_eval(const ex & x)
941 {
942         if (x.info(info_flags::numeric)) {
943                 // acosh(0) -> Pi*I/2
944                 if (x.is_zero())
945                         return Pi*I*numeric(1,2);
946                 // acosh(1) -> 0
947                 if (x.is_equal(_ex1()))
948                         return _ex0();
949                 // acosh(-1) -> Pi*I
950                 if (x.is_equal(_ex_1()))
951                         return Pi*I;
952                 // acosh(float) -> float
953                 if (!x.info(info_flags::crational))
954                         return acosh(ex_to<numeric>(x));
955         }
956         
957         return acosh(x).hold();
958 }
959
960 static ex acosh_deriv(const ex & x, unsigned deriv_param)
961 {
962         GINAC_ASSERT(deriv_param==0);
963         
964         // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
965         return power(x+_ex_1(),_ex_1_2())*power(x+_ex1(),_ex_1_2());
966 }
967
968 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
969                          evalf_func(acosh_evalf).
970                          derivative_func(acosh_deriv));
971
972 //////////
973 // inverse hyperbolic tangent (trigonometric function)
974 //////////
975
976 static ex atanh_evalf(const ex & x)
977 {
978         if (is_exactly_a<numeric>(x))
979                 return atanh(ex_to<numeric>(x));
980         
981         return atanh(x).hold();
982 }
983
984 static ex atanh_eval(const ex & x)
985 {
986         if (x.info(info_flags::numeric)) {
987                 // atanh(0) -> 0
988                 if (x.is_zero())
989                         return _ex0();
990                 // atanh({+|-}1) -> throw
991                 if (x.is_equal(_ex1()) || x.is_equal(_ex_1()))
992                         throw (pole_error("atanh_eval(): logarithmic pole",0));
993                 // atanh(float) -> float
994                 if (!x.info(info_flags::crational))
995                         return atanh(ex_to<numeric>(x));
996         }
997         
998         return atanh(x).hold();
999 }
1000
1001 static ex atanh_deriv(const ex & x, unsigned deriv_param)
1002 {
1003         GINAC_ASSERT(deriv_param==0);
1004         
1005         // d/dx atanh(x) -> 1/(1-x^2)
1006         return power(_ex1()-power(x,_ex2()),_ex_1());
1007 }
1008
1009 static ex atanh_series(const ex &arg,
1010                        const relational &rel,
1011                        int order,
1012                        unsigned options)
1013 {
1014         GINAC_ASSERT(is_exactly_a<symbol>(rel.lhs()));
1015         // method:
1016         // Taylor series where there is no pole or cut falls back to atanh_deriv.
1017         // There are two branch cuts, one runnig from 1 up the real axis and one
1018         // one running from -1 down the real axis.  The points 1 and -1 are poles
1019         // On the branch cuts and the poles series expand
1020         //     (log(1+x)-log(1-x))/2
1021         // instead.
1022         const ex arg_pt = arg.subs(rel);
1023         if (!(arg_pt).info(info_flags::real))
1024                 throw do_taylor();     // Im(x) != 0
1025         if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1())
1026                 throw do_taylor();     // Im(x) == 0, but abs(x)<1
1027         // care for the poles, using the defining formula for atanh()...
1028         if (arg_pt.is_equal(_ex1()) || arg_pt.is_equal(_ex_1()))
1029                 return ((log(_ex1()+arg)-log(_ex1()-arg))*_ex1_2()).series(rel, order, options);
1030         // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
1031         if (!(options & series_options::suppress_branchcut)) {
1032                 // method:
1033                 // This is the branch cut: assemble the primitive series manually and
1034                 // then add the corresponding complex step function.
1035                 const symbol &s = ex_to<symbol>(rel.lhs());
1036                 const ex point = rel.rhs();
1037                 const symbol foo;
1038                 const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point);
1039                 ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2();
1040                 if (arg_pt<_ex0())
1041                         Order0correction += log((arg_pt+_ex_1())/(arg_pt+_ex1()))*_ex1_2();
1042                 else
1043                         Order0correction += log((arg_pt+_ex1())/(arg_pt+_ex_1()))*_ex_1_2();
1044                 epvector seq;
1045                 seq.push_back(expair(Order0correction, _ex0()));
1046                 seq.push_back(expair(Order(_ex1()), order));
1047                 return series(replarg - pseries(rel, seq), rel, order);
1048         }
1049         throw do_taylor();
1050 }
1051
1052 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
1053                          evalf_func(atanh_evalf).
1054                          derivative_func(atanh_deriv).
1055                          series_func(atanh_series));
1056
1057
1058 } // namespace GiNaC