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