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