GiNaC  1.8.3
inifcns_trans.cpp
Go to the documentation of this file.
1 
6 /*
7  * GiNaC Copyright (C) 1999-2022 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22  */
23 
24 #include "inifcns.h"
25 #include "ex.h"
26 #include "constant.h"
27 #include "add.h"
28 #include "mul.h"
29 #include "numeric.h"
30 #include "power.h"
31 #include "operators.h"
32 #include "relational.h"
33 #include "symbol.h"
34 #include "pseries.h"
35 #include "utils.h"
36 
37 #include <stdexcept>
38 #include <vector>
39 
40 namespace GiNaC {
41 
43 // exponential function
45 
46 static ex exp_evalf(const ex & x)
47 {
48  if (is_exactly_a<numeric>(x))
49  return exp(ex_to<numeric>(x));
50 
51  return exp(x).hold();
52 }
53 
54 static ex exp_eval(const ex & x)
55 {
56  // exp(0) -> 1
57  if (x.is_zero()) {
58  return _ex1;
59  }
60 
61  // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
62  const ex TwoExOverPiI=(_ex2*x)/(Pi*I);
63  if (TwoExOverPiI.info(info_flags::integer)) {
64  const numeric z = mod(ex_to<numeric>(TwoExOverPiI),*_num4_p);
65  if (z.is_equal(*_num0_p))
66  return _ex1;
67  if (z.is_equal(*_num1_p))
68  return ex(I);
69  if (z.is_equal(*_num2_p))
70  return _ex_1;
71  if (z.is_equal(*_num3_p))
72  return ex(-I);
73  }
74 
75  // exp(log(x)) -> x
76  if (is_ex_the_function(x, log))
77  return x.op(0);
78 
79  // exp(float) -> float
81  return exp(ex_to<numeric>(x));
82 
83  return exp(x).hold();
84 }
85 
86 static ex exp_expand(const ex & arg, unsigned options)
87 {
88  ex exp_arg;
90  exp_arg = arg.expand(options);
91  else
92  exp_arg=arg;
93 
95  && is_exactly_a<add>(exp_arg)) {
96  exvector prodseq;
97  prodseq.reserve(exp_arg.nops());
98  for (const_iterator i = exp_arg.begin(); i != exp_arg.end(); ++i)
99  prodseq.push_back(exp(*i));
100 
101  return dynallocate<mul>(prodseq).setflag(status_flags::expanded);
102  }
103 
104  return exp(exp_arg).hold();
105 }
106 
107 static ex exp_deriv(const ex & x, unsigned deriv_param)
108 {
109  GINAC_ASSERT(deriv_param==0);
110 
111  // d/dx exp(x) -> exp(x)
112  return exp(x);
113 }
114 
115 static ex exp_real_part(const ex & x)
116 {
118 }
119 
120 static ex exp_imag_part(const ex & x)
121 {
123 }
124 
125 static ex exp_conjugate(const ex & x)
126 {
127  // conjugate(exp(x))==exp(conjugate(x))
128  return exp(x.conjugate());
129 }
130 
131 static ex exp_power(const ex & x, const ex & a)
132 {
133  /*
134  * The power law (e^x)^a=e^(x*a) is used in two cases:
135  * a) a is an integer and x may be complex;
136  * b) both x and a are reals.
137  * Negative a is excluded to keep automatic simplifications like exp(x)/exp(x)=1.
138  */
141  return exp(x*a);
142  else if (a.info(info_flags::negative)
144  return power(exp(-x*a), _ex_1).hold();
145 
146  return power(exp(x), a).hold();
147 }
148 
150  evalf_func(exp_evalf).
151  expand_func(exp_expand).
152  derivative_func(exp_deriv).
153  real_part_func(exp_real_part).
154  imag_part_func(exp_imag_part).
155  conjugate_func(exp_conjugate).
156  power_func(exp_power).
157  latex_name("\\exp"));
158 
160 // natural logarithm
162 
163 static ex log_evalf(const ex & x)
164 {
165  if (is_exactly_a<numeric>(x))
166  return log(ex_to<numeric>(x));
167 
168  return log(x).hold();
169 }
170 
171 static ex log_eval(const ex & x)
172 {
173  if (x.info(info_flags::numeric)) {
174  if (x.is_zero()) // log(0) -> infinity
175  throw(pole_error("log_eval(): log(0)",0));
177  return (log(-x)+I*Pi);
178  if (x.is_equal(_ex1)) // log(1) -> 0
179  return _ex0;
180  if (x.is_equal(I)) // log(I) -> Pi*I/2
181  return (Pi*I*_ex1_2);
182  if (x.is_equal(-I)) // log(-I) -> -Pi*I/2
183  return (Pi*I*_ex_1_2);
184 
185  // log(float) -> float
187  return log(ex_to<numeric>(x));
188  }
189 
190  // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
191  if (is_ex_the_function(x, exp)) {
192  const ex &t = x.op(0);
193  if (t.info(info_flags::real))
194  return t;
195  }
196 
197  return log(x).hold();
198 }
199 
200 static ex log_deriv(const ex & x, unsigned deriv_param)
201 {
202  GINAC_ASSERT(deriv_param==0);
203 
204  // d/dx log(x) -> 1/x
205  return power(x, _ex_1);
206 }
207 
208 static ex log_series(const ex &arg,
209  const relational &rel,
210  int order,
211  unsigned options)
212 {
213  GINAC_ASSERT(is_a<symbol>(rel.lhs()));
214  ex arg_pt;
215  bool must_expand_arg = false;
216  // maybe substitution of rel into arg fails because of a pole
217  try {
218  arg_pt = arg.subs(rel, subs_options::no_pattern);
219  } catch (pole_error &) {
220  must_expand_arg = true;
221  }
222  // or we are at the branch point anyways
223  if (arg_pt.is_zero())
224  must_expand_arg = true;
225 
226  if (arg.diff(ex_to<symbol>(rel.lhs())).is_zero()) {
227  throw do_taylor();
228  }
229 
230  if (must_expand_arg) {
231  // method:
232  // This is the branch point: Series expand the argument first, then
233  // trivially factorize it to isolate that part which has constant
234  // leading coefficient in this fashion:
235  // x^n + x^(n+1) +...+ Order(x^(n+m)) -> x^n * (1 + x +...+ Order(x^m)).
236  // Return a plain n*log(x) for the x^n part and series expand the
237  // other part. Add them together and reexpand again in order to have
238  // one unnested pseries object. All this also works for negative n.
239  pseries argser; // series expansion of log's argument
240  unsigned extra_ord = 0; // extra expansion order
241  do {
242  // oops, the argument expanded to a pure Order(x^something)...
243  argser = ex_to<pseries>(arg.series(rel, order+extra_ord, options));
244  ++extra_ord;
245  } while (!argser.is_terminating() && argser.nops()==1);
246 
247  const symbol &s = ex_to<symbol>(rel.lhs());
248  const ex &point = rel.rhs();
249  const int n = argser.ldegree(s);
250  epvector seq;
251  // construct what we carelessly called the n*log(x) term above
252  const ex coeff = argser.coeff(s, n);
253  // expand the log, but only if coeff is real and > 0, since otherwise
254  // it would make the branch cut run into the wrong direction
256  seq.push_back(expair(n*log(s-point)+log(coeff), _ex0));
257  else
258  seq.push_back(expair(log(coeff*pow(s-point, n)), _ex0));
259 
260  if (!argser.is_terminating() || argser.nops()!=1) {
261  // in this case n more (or less) terms are needed
262  // (sadly, to generate them, we have to start from the beginning)
263  if (n == 0 && coeff == 1) {
264  ex rest = pseries(rel, epvector{expair(-1, _ex0), expair(Order(_ex1), order)}).add_series(argser);
265  ex acc = dynallocate<pseries>(rel, epvector());
266  for (int i = order-1; i>0; --i) {
267  epvector cterm { expair(i%2 ? _ex1/i : _ex_1/i, _ex0) };
268  acc = pseries(rel, std::move(cterm)).add_series(ex_to<pseries>(acc));
269  acc = (ex_to<pseries>(rest)).mul_series(ex_to<pseries>(acc));
270  }
271  return acc;
272  }
273  const ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
274  return pseries(rel, std::move(seq)).add_series(ex_to<pseries>(log(newarg).series(rel, order, options)));
275  } else // it was a monomial
276  return pseries(rel, std::move(seq));
277  }
279  arg_pt.info(info_flags::negative)) {
280  // method:
281  // This is the branch cut: assemble the primitive series manually and
282  // then add the corresponding complex step function.
283  const symbol &s = ex_to<symbol>(rel.lhs());
284  const ex &point = rel.rhs();
285  const symbol foo;
286  const ex replarg = series(log(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
287  epvector seq;
288  if (order > 0) {
289  seq.reserve(2);
290  seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0));
291  }
292  seq.push_back(expair(Order(_ex1), order));
293  return series(replarg - I*Pi + pseries(rel, std::move(seq)), rel, order);
294  }
295  throw do_taylor(); // caught by function::series()
296 }
297 
298 static ex log_real_part(const ex & x)
299 {
301  return log(x).hold();
302  return log(abs(x));
303 }
304 
305 static ex log_imag_part(const ex & x)
306 {
308  return 0;
309  return atan2(GiNaC::imag_part(x), GiNaC::real_part(x));
310 }
311 
312 static ex log_expand(const ex & arg, unsigned options)
313 {
315  && is_exactly_a<mul>(arg) && !arg.info(info_flags::indefinite)) {
316  exvector sumseq;
317  exvector prodseq;
318  sumseq.reserve(arg.nops());
319  prodseq.reserve(arg.nops());
320  bool possign=true;
321 
322  // searching for positive/negative factors
323  for (const_iterator i = arg.begin(); i != arg.end(); ++i) {
324  ex e;
326  e=i->expand(options);
327  else
328  e=*i;
329  if (e.info(info_flags::positive))
330  sumseq.push_back(log(e));
331  else if (e.info(info_flags::negative)) {
332  sumseq.push_back(log(-e));
333  possign = !possign;
334  } else
335  prodseq.push_back(e);
336  }
337 
338  if (sumseq.size() > 0) {
339  ex newarg;
341  newarg=((possign?_ex1:_ex_1)*mul(prodseq)).expand(options);
342  else {
343  newarg=(possign?_ex1:_ex_1)*mul(prodseq);
344  ex_to<basic>(newarg).setflag(status_flags::purely_indefinite);
345  }
346  return add(sumseq)+log(newarg);
347  } else {
349  ex_to<basic>(arg).setflag(status_flags::purely_indefinite);
350  }
351  }
352 
354  return log(arg.expand(options)).hold();
355  else
356  return log(arg).hold();
357 }
358 
359 static ex log_conjugate(const ex & x)
360 {
361  // conjugate(log(x))==log(conjugate(x)) unless on the branch cut which
362  // runs along the negative real axis.
363  if (x.info(info_flags::positive)) {
364  return log(x);
365  }
366  if (is_exactly_a<numeric>(x) &&
367  !x.imag_part().is_zero()) {
368  return log(x.conjugate());
369  }
370  return conjugate_function(log(x)).hold();
371 }
372 
374  evalf_func(log_evalf).
375  expand_func(log_expand).
376  derivative_func(log_deriv).
377  series_func(log_series).
378  real_part_func(log_real_part).
379  imag_part_func(log_imag_part).
380  conjugate_func(log_conjugate).
381  latex_name("\\ln"));
382 
384 // sine (trigonometric function)
386 
387 static ex sin_evalf(const ex & x)
388 {
389  if (is_exactly_a<numeric>(x))
390  return sin(ex_to<numeric>(x));
391 
392  return sin(x).hold();
393 }
394 
395 static ex sin_eval(const ex & x)
396 {
397  // sin(n/d*Pi) -> { all known non-nested radicals }
398  const ex SixtyExOverPi = _ex60*x/Pi;
399  ex sign = _ex1;
400  if (SixtyExOverPi.info(info_flags::integer)) {
401  numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p);
402  if (z>=*_num60_p) {
403  // wrap to interval [0, Pi)
404  z -= *_num60_p;
405  sign = _ex_1;
406  }
407  if (z>*_num30_p) {
408  // wrap to interval [0, Pi/2)
409  z = *_num60_p-z;
410  }
411  if (z.is_equal(*_num0_p)) // sin(0) -> 0
412  return _ex0;
413  if (z.is_equal(*_num5_p)) // sin(Pi/12) -> sqrt(6)/4*(1-sqrt(3)/3)
414  return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
415  if (z.is_equal(*_num6_p)) // sin(Pi/10) -> sqrt(5)/4-1/4
416  return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
417  if (z.is_equal(*_num10_p)) // sin(Pi/6) -> 1/2
418  return sign*_ex1_2;
419  if (z.is_equal(*_num15_p)) // sin(Pi/4) -> sqrt(2)/2
420  return sign*_ex1_2*sqrt(_ex2);
421  if (z.is_equal(*_num18_p)) // sin(3/10*Pi) -> sqrt(5)/4+1/4
422  return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
423  if (z.is_equal(*_num20_p)) // sin(Pi/3) -> sqrt(3)/2
424  return sign*_ex1_2*sqrt(_ex3);
425  if (z.is_equal(*_num25_p)) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
426  return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
427  if (z.is_equal(*_num30_p)) // sin(Pi/2) -> 1
428  return sign;
429  }
430 
431  if (is_exactly_a<function>(x)) {
432  const ex &t = x.op(0);
433 
434  // sin(asin(x)) -> x
435  if (is_ex_the_function(x, asin))
436  return t;
437 
438  // sin(acos(x)) -> sqrt(1-x^2)
439  if (is_ex_the_function(x, acos))
440  return sqrt(_ex1-power(t,_ex2));
441 
442  // sin(atan(x)) -> x/sqrt(1+x^2)
443  if (is_ex_the_function(x, atan))
444  return t*power(_ex1+power(t,_ex2),_ex_1_2);
445  }
446 
447  // sin(float) -> float
449  return sin(ex_to<numeric>(x));
450 
451  // sin() is odd
453  return -sin(-x);
454 
455  return sin(x).hold();
456 }
457 
458 static ex sin_deriv(const ex & x, unsigned deriv_param)
459 {
460  GINAC_ASSERT(deriv_param==0);
461 
462  // d/dx sin(x) -> cos(x)
463  return cos(x);
464 }
465 
466 static ex sin_real_part(const ex & x)
467 {
469 }
470 
471 static ex sin_imag_part(const ex & x)
472 {
474 }
475 
476 static ex sin_conjugate(const ex & x)
477 {
478  // conjugate(sin(x))==sin(conjugate(x))
479  return sin(x.conjugate());
480 }
481 
483  evalf_func(sin_evalf).
484  derivative_func(sin_deriv).
485  real_part_func(sin_real_part).
486  imag_part_func(sin_imag_part).
487  conjugate_func(sin_conjugate).
488  latex_name("\\sin"));
489 
491 // cosine (trigonometric function)
493 
494 static ex cos_evalf(const ex & x)
495 {
496  if (is_exactly_a<numeric>(x))
497  return cos(ex_to<numeric>(x));
498 
499  return cos(x).hold();
500 }
501 
502 static ex cos_eval(const ex & x)
503 {
504  // cos(n/d*Pi) -> { all known non-nested radicals }
505  const ex SixtyExOverPi = _ex60*x/Pi;
506  ex sign = _ex1;
507  if (SixtyExOverPi.info(info_flags::integer)) {
508  numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p);
509  if (z>=*_num60_p) {
510  // wrap to interval [0, Pi)
511  z = *_num120_p-z;
512  }
513  if (z>=*_num30_p) {
514  // wrap to interval [0, Pi/2)
515  z = *_num60_p-z;
516  sign = _ex_1;
517  }
518  if (z.is_equal(*_num0_p)) // cos(0) -> 1
519  return sign;
520  if (z.is_equal(*_num5_p)) // cos(Pi/12) -> sqrt(6)/4*(1+sqrt(3)/3)
521  return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
522  if (z.is_equal(*_num10_p)) // cos(Pi/6) -> sqrt(3)/2
523  return sign*_ex1_2*sqrt(_ex3);
524  if (z.is_equal(*_num12_p)) // cos(Pi/5) -> sqrt(5)/4+1/4
525  return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
526  if (z.is_equal(*_num15_p)) // cos(Pi/4) -> sqrt(2)/2
527  return sign*_ex1_2*sqrt(_ex2);
528  if (z.is_equal(*_num20_p)) // cos(Pi/3) -> 1/2
529  return sign*_ex1_2;
530  if (z.is_equal(*_num24_p)) // cos(2/5*Pi) -> sqrt(5)/4-1/4x
531  return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
532  if (z.is_equal(*_num25_p)) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
533  return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
534  if (z.is_equal(*_num30_p)) // cos(Pi/2) -> 0
535  return _ex0;
536  }
537 
538  if (is_exactly_a<function>(x)) {
539  const ex &t = x.op(0);
540 
541  // cos(acos(x)) -> x
542  if (is_ex_the_function(x, acos))
543  return t;
544 
545  // cos(asin(x)) -> sqrt(1-x^2)
546  if (is_ex_the_function(x, asin))
547  return sqrt(_ex1-power(t,_ex2));
548 
549  // cos(atan(x)) -> 1/sqrt(1+x^2)
550  if (is_ex_the_function(x, atan))
551  return power(_ex1+power(t,_ex2),_ex_1_2);
552  }
553 
554  // cos(float) -> float
556  return cos(ex_to<numeric>(x));
557 
558  // cos() is even
560  return cos(-x);
561 
562  return cos(x).hold();
563 }
564 
565 static ex cos_deriv(const ex & x, unsigned deriv_param)
566 {
567  GINAC_ASSERT(deriv_param==0);
568 
569  // d/dx cos(x) -> -sin(x)
570  return -sin(x);
571 }
572 
573 static ex cos_real_part(const ex & x)
574 {
576 }
577 
578 static ex cos_imag_part(const ex & x)
579 {
581 }
582 
583 static ex cos_conjugate(const ex & x)
584 {
585  // conjugate(cos(x))==cos(conjugate(x))
586  return cos(x.conjugate());
587 }
588 
590  evalf_func(cos_evalf).
591  derivative_func(cos_deriv).
592  real_part_func(cos_real_part).
593  imag_part_func(cos_imag_part).
594  conjugate_func(cos_conjugate).
595  latex_name("\\cos"));
596 
598 // tangent (trigonometric function)
600 
601 static ex tan_evalf(const ex & x)
602 {
603  if (is_exactly_a<numeric>(x))
604  return tan(ex_to<numeric>(x));
605 
606  return tan(x).hold();
607 }
608 
609 static ex tan_eval(const ex & x)
610 {
611  // tan(n/d*Pi) -> { all known non-nested radicals }
612  const ex SixtyExOverPi = _ex60*x/Pi;
613  ex sign = _ex1;
614  if (SixtyExOverPi.info(info_flags::integer)) {
615  numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num60_p);
616  if (z>=*_num60_p) {
617  // wrap to interval [0, Pi)
618  z -= *_num60_p;
619  }
620  if (z>=*_num30_p) {
621  // wrap to interval [0, Pi/2)
622  z = *_num60_p-z;
623  sign = _ex_1;
624  }
625  if (z.is_equal(*_num0_p)) // tan(0) -> 0
626  return _ex0;
627  if (z.is_equal(*_num5_p)) // tan(Pi/12) -> 2-sqrt(3)
628  return sign*(_ex2-sqrt(_ex3));
629  if (z.is_equal(*_num10_p)) // tan(Pi/6) -> sqrt(3)/3
630  return sign*_ex1_3*sqrt(_ex3);
631  if (z.is_equal(*_num15_p)) // tan(Pi/4) -> 1
632  return sign;
633  if (z.is_equal(*_num20_p)) // tan(Pi/3) -> sqrt(3)
634  return sign*sqrt(_ex3);
635  if (z.is_equal(*_num25_p)) // tan(5/12*Pi) -> 2+sqrt(3)
636  return sign*(sqrt(_ex3)+_ex2);
637  if (z.is_equal(*_num30_p)) // tan(Pi/2) -> infinity
638  throw (pole_error("tan_eval(): simple pole",1));
639  }
640 
641  if (is_exactly_a<function>(x)) {
642  const ex &t = x.op(0);
643 
644  // tan(atan(x)) -> x
645  if (is_ex_the_function(x, atan))
646  return t;
647 
648  // tan(asin(x)) -> x/sqrt(1+x^2)
649  if (is_ex_the_function(x, asin))
650  return t*power(_ex1-power(t,_ex2),_ex_1_2);
651 
652  // tan(acos(x)) -> sqrt(1-x^2)/x
653  if (is_ex_the_function(x, acos))
654  return power(t,_ex_1)*sqrt(_ex1-power(t,_ex2));
655  }
656 
657  // tan(float) -> float
659  return tan(ex_to<numeric>(x));
660  }
661 
662  // tan() is odd
664  return -tan(-x);
665 
666  return tan(x).hold();
667 }
668 
669 static ex tan_deriv(const ex & x, unsigned deriv_param)
670 {
671  GINAC_ASSERT(deriv_param==0);
672 
673  // d/dx tan(x) -> 1+tan(x)^2;
674  return (_ex1+power(tan(x),_ex2));
675 }
676 
677 static ex tan_real_part(const ex & x)
678 {
679  ex a = GiNaC::real_part(x);
680  ex b = GiNaC::imag_part(x);
681  return tan(a)/(1+power(tan(a),2)*power(tan(b),2));
682 }
683 
684 static ex tan_imag_part(const ex & x)
685 {
686  ex a = GiNaC::real_part(x);
687  ex b = GiNaC::imag_part(x);
688  return tanh(b)/(1+power(tan(a),2)*power(tan(b),2));
689 }
690 
691 static ex tan_series(const ex &x,
692  const relational &rel,
693  int order,
694  unsigned options)
695 {
696  GINAC_ASSERT(is_a<symbol>(rel.lhs()));
697  // method:
698  // Taylor series where there is no pole falls back to tan_deriv.
699  // On a pole simply expand sin(x)/cos(x).
700  const ex x_pt = x.subs(rel, subs_options::no_pattern);
701  if (!(2*x_pt/Pi).info(info_flags::odd))
702  throw do_taylor(); // caught by function::series()
703  // if we got here we have to care for a simple pole
704  return (sin(x)/cos(x)).series(rel, order, options);
705 }
706 
707 static ex tan_conjugate(const ex & x)
708 {
709  // conjugate(tan(x))==tan(conjugate(x))
710  return tan(x.conjugate());
711 }
712 
714  evalf_func(tan_evalf).
715  derivative_func(tan_deriv).
716  series_func(tan_series).
717  real_part_func(tan_real_part).
718  imag_part_func(tan_imag_part).
719  conjugate_func(tan_conjugate).
720  latex_name("\\tan"));
721 
723 // inverse sine (arc sine)
725 
726 static ex asin_evalf(const ex & x)
727 {
728  if (is_exactly_a<numeric>(x))
729  return asin(ex_to<numeric>(x));
730 
731  return asin(x).hold();
732 }
733 
734 static ex asin_eval(const ex & x)
735 {
736  if (x.info(info_flags::numeric)) {
737 
738  // asin(0) -> 0
739  if (x.is_zero())
740  return x;
741 
742  // asin(1/2) -> Pi/6
743  if (x.is_equal(_ex1_2))
744  return numeric(1,6)*Pi;
745 
746  // asin(1) -> Pi/2
747  if (x.is_equal(_ex1))
748  return _ex1_2*Pi;
749 
750  // asin(-1/2) -> -Pi/6
751  if (x.is_equal(_ex_1_2))
752  return numeric(-1,6)*Pi;
753 
754  // asin(-1) -> -Pi/2
755  if (x.is_equal(_ex_1))
756  return _ex_1_2*Pi;
757 
758  // asin(float) -> float
760  return asin(ex_to<numeric>(x));
761 
762  // asin() is odd
764  return -asin(-x);
765  }
766 
767  return asin(x).hold();
768 }
769 
770 static ex asin_deriv(const ex & x, unsigned deriv_param)
771 {
772  GINAC_ASSERT(deriv_param==0);
773 
774  // d/dx asin(x) -> 1/sqrt(1-x^2)
775  return power(1-power(x,_ex2),_ex_1_2);
776 }
777 
778 static ex asin_conjugate(const ex & x)
779 {
780  // conjugate(asin(x))==asin(conjugate(x)) unless on the branch cuts which
781  // run along the real axis outside the interval [-1, +1].
782  if (is_exactly_a<numeric>(x) &&
783  (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) {
784  return asin(x.conjugate());
785  }
786  return conjugate_function(asin(x)).hold();
787 }
788 
790  evalf_func(asin_evalf).
791  derivative_func(asin_deriv).
792  conjugate_func(asin_conjugate).
793  latex_name("\\arcsin"));
794 
796 // inverse cosine (arc cosine)
798 
799 static ex acos_evalf(const ex & x)
800 {
801  if (is_exactly_a<numeric>(x))
802  return acos(ex_to<numeric>(x));
803 
804  return acos(x).hold();
805 }
806 
807 static ex acos_eval(const ex & x)
808 {
809  if (x.info(info_flags::numeric)) {
810 
811  // acos(1) -> 0
812  if (x.is_equal(_ex1))
813  return _ex0;
814 
815  // acos(1/2) -> Pi/3
816  if (x.is_equal(_ex1_2))
817  return _ex1_3*Pi;
818 
819  // acos(0) -> Pi/2
820  if (x.is_zero())
821  return _ex1_2*Pi;
822 
823  // acos(-1/2) -> 2/3*Pi
824  if (x.is_equal(_ex_1_2))
825  return numeric(2,3)*Pi;
826 
827  // acos(-1) -> Pi
828  if (x.is_equal(_ex_1))
829  return Pi;
830 
831  // acos(float) -> float
833  return acos(ex_to<numeric>(x));
834 
835  // acos(-x) -> Pi-acos(x)
837  return Pi-acos(-x);
838  }
839 
840  return acos(x).hold();
841 }
842 
843 static ex acos_deriv(const ex & x, unsigned deriv_param)
844 {
845  GINAC_ASSERT(deriv_param==0);
846 
847  // d/dx acos(x) -> -1/sqrt(1-x^2)
848  return -power(1-power(x,_ex2),_ex_1_2);
849 }
850 
851 static ex acos_conjugate(const ex & x)
852 {
853  // conjugate(acos(x))==acos(conjugate(x)) unless on the branch cuts which
854  // run along the real axis outside the interval [-1, +1].
855  if (is_exactly_a<numeric>(x) &&
856  (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) {
857  return acos(x.conjugate());
858  }
859  return conjugate_function(acos(x)).hold();
860 }
861 
863  evalf_func(acos_evalf).
864  derivative_func(acos_deriv).
865  conjugate_func(acos_conjugate).
866  latex_name("\\arccos"));
867 
869 // inverse tangent (arc tangent)
871 
872 static ex atan_evalf(const ex & x)
873 {
874  if (is_exactly_a<numeric>(x))
875  return atan(ex_to<numeric>(x));
876 
877  return atan(x).hold();
878 }
879 
880 static ex atan_eval(const ex & x)
881 {
882  if (x.info(info_flags::numeric)) {
883 
884  // atan(0) -> 0
885  if (x.is_zero())
886  return _ex0;
887 
888  // atan(1) -> Pi/4
889  if (x.is_equal(_ex1))
890  return _ex1_4*Pi;
891 
892  // atan(-1) -> -Pi/4
893  if (x.is_equal(_ex_1))
894  return _ex_1_4*Pi;
895 
896  if (x.is_equal(I) || x.is_equal(-I))
897  throw (pole_error("atan_eval(): logarithmic pole",0));
898 
899  // atan(float) -> float
901  return atan(ex_to<numeric>(x));
902 
903  // atan() is odd
905  return -atan(-x);
906  }
907 
908  return atan(x).hold();
909 }
910 
911 static ex atan_deriv(const ex & x, unsigned deriv_param)
912 {
913  GINAC_ASSERT(deriv_param==0);
914 
915  // d/dx atan(x) -> 1/(1+x^2)
916  return power(_ex1+power(x,_ex2), _ex_1);
917 }
918 
919 static ex atan_series(const ex &arg,
920  const relational &rel,
921  int order,
922  unsigned options)
923 {
924  GINAC_ASSERT(is_a<symbol>(rel.lhs()));
925  // method:
926  // Taylor series where there is no pole or cut falls back to atan_deriv.
927  // There are two branch cuts, one runnig from I up the imaginary axis and
928  // one running from -I down the imaginary axis. The points I and -I are
929  // poles.
930  // On the branch cuts and the poles series expand
931  // (log(1+I*x)-log(1-I*x))/(2*I)
932  // instead.
933  const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
934  if (!(I*arg_pt).info(info_flags::real))
935  throw do_taylor(); // Re(x) != 0
936  if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1)
937  throw do_taylor(); // Re(x) == 0, but abs(x)<1
938  // care for the poles, using the defining formula for atan()...
939  if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
940  return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
942  // method:
943  // This is the branch cut: assemble the primitive series manually and
944  // then add the corresponding complex step function.
945  const symbol &s = ex_to<symbol>(rel.lhs());
946  const ex &point = rel.rhs();
947  const symbol foo;
948  const ex replarg = series(atan(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
949  ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2;
950  if ((I*arg_pt)<_ex0)
951  Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2;
952  else
953  Order0correction += log((I*arg_pt+_ex1)/(I*arg_pt+_ex_1))*I*_ex1_2;
954  epvector seq;
955  if (order > 0) {
956  seq.reserve(2);
957  seq.push_back(expair(Order0correction, _ex0));
958  }
959  seq.push_back(expair(Order(_ex1), order));
960  return series(replarg - pseries(rel, std::move(seq)), rel, order);
961  }
962  throw do_taylor();
963 }
964 
965 static ex atan_conjugate(const ex & x)
966 {
967  // conjugate(atan(x))==atan(conjugate(x)) unless on the branch cuts which
968  // run along the imaginary axis outside the interval [-I, +I].
969  if (x.info(info_flags::real))
970  return atan(x);
971  if (is_exactly_a<numeric>(x)) {
972  const numeric x_re = ex_to<numeric>(x.real_part());
973  const numeric x_im = ex_to<numeric>(x.imag_part());
974  if (!x_re.is_zero() ||
975  (x_im > *_num_1_p && x_im < *_num1_p))
976  return atan(x.conjugate());
977  }
978  return conjugate_function(atan(x)).hold();
979 }
980 
982  evalf_func(atan_evalf).
983  derivative_func(atan_deriv).
984  series_func(atan_series).
985  conjugate_func(atan_conjugate).
986  latex_name("\\arctan"));
987 
989 // inverse tangent (atan2(y,x))
991 
992 static ex atan2_evalf(const ex &y, const ex &x)
993 {
994  if (is_exactly_a<numeric>(y) && is_exactly_a<numeric>(x))
995  return atan(ex_to<numeric>(y), ex_to<numeric>(x));
996 
997  return atan2(y, x).hold();
998 }
999 
1000 static ex atan2_eval(const ex & y, const ex & x)
1001 {
1002  if (y.is_zero()) {
1003 
1004  // atan2(0, 0) -> 0
1005  if (x.is_zero())
1006  return _ex0;
1007 
1008  // atan2(0, x), x real and positive -> 0
1010  return _ex0;
1011 
1012  // atan2(0, x), x real and negative -> Pi
1014  return Pi;
1015  }
1016 
1017  if (x.is_zero()) {
1018 
1019  // atan2(y, 0), y real and positive -> Pi/2
1020  if (y.info(info_flags::positive))
1021  return _ex1_2*Pi;
1022 
1023  // atan2(y, 0), y real and negative -> -Pi/2
1024  if (y.info(info_flags::negative))
1025  return _ex_1_2*Pi;
1026  }
1027 
1028  if (y.is_equal(x)) {
1029 
1030  // atan2(y, y), y real and positive -> Pi/4
1031  if (y.info(info_flags::positive))
1032  return _ex1_4*Pi;
1033 
1034  // atan2(y, y), y real and negative -> -3/4*Pi
1035  if (y.info(info_flags::negative))
1036  return numeric(-3, 4)*Pi;
1037  }
1038 
1039  if (y.is_equal(-x)) {
1040 
1041  // atan2(y, -y), y real and positive -> 3*Pi/4
1042  if (y.info(info_flags::positive))
1043  return numeric(3, 4)*Pi;
1044 
1045  // atan2(y, -y), y real and negative -> -Pi/4
1046  if (y.info(info_flags::negative))
1047  return _ex_1_4*Pi;
1048  }
1049 
1050  // atan2(float, float) -> float
1051  if (is_a<numeric>(y) && !y.info(info_flags::crational) &&
1052  is_a<numeric>(x) && !x.info(info_flags::crational))
1053  return atan(ex_to<numeric>(y), ex_to<numeric>(x));
1054 
1055  // atan2(real, real) -> atan(y/x) +/- Pi
1058  return atan(y/x);
1059 
1060  if (x.info(info_flags::negative)) {
1061  if (y.info(info_flags::positive))
1062  return atan(y/x)+Pi;
1063  if (y.info(info_flags::negative))
1064  return atan(y/x)-Pi;
1065  }
1066  }
1067 
1068  return atan2(y, x).hold();
1069 }
1070 
1071 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
1072 {
1073  GINAC_ASSERT(deriv_param<2);
1074 
1075  if (deriv_param==0) {
1076  // d/dy atan2(y,x)
1077  return x*power(power(x,_ex2)+power(y,_ex2),_ex_1);
1078  }
1079  // d/dx atan2(y,x)
1080  return -y*power(power(x,_ex2)+power(y,_ex2),_ex_1);
1081 }
1082 
1083 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
1084  evalf_func(atan2_evalf).
1085  derivative_func(atan2_deriv));
1086 
1088 // hyperbolic sine (trigonometric function)
1090 
1091 static ex sinh_evalf(const ex & x)
1092 {
1093  if (is_exactly_a<numeric>(x))
1094  return sinh(ex_to<numeric>(x));
1095 
1096  return sinh(x).hold();
1097 }
1098 
1099 static ex sinh_eval(const ex & x)
1100 {
1101  if (x.info(info_flags::numeric)) {
1102 
1103  // sinh(0) -> 0
1104  if (x.is_zero())
1105  return _ex0;
1106 
1107  // sinh(float) -> float
1109  return sinh(ex_to<numeric>(x));
1110 
1111  // sinh() is odd
1113  return -sinh(-x);
1114  }
1115 
1116  if ((x/Pi).info(info_flags::numeric) &&
1117  ex_to<numeric>(x/Pi).real().is_zero()) // sinh(I*x) -> I*sin(x)
1118  return I*sin(x/I);
1119 
1120  if (is_exactly_a<function>(x)) {
1121  const ex &t = x.op(0);
1122 
1123  // sinh(asinh(x)) -> x
1124  if (is_ex_the_function(x, asinh))
1125  return t;
1126 
1127  // sinh(acosh(x)) -> sqrt(x-1) * sqrt(x+1)
1128  if (is_ex_the_function(x, acosh))
1129  return sqrt(t-_ex1)*sqrt(t+_ex1);
1130 
1131  // sinh(atanh(x)) -> x/sqrt(1-x^2)
1132  if (is_ex_the_function(x, atanh))
1133  return t*power(_ex1-power(t,_ex2),_ex_1_2);
1134  }
1135 
1136  return sinh(x).hold();
1137 }
1138 
1139 static ex sinh_deriv(const ex & x, unsigned deriv_param)
1140 {
1141  GINAC_ASSERT(deriv_param==0);
1142 
1143  // d/dx sinh(x) -> cosh(x)
1144  return cosh(x);
1145 }
1146 
1147 static ex sinh_real_part(const ex & x)
1148 {
1150 }
1151 
1152 static ex sinh_imag_part(const ex & x)
1153 {
1155 }
1156 
1157 static ex sinh_conjugate(const ex & x)
1158 {
1159  // conjugate(sinh(x))==sinh(conjugate(x))
1160  return sinh(x.conjugate());
1161 }
1162 
1164  evalf_func(sinh_evalf).
1165  derivative_func(sinh_deriv).
1166  real_part_func(sinh_real_part).
1167  imag_part_func(sinh_imag_part).
1168  conjugate_func(sinh_conjugate).
1169  latex_name("\\sinh"));
1170 
1172 // hyperbolic cosine (trigonometric function)
1174 
1175 static ex cosh_evalf(const ex & x)
1176 {
1177  if (is_exactly_a<numeric>(x))
1178  return cosh(ex_to<numeric>(x));
1179 
1180  return cosh(x).hold();
1181 }
1182 
1183 static ex cosh_eval(const ex & x)
1184 {
1185  if (x.info(info_flags::numeric)) {
1186 
1187  // cosh(0) -> 1
1188  if (x.is_zero())
1189  return _ex1;
1190 
1191  // cosh(float) -> float
1193  return cosh(ex_to<numeric>(x));
1194 
1195  // cosh() is even
1197  return cosh(-x);
1198  }
1199 
1200  if ((x/Pi).info(info_flags::numeric) &&
1201  ex_to<numeric>(x/Pi).real().is_zero()) // cosh(I*x) -> cos(x)
1202  return cos(x/I);
1203 
1204  if (is_exactly_a<function>(x)) {
1205  const ex &t = x.op(0);
1206 
1207  // cosh(acosh(x)) -> x
1208  if (is_ex_the_function(x, acosh))
1209  return t;
1210 
1211  // cosh(asinh(x)) -> sqrt(1+x^2)
1212  if (is_ex_the_function(x, asinh))
1213  return sqrt(_ex1+power(t,_ex2));
1214 
1215  // cosh(atanh(x)) -> 1/sqrt(1-x^2)
1216  if (is_ex_the_function(x, atanh))
1217  return power(_ex1-power(t,_ex2),_ex_1_2);
1218  }
1219 
1220  return cosh(x).hold();
1221 }
1222 
1223 static ex cosh_deriv(const ex & x, unsigned deriv_param)
1224 {
1225  GINAC_ASSERT(deriv_param==0);
1226 
1227  // d/dx cosh(x) -> sinh(x)
1228  return sinh(x);
1229 }
1230 
1231 static ex cosh_real_part(const ex & x)
1232 {
1234 }
1235 
1236 static ex cosh_imag_part(const ex & x)
1237 {
1239 }
1240 
1241 static ex cosh_conjugate(const ex & x)
1242 {
1243  // conjugate(cosh(x))==cosh(conjugate(x))
1244  return cosh(x.conjugate());
1245 }
1246 
1248  evalf_func(cosh_evalf).
1249  derivative_func(cosh_deriv).
1250  real_part_func(cosh_real_part).
1251  imag_part_func(cosh_imag_part).
1252  conjugate_func(cosh_conjugate).
1253  latex_name("\\cosh"));
1254 
1256 // hyperbolic tangent (trigonometric function)
1258 
1259 static ex tanh_evalf(const ex & x)
1260 {
1261  if (is_exactly_a<numeric>(x))
1262  return tanh(ex_to<numeric>(x));
1263 
1264  return tanh(x).hold();
1265 }
1266 
1267 static ex tanh_eval(const ex & x)
1268 {
1269  if (x.info(info_flags::numeric)) {
1270 
1271  // tanh(0) -> 0
1272  if (x.is_zero())
1273  return _ex0;
1274 
1275  // tanh(float) -> float
1277  return tanh(ex_to<numeric>(x));
1278 
1279  // tanh() is odd
1281  return -tanh(-x);
1282  }
1283 
1284  if ((x/Pi).info(info_flags::numeric) &&
1285  ex_to<numeric>(x/Pi).real().is_zero()) // tanh(I*x) -> I*tan(x);
1286  return I*tan(x/I);
1287 
1288  if (is_exactly_a<function>(x)) {
1289  const ex &t = x.op(0);
1290 
1291  // tanh(atanh(x)) -> x
1292  if (is_ex_the_function(x, atanh))
1293  return t;
1294 
1295  // tanh(asinh(x)) -> x/sqrt(1+x^2)
1296  if (is_ex_the_function(x, asinh))
1297  return t*power(_ex1+power(t,_ex2),_ex_1_2);
1298 
1299  // tanh(acosh(x)) -> sqrt(x-1)*sqrt(x+1)/x
1300  if (is_ex_the_function(x, acosh))
1301  return sqrt(t-_ex1)*sqrt(t+_ex1)*power(t,_ex_1);
1302  }
1303 
1304  return tanh(x).hold();
1305 }
1306 
1307 static ex tanh_deriv(const ex & x, unsigned deriv_param)
1308 {
1309  GINAC_ASSERT(deriv_param==0);
1310 
1311  // d/dx tanh(x) -> 1-tanh(x)^2
1312  return _ex1-power(tanh(x),_ex2);
1313 }
1314 
1315 static ex tanh_series(const ex &x,
1316  const relational &rel,
1317  int order,
1318  unsigned options)
1319 {
1320  GINAC_ASSERT(is_a<symbol>(rel.lhs()));
1321  // method:
1322  // Taylor series where there is no pole falls back to tanh_deriv.
1323  // On a pole simply expand sinh(x)/cosh(x).
1324  const ex x_pt = x.subs(rel, subs_options::no_pattern);
1325  if (!(2*I*x_pt/Pi).info(info_flags::odd))
1326  throw do_taylor(); // caught by function::series()
1327  // if we got here we have to care for a simple pole
1328  return (sinh(x)/cosh(x)).series(rel, order, options);
1329 }
1330 
1331 static ex tanh_real_part(const ex & x)
1332 {
1333  ex a = GiNaC::real_part(x);
1334  ex b = GiNaC::imag_part(x);
1335  return tanh(a)/(1+power(tanh(a),2)*power(tan(b),2));
1336 }
1337 
1338 static ex tanh_imag_part(const ex & x)
1339 {
1340  ex a = GiNaC::real_part(x);
1341  ex b = GiNaC::imag_part(x);
1342  return tan(b)/(1+power(tanh(a),2)*power(tan(b),2));
1343 }
1344 
1345 static ex tanh_conjugate(const ex & x)
1346 {
1347  // conjugate(tanh(x))==tanh(conjugate(x))
1348  return tanh(x.conjugate());
1349 }
1350 
1352  evalf_func(tanh_evalf).
1353  derivative_func(tanh_deriv).
1354  series_func(tanh_series).
1355  real_part_func(tanh_real_part).
1356  imag_part_func(tanh_imag_part).
1357  conjugate_func(tanh_conjugate).
1358  latex_name("\\tanh"));
1359 
1361 // inverse hyperbolic sine (trigonometric function)
1363 
1364 static ex asinh_evalf(const ex & x)
1365 {
1366  if (is_exactly_a<numeric>(x))
1367  return asinh(ex_to<numeric>(x));
1368 
1369  return asinh(x).hold();
1370 }
1371 
1372 static ex asinh_eval(const ex & x)
1373 {
1374  if (x.info(info_flags::numeric)) {
1375 
1376  // asinh(0) -> 0
1377  if (x.is_zero())
1378  return _ex0;
1379 
1380  // asinh(float) -> float
1382  return asinh(ex_to<numeric>(x));
1383 
1384  // asinh() is odd
1386  return -asinh(-x);
1387  }
1388 
1389  return asinh(x).hold();
1390 }
1391 
1392 static ex asinh_deriv(const ex & x, unsigned deriv_param)
1393 {
1394  GINAC_ASSERT(deriv_param==0);
1395 
1396  // d/dx asinh(x) -> 1/sqrt(1+x^2)
1397  return power(_ex1+power(x,_ex2),_ex_1_2);
1398 }
1399 
1400 static ex asinh_conjugate(const ex & x)
1401 {
1402  // conjugate(asinh(x))==asinh(conjugate(x)) unless on the branch cuts which
1403  // run along the imaginary axis outside the interval [-I, +I].
1404  if (x.info(info_flags::real))
1405  return asinh(x);
1406  if (is_exactly_a<numeric>(x)) {
1407  const numeric x_re = ex_to<numeric>(x.real_part());
1408  const numeric x_im = ex_to<numeric>(x.imag_part());
1409  if (!x_re.is_zero() ||
1410  (x_im > *_num_1_p && x_im < *_num1_p))
1411  return asinh(x.conjugate());
1412  }
1413  return conjugate_function(asinh(x)).hold();
1414 }
1415 
1417  evalf_func(asinh_evalf).
1418  derivative_func(asinh_deriv).
1419  conjugate_func(asinh_conjugate));
1420 
1422 // inverse hyperbolic cosine (trigonometric function)
1424 
1425 static ex acosh_evalf(const ex & x)
1426 {
1427  if (is_exactly_a<numeric>(x))
1428  return acosh(ex_to<numeric>(x));
1429 
1430  return acosh(x).hold();
1431 }
1432 
1433 static ex acosh_eval(const ex & x)
1434 {
1435  if (x.info(info_flags::numeric)) {
1436 
1437  // acosh(0) -> Pi*I/2
1438  if (x.is_zero())
1439  return Pi*I*numeric(1,2);
1440 
1441  // acosh(1) -> 0
1442  if (x.is_equal(_ex1))
1443  return _ex0;
1444 
1445  // acosh(-1) -> Pi*I
1446  if (x.is_equal(_ex_1))
1447  return Pi*I;
1448 
1449  // acosh(float) -> float
1451  return acosh(ex_to<numeric>(x));
1452 
1453  // acosh(-x) -> Pi*I-acosh(x)
1455  return Pi*I-acosh(-x);
1456  }
1457 
1458  return acosh(x).hold();
1459 }
1460 
1461 static ex acosh_deriv(const ex & x, unsigned deriv_param)
1462 {
1463  GINAC_ASSERT(deriv_param==0);
1464 
1465  // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
1466  return power(x+_ex_1,_ex_1_2)*power(x+_ex1,_ex_1_2);
1467 }
1468 
1469 static ex acosh_conjugate(const ex & x)
1470 {
1471  // conjugate(acosh(x))==acosh(conjugate(x)) unless on the branch cut
1472  // which runs along the real axis from +1 to -inf.
1473  if (is_exactly_a<numeric>(x) &&
1474  (!x.imag_part().is_zero() || x > *_num1_p)) {
1475  return acosh(x.conjugate());
1476  }
1477  return conjugate_function(acosh(x)).hold();
1478 }
1479 
1481  evalf_func(acosh_evalf).
1482  derivative_func(acosh_deriv).
1483  conjugate_func(acosh_conjugate));
1484 
1486 // inverse hyperbolic tangent (trigonometric function)
1488 
1489 static ex atanh_evalf(const ex & x)
1490 {
1491  if (is_exactly_a<numeric>(x))
1492  return atanh(ex_to<numeric>(x));
1493 
1494  return atanh(x).hold();
1495 }
1496 
1497 static ex atanh_eval(const ex & x)
1498 {
1499  if (x.info(info_flags::numeric)) {
1500 
1501  // atanh(0) -> 0
1502  if (x.is_zero())
1503  return _ex0;
1504 
1505  // atanh({+|-}1) -> throw
1506  if (x.is_equal(_ex1) || x.is_equal(_ex_1))
1507  throw (pole_error("atanh_eval(): logarithmic pole",0));
1508 
1509  // atanh(float) -> float
1511  return atanh(ex_to<numeric>(x));
1512 
1513  // atanh() is odd
1515  return -atanh(-x);
1516  }
1517 
1518  return atanh(x).hold();
1519 }
1520 
1521 static ex atanh_deriv(const ex & x, unsigned deriv_param)
1522 {
1523  GINAC_ASSERT(deriv_param==0);
1524 
1525  // d/dx atanh(x) -> 1/(1-x^2)
1526  return power(_ex1-power(x,_ex2),_ex_1);
1527 }
1528 
1529 static ex atanh_series(const ex &arg,
1530  const relational &rel,
1531  int order,
1532  unsigned options)
1533 {
1534  GINAC_ASSERT(is_a<symbol>(rel.lhs()));
1535  // method:
1536  // Taylor series where there is no pole or cut falls back to atanh_deriv.
1537  // There are two branch cuts, one runnig from 1 up the real axis and one
1538  // one running from -1 down the real axis. The points 1 and -1 are poles
1539  // On the branch cuts and the poles series expand
1540  // (log(1+x)-log(1-x))/2
1541  // instead.
1542  const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
1543  if (!(arg_pt).info(info_flags::real))
1544  throw do_taylor(); // Im(x) != 0
1545  if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1)
1546  throw do_taylor(); // Im(x) == 0, but abs(x)<1
1547  // care for the poles, using the defining formula for atanh()...
1548  if (arg_pt.is_equal(_ex1) || arg_pt.is_equal(_ex_1))
1549  return ((log(_ex1+arg)-log(_ex1-arg))*_ex1_2).series(rel, order, options);
1550  // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
1552  // method:
1553  // This is the branch cut: assemble the primitive series manually and
1554  // then add the corresponding complex step function.
1555  const symbol &s = ex_to<symbol>(rel.lhs());
1556  const ex &point = rel.rhs();
1557  const symbol foo;
1558  const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
1559  ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2;
1560  if (arg_pt<_ex0)
1561  Order0correction += log((arg_pt+_ex_1)/(arg_pt+_ex1))*_ex1_2;
1562  else
1563  Order0correction += log((arg_pt+_ex1)/(arg_pt+_ex_1))*_ex_1_2;
1564  epvector seq;
1565  if (order > 0) {
1566  seq.reserve(2);
1567  seq.push_back(expair(Order0correction, _ex0));
1568  }
1569  seq.push_back(expair(Order(_ex1), order));
1570  return series(replarg - pseries(rel, std::move(seq)), rel, order);
1571  }
1572  throw do_taylor();
1573 }
1574 
1575 static ex atanh_conjugate(const ex & x)
1576 {
1577  // conjugate(atanh(x))==atanh(conjugate(x)) unless on the branch cuts which
1578  // run along the real axis outside the interval [-1, +1].
1579  if (is_exactly_a<numeric>(x) &&
1580  (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) {
1581  return atanh(x.conjugate());
1582  }
1583  return conjugate_function(atanh(x)).hold();
1584 }
1585 
1587  evalf_func(atanh_evalf).
1588  derivative_func(atanh_deriv).
1589  series_func(atanh_series).
1590  conjugate_func(atanh_conjugate));
1591 
1592 
1593 } // namespace GiNaC
Interface to GiNaC's sums of expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
Sum of expressions.
Definition: add.h:32
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
Exception class thrown by classes which provide their own series expansion to signal that ordinary Ta...
Definition: function.h:668
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
const_iterator begin() const noexcept
Definition: ex.h:647
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
Definition: ex.cpp:86
ex expand(unsigned options=0) const
Definition: ex.cpp:73
bool is_equal(const ex &other) const
Definition: ex.h:345
ex conjugate() const
Definition: ex.h:146
const_iterator end() const noexcept
Definition: ex.h:652
size_t nops() const
Definition: ex.h:135
ex series(const ex &r, int order, unsigned options=0) const
Compute the truncated series expansion of an expression.
Definition: pseries.cpp:1259
ex imag_part() const
Definition: ex.h:148
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
bool info(unsigned inf) const
Definition: ex.h:132
bool is_zero() const
Definition: ex.h:213
ex op(size_t i) const
Definition: ex.h:136
ex real_part() const
Definition: ex.h:147
A pair of expressions.
Definition: expair.h:38
@ expand_transcendental
expands transcendental functions like log and exp
Definition: flags.h:35
@ expand_function_args
expands the arguments of functions
Definition: flags.h:33
Product of expressions.
Definition: mul.h:32
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
bool is_equal(const numeric &other) const
Definition: numeric.cpp:1122
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
Exception class thrown when a singularity is encountered.
Definition: numeric.h:70
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
This class holds a extended truncated power series (positive and negative integer powers).
Definition: pseries.h:36
int ldegree(const ex &s) const override
Return degree of lowest power of the series.
Definition: pseries.cpp:335
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in power series if s is the expansion variable.
Definition: pseries.cpp:357
size_t nops() const override
Return the number of operands including a possible order term.
Definition: pseries.cpp:296
bool is_terminating() const
Returns true if there is no order term, i.e.
Definition: pseries.cpp:586
ex add_series(const pseries &other) const
Add one series object to another, producing a pseries object that represents the sum.
Definition: pseries.cpp:681
This class holds a relation consisting of two expressions and a logical relation between them.
Definition: relational.h:35
ex rhs() const
Definition: relational.h:82
ex lhs() const
Definition: relational.h:81
@ suppress_branchcut
Suppress branch cuts in series expansion.
Definition: flags.h:83
@ expanded
.expand(0) has already done its job (other expand() options ignore this flag)
Definition: flags.h:204
@ no_pattern
disable pattern matching
Definition: flags.h:51
Basic CAS symbol.
Definition: symbol.h:39
Interface to GiNaC's constant types and some special constants.
Interface to GiNaC's light-weight expression handles.
unsigned options
Definition: factor.cpp:2480
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
#define is_ex_the_function(OBJ, FUNCNAME)
Definition: function.h:765
Interface to GiNaC's initially known functions.
int order
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
static ex atan_eval(const ex &x)
bool is_zero(const ex &thisex)
Definition: ex.h:820
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
const numeric atan(const numeric &x)
Numeric arcustangent.
Definition: numeric.cpp:1508
ex real_part(const ex &thisex)
Definition: ex.h:721
static ex acosh_eval(const ex &x)
const numeric * _num_1_p
Definition: utils.cpp:351
static ex atanh_series(const ex &arg, const relational &rel, int order, unsigned options)
const ex _ex2
Definition: utils.cpp:389
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
const ex _ex_1_2
Definition: utils.cpp:356
const numeric * _num6_p
Definition: utils.cpp:404
static ex atanh_deriv(const ex &x, unsigned deriv_param)
static ex log_eval(const ex &x)
static ex tanh_eval(const ex &x)
const ex _ex1_2
Definition: utils.cpp:381
static ex asinh_eval(const ex &x)
static ex atanh_conjugate(const ex &x)
static ex asin_evalf(const ex &x)
const numeric cosh(const numeric &x)
Numeric hyperbolic cosine (trigonometric function).
Definition: numeric.cpp:1563
static ex atan_deriv(const ex &x, unsigned deriv_param)
static ex tan_eval(const ex &x)
const numeric * _num3_p
Definition: utils.cpp:392
std::vector< expair > epvector
expair-vector
Definition: expairseq.h:33
static ex cosh_eval(const ex &x)
const numeric mod(const numeric &a, const numeric &b)
Modulus (in positive representation).
Definition: numeric.cpp:2328
static ex log_deriv(const ex &x, unsigned deriv_param)
static ex asin_deriv(const ex &x, unsigned deriv_param)
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
static ex sinh_conjugate(const ex &x)
const numeric asin(const numeric &x)
Numeric inverse sine (trigonometric function).
Definition: numeric.cpp:1488
const numeric * _num25_p
Definition: utils.cpp:448
const ex _ex1
Definition: utils.cpp:385
const numeric tanh(const numeric &x)
Numeric hyperbolic tangent (trigonometric function).
Definition: numeric.cpp:1572
static ex tan_evalf(const ex &x)
int csgn(const numeric &x)
Definition: numeric.h:260
const numeric acos(const numeric &x)
Numeric inverse cosine (trigonometric function).
Definition: numeric.cpp:1497
static ex tan_real_part(const ex &x)
static ex atan_conjugate(const ex &x)
const numeric * _num24_p
Definition: utils.cpp:444
static ex tan_deriv(const ex &x, unsigned deriv_param)
const numeric sqrt(const numeric &x)
Numeric square root.
Definition: numeric.cpp:2475
static ex acos_conjugate(const ex &x)
const numeric * _num30_p
Definition: utils.cpp:452
const ex _ex3
Definition: utils.cpp:393
const numeric * _num10_p
Definition: utils.cpp:420
ex series(const ex &thisex, const ex &r, int order, unsigned options=0)
Definition: ex.h:781
const ex _ex6
Definition: utils.cpp:405
const numeric sinh(const numeric &x)
Numeric hyperbolic sine (trigonometric function).
Definition: numeric.cpp:1554
const numeric * _num4_p
Definition: utils.cpp:396
static ex asin_eval(const ex &x)
static ex atan_evalf(const ex &x)
const numeric * _num2_p
Definition: utils.cpp:388
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
static ex cosh_conjugate(const ex &x)
static ex atan2_eval(const ex &y, const ex &x)
static ex sinh_deriv(const ex &x, unsigned deriv_param)
static ex sin_conjugate(const ex &x)
static ex exp_real_part(const ex &x)
static ex tanh_deriv(const ex &x, unsigned deriv_param)
static ex acosh_evalf(const ex &x)
const numeric acosh(const numeric &x)
Numeric inverse hyperbolic cosine (trigonometric function).
Definition: numeric.cpp:1590
const numeric cos(const numeric &x)
Numeric cosine (trigonometric function).
Definition: numeric.cpp:1470
const ex _ex_1
Definition: utils.cpp:352
static ex cos_imag_part(const ex &x)
static ex atanh_evalf(const ex &x)
const numeric * _num12_p
Definition: utils.cpp:428
static ex cos_conjugate(const ex &x)
static ex acosh_deriv(const ex &x, unsigned deriv_param)
static ex acos_evalf(const ex &x)
const constant Pi("Pi", PiEvalf, "\\pi", domain::positive)
Pi.
Definition: constant.h:82
static ex exp_eval(const ex &x)
static ex cosh_evalf(const ex &x)
const numeric atanh(const numeric &x)
Numeric inverse hyperbolic tangent (trigonometric function).
Definition: numeric.cpp:1599
static ex exp_expand(const ex &arg, unsigned options)
static ex cosh_real_part(const ex &x)
static ex atan_series(const ex &arg, const relational &rel, int order, unsigned options)
static ex acos_eval(const ex &x)
static ex sinh_imag_part(const ex &x)
const ex _ex5
Definition: utils.cpp:401
static ex cos_real_part(const ex &x)
static ex tanh_conjugate(const ex &x)
const numeric log(const numeric &x)
Natural logarithm.
Definition: numeric.cpp:1450
const numeric real(const numeric &x)
Definition: numeric.h:311
const numeric sin(const numeric &x)
Numeric sine (trigonometric function).
Definition: numeric.cpp:1461
static ex cos_evalf(const ex &x)
static ex atan2_evalf(const ex &y, const ex &x)
static ex atan2_deriv(const ex &y, const ex &x, unsigned deriv_param)
const numeric * _num1_p
Definition: utils.cpp:384
static ex sin_evalf(const ex &x)
static ex log_real_part(const ex &x)
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:742
static ex exp_power(const ex &x, const ex &a)
static ex sin_deriv(const ex &x, unsigned deriv_param)
static ex cos_eval(const ex &x)
static ex asinh_evalf(const ex &x)
static ex sinh_eval(const ex &x)
static ex tanh_evalf(const ex &x)
const numeric * _num60_p
Definition: utils.cpp:460
const ex _ex_1_4
Definition: utils.cpp:364
static ex sinh_evalf(const ex &x)
static ex exp_imag_part(const ex &x)
const numeric asinh(const numeric &x)
Numeric inverse hyperbolic sine (trigonometric function).
Definition: numeric.cpp:1581
static ex sin_imag_part(const ex &x)
const ex _ex60
Definition: utils.cpp:461
const numeric tan(const numeric &x)
Numeric tangent (trigonometric function).
Definition: numeric.cpp:1479
static ex tan_series(const ex &x, const relational &rel, int order, unsigned options)
const ex _ex1_4
Definition: utils.cpp:373
static ex tan_conjugate(const ex &x)
static ex tan_imag_part(const ex &x)
static ex log_expand(const ex &arg, unsigned options)
static ex log_evalf(const ex &x)
static ex cosh_imag_part(const ex &x)
static ex exp_evalf(const ex &x)
static ex acos_deriv(const ex &x, unsigned deriv_param)
static ex log_conjugate(const ex &x)
REGISTER_FUNCTION(conjugate_function, eval_func(conjugate_eval). evalf_func(conjugate_evalf). expl_derivative_func(conjugate_expl_derivative). info_func(conjugate_info). print_func< print_latex >(conjugate_print_latex). conjugate_func(conjugate_conjugate). real_part_func(conjugate_real_part). imag_part_func(conjugate_imag_part). set_name("conjugate","conjugate"))
static ex sinh_real_part(const ex &x)
const numeric * _num15_p
Definition: utils.cpp:432
static ex log_series(const ex &arg, const relational &rel, int order, unsigned options)
static ex sin_eval(const ex &x)
static ex asinh_conjugate(const ex &x)
static ex asin_conjugate(const ex &x)
static ex cosh_deriv(const ex &x, unsigned deriv_param)
static ex acosh_conjugate(const ex &x)
const ex _ex_1_3
Definition: utils.cpp:360
const ex _ex0
Definition: utils.cpp:369
const numeric * _num120_p
Definition: utils.cpp:464
std::vector< ex > exvector
Definition: basic.h:46
static ex log_imag_part(const ex &x)
static ex tanh_series(const ex &x, const relational &rel, int order, unsigned options)
static ex tanh_real_part(const ex &x)
ex imag_part(const ex &thisex)
Definition: ex.h:724
static ex tanh_imag_part(const ex &x)
static ex cos_deriv(const ex &x, unsigned deriv_param)
static ex sin_real_part(const ex &x)
static ex exp_deriv(const ex &x, unsigned deriv_param)
static ex exp_conjugate(const ex &x)
const numeric * _num20_p
Definition: utils.cpp:440
static ex atanh_eval(const ex &x)
const numeric * _num18_p
Definition: utils.cpp:436
const numeric * _num0_p
Definition: utils.cpp:367
const numeric * _num5_p
Definition: utils.cpp:400
static ex asinh_deriv(const ex &x, unsigned deriv_param)
const ex _ex1_3
Definition: utils.cpp:377
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to class for extended truncated power series.
Interface to relations between expressions.
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.