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-2005 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 <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
  59         // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
  60         const ex TwoExOverPiI=(_ex2*x)/(Pi*I);
  61         if (TwoExOverPiI.info(info_flags::integer)) {
  62                 const numeric z = mod(ex_to<numeric>(TwoExOverPiI),*_num4_p);
  63                 if (z.is_equal(*_num0_p))
  64                         return _ex1;
  65                 if (z.is_equal(*_num1_p))
  66                         return ex(I);
  67                 if (z.is_equal(*_num2_p))
  68                         return _ex_1;
  69                 if (z.is_equal(*_num3_p))
  70                         return ex(-I);
  71         }
  72
  73         // exp(log(x)) -> x
  74         if (is_ex_the_function(x, log))
  75                 return x.op(0);
  76
  77         // exp(float) -> float
  78         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
  79                 return exp(ex_to<numeric>(x));
  80
  81         return exp(x).hold();
  82 }
  83
  84 static ex exp_deriv(const ex & x, unsigned deriv_param)
  85 {
  86         GINAC_ASSERT(deriv_param==0);
  87
  88         // d/dx exp(x) -> exp(x)
  89         return exp(x);
  90 }
  91
  92 REGISTER_FUNCTION(exp, eval_func(exp_eval).
  93                        evalf_func(exp_evalf).
  94                        derivative_func(exp_deriv).
  95                        latex_name("\\exp"));
  96
  97 //////////
  98 // natural logarithm
  99 //////////
 100
 101 static ex log_evalf(const ex & x)
 102 {
 103         if (is_exactly_a<numeric>(x))
 104                 return log(ex_to<numeric>(x));
 105
 106         return log(x).hold();
 107 }
 108
 109 static ex log_eval(const ex & x)
 110 {
 111         if (x.info(info_flags::numeric)) {
 112                 if (x.is_zero())         // log(0) -> infinity
 113                         throw(pole_error("log_eval(): log(0)",0));
 114                 if (x.info(info_flags::rational) && x.info(info_flags::negative))
 115                         return (log(-x)+I*Pi);
 116                 if (x.is_equal(_ex1))  // log(1) -> 0
 117                         return _ex0;
 118                 if (x.is_equal(I))       // log(I) -> Pi*I/2
 119                         return (Pi*I*_ex1_2);
 120                 if (x.is_equal(-I))      // log(-I) -> -Pi*I/2
 121                         return (Pi*I*_ex_1_2);
 122
 123                 // log(float) -> float
 124                 if (!x.info(info_flags::crational))
 125                         return log(ex_to<numeric>(x));
 126         }
 127
 128         // log(exp(t)) -> t (if -Pi < t.imag() <= Pi):
 129         if (is_ex_the_function(x, exp)) {
 130                 const ex &t = x.op(0);
 131                 if (is_a<symbol>(t) && t.info(info_flags::real)) {
 132                         return t;
 133                 }
 134                 if (t.info(info_flags::numeric)) {
 135                         const numeric &nt = ex_to<numeric>(t);
 136                         if (nt.is_real())
 137                                 return t;
 138                 }
 139         }
 140
 141         return log(x).hold();
 142 }
 143
 144 static ex log_deriv(const ex & x, unsigned deriv_param)
 145 {
 146         GINAC_ASSERT(deriv_param==0);
 147
 148         // d/dx log(x) -> 1/x
 149         return power(x, _ex_1);
 150 }
 151
 152 static ex log_series(const ex &arg,
 153                      const relational &rel,
 154                      int order,
 155                      unsigned options)
 156 {
 157         GINAC_ASSERT(is_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, subs_options::no_pattern);
 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                         if (n == 0 && coeff == 1) {
 204                                 epvector epv;
 205                                 ex acc = (new pseries(rel, epv))->setflag(status_flags::dynallocated);
 206                                 epv.reserve(2);
 207                                 epv.push_back(expair(-1, _ex0));
 208                                 epv.push_back(expair(Order(_ex1), order));
 209                                 ex rest = pseries(rel, epv).add_series(argser);
 210                                 for (int i = order-1; i>0; --i) {
 211                                         epvector cterm;
 212                                         cterm.reserve(1);
 213                                         cterm.push_back(expair(i%2 ? _ex1/i : _ex_1/i, _ex0));
 214                                         acc = pseries(rel, cterm).add_series(ex_to<pseries>(acc));
 215                                         acc = (ex_to<pseries>(rest)).mul_series(ex_to<pseries>(acc));
 216                                 }
 217                                 return acc;
 218                         }
 219                         const ex newarg = ex_to<pseries>((arg/coeff).series(rel, order+n, options)).shift_exponents(-n).convert_to_poly(true);
 220                         return pseries(rel, seq).add_series(ex_to<pseries>(log(newarg).series(rel, order, options)));
 221                 } else  // it was a monomial
 222                         return pseries(rel, seq);
 223         }
 224         if (!(options & series_options::suppress_branchcut) &&
 225              arg_pt.info(info_flags::negative)) {
 226                 // method:
 227                 // This is the branch cut: assemble the primitive series manually and
 228                 // then add the corresponding complex step function.
 229                 const symbol &s = ex_to<symbol>(rel.lhs());
 230                 const ex &point = rel.rhs();
 231                 const symbol foo;
 232                 const ex replarg = series(log(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
 233                 epvector seq;
 234                 seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0));
 235                 seq.push_back(expair(Order(_ex1), order));
 236                 return series(replarg - I*Pi + pseries(rel, seq), rel, order);
 237         }
 238         throw do_taylor();  // caught by function::series()
 239 }
 240
 241 REGISTER_FUNCTION(log, eval_func(log_eval).
 242                        evalf_func(log_evalf).
 243                        derivative_func(log_deriv).
 244                        series_func(log_series).
 245                        latex_name("\\ln"));
 246
 247 //////////
 248 // sine (trigonometric function)
 249 //////////
 250
 251 static ex sin_evalf(const ex & x)
 252 {
 253         if (is_exactly_a<numeric>(x))
 254                 return sin(ex_to<numeric>(x));
 255
 256         return sin(x).hold();
 257 }
 258
 259 static ex sin_eval(const ex & x)
 260 {
 261         // sin(n/d*Pi) -> { all known non-nested radicals }
 262         const ex SixtyExOverPi = _ex60*x/Pi;
 263         ex sign = _ex1;
 264         if (SixtyExOverPi.info(info_flags::integer)) {
 265                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p);
 266                 if (z>=*_num60_p) {
 267                         // wrap to interval [0, Pi)
 268                         z -= *_num60_p;
 269                         sign = _ex_1;
 270                 }
 271                 if (z>*_num30_p) {
 272                         // wrap to interval [0, Pi/2)
 273                         z = *_num60_p-z;
 274                 }
 275                 if (z.is_equal(*_num0_p))  // sin(0)       -> 0
 276                         return _ex0;
 277                 if (z.is_equal(*_num5_p))  // sin(Pi/12)   -> sqrt(6)/4*(1-sqrt(3)/3)
 278                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
 279                 if (z.is_equal(*_num6_p))  // sin(Pi/10)   -> sqrt(5)/4-1/4
 280                         return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
 281                 if (z.is_equal(*_num10_p)) // sin(Pi/6)    -> 1/2
 282                         return sign*_ex1_2;
 283                 if (z.is_equal(*_num15_p)) // sin(Pi/4)    -> sqrt(2)/2
 284                         return sign*_ex1_2*sqrt(_ex2);
 285                 if (z.is_equal(*_num18_p)) // sin(3/10*Pi) -> sqrt(5)/4+1/4
 286                         return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
 287                 if (z.is_equal(*_num20_p)) // sin(Pi/3)    -> sqrt(3)/2
 288                         return sign*_ex1_2*sqrt(_ex3);
 289                 if (z.is_equal(*_num25_p)) // sin(5/12*Pi) -> sqrt(6)/4*(1+sqrt(3)/3)
 290                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
 291                 if (z.is_equal(*_num30_p)) // sin(Pi/2)    -> 1
 292                         return sign;
 293         }
 294
 295         if (is_exactly_a<function>(x)) {
 296                 const ex &t = x.op(0);
 297
 298                 // sin(asin(x)) -> x
 299                 if (is_ex_the_function(x, asin))
 300                         return t;
 301
 302                 // sin(acos(x)) -> sqrt(1-x^2)
 303                 if (is_ex_the_function(x, acos))
 304                         return sqrt(_ex1-power(t,_ex2));
 305
 306                 // sin(atan(x)) -> x/sqrt(1+x^2)
 307                 if (is_ex_the_function(x, atan))
 308                         return t*power(_ex1+power(t,_ex2),_ex_1_2);
 309         }
 310
 311         // sin(float) -> float
 312         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
 313                 return sin(ex_to<numeric>(x));
 314
 315         // sin() is odd
 316         if (x.info(info_flags::negative))
 317                 return -sin(-x);
 318
 319         return sin(x).hold();
 320 }
 321
 322 static ex sin_deriv(const ex & x, unsigned deriv_param)
 323 {
 324         GINAC_ASSERT(deriv_param==0);
 325
 326         // d/dx sin(x) -> cos(x)
 327         return cos(x);
 328 }
 329
 330 REGISTER_FUNCTION(sin, eval_func(sin_eval).
 331                        evalf_func(sin_evalf).
 332                        derivative_func(sin_deriv).
 333                        latex_name("\\sin"));
 334
 335 //////////
 336 // cosine (trigonometric function)
 337 //////////
 338
 339 static ex cos_evalf(const ex & x)
 340 {
 341         if (is_exactly_a<numeric>(x))
 342                 return cos(ex_to<numeric>(x));
 343
 344         return cos(x).hold();
 345 }
 346
 347 static ex cos_eval(const ex & x)
 348 {
 349         // cos(n/d*Pi) -> { all known non-nested radicals }
 350         const ex SixtyExOverPi = _ex60*x/Pi;
 351         ex sign = _ex1;
 352         if (SixtyExOverPi.info(info_flags::integer)) {
 353                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p);
 354                 if (z>=*_num60_p) {
 355                         // wrap to interval [0, Pi)
 356                         z = *_num120_p-z;
 357                 }
 358                 if (z>=*_num30_p) {
 359                         // wrap to interval [0, Pi/2)
 360                         z = *_num60_p-z;
 361                         sign = _ex_1;
 362                 }
 363                 if (z.is_equal(*_num0_p))  // cos(0)       -> 1
 364                         return sign;
 365                 if (z.is_equal(*_num5_p))  // cos(Pi/12)   -> sqrt(6)/4*(1+sqrt(3)/3)
 366                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex1_3*sqrt(_ex3));
 367                 if (z.is_equal(*_num10_p)) // cos(Pi/6)    -> sqrt(3)/2
 368                         return sign*_ex1_2*sqrt(_ex3);
 369                 if (z.is_equal(*_num12_p)) // cos(Pi/5)    -> sqrt(5)/4+1/4
 370                         return sign*(_ex1_4*sqrt(_ex5)+_ex1_4);
 371                 if (z.is_equal(*_num15_p)) // cos(Pi/4)    -> sqrt(2)/2
 372                         return sign*_ex1_2*sqrt(_ex2);
 373                 if (z.is_equal(*_num20_p)) // cos(Pi/3)    -> 1/2
 374                         return sign*_ex1_2;
 375                 if (z.is_equal(*_num24_p)) // cos(2/5*Pi)  -> sqrt(5)/4-1/4x
 376                         return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4);
 377                 if (z.is_equal(*_num25_p)) // cos(5/12*Pi) -> sqrt(6)/4*(1-sqrt(3)/3)
 378                         return sign*_ex1_4*sqrt(_ex6)*(_ex1+_ex_1_3*sqrt(_ex3));
 379                 if (z.is_equal(*_num30_p)) // cos(Pi/2)    -> 0
 380                         return _ex0;
 381         }
 382
 383         if (is_exactly_a<function>(x)) {
 384                 const ex &t = x.op(0);
 385
 386                 // cos(acos(x)) -> x
 387                 if (is_ex_the_function(x, acos))
 388                         return t;
 389
 390                 // cos(asin(x)) -> sqrt(1-x^2)
 391                 if (is_ex_the_function(x, asin))
 392                         return sqrt(_ex1-power(t,_ex2));
 393
 394                 // cos(atan(x)) -> 1/sqrt(1+x^2)
 395                 if (is_ex_the_function(x, atan))
 396                         return power(_ex1+power(t,_ex2),_ex_1_2);
 397         }
 398
 399         // cos(float) -> float
 400         if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
 401                 return cos(ex_to<numeric>(x));
 402
 403         // cos() is even
 404         if (x.info(info_flags::negative))
 405                 return cos(-x);
 406
 407         return cos(x).hold();
 408 }
 409
 410 static ex cos_deriv(const ex & x, unsigned deriv_param)
 411 {
 412         GINAC_ASSERT(deriv_param==0);
 413
 414         // d/dx cos(x) -> -sin(x)
 415         return -sin(x);
 416 }
 417
 418 REGISTER_FUNCTION(cos, eval_func(cos_eval).
 419                        evalf_func(cos_evalf).
 420                        derivative_func(cos_deriv).
 421                        latex_name("\\cos"));
 422
 423 //////////
 424 // tangent (trigonometric function)
 425 //////////
 426
 427 static ex tan_evalf(const ex & x)
 428 {
 429         if (is_exactly_a<numeric>(x))
 430                 return tan(ex_to<numeric>(x));
 431
 432         return tan(x).hold();
 433 }
 434
 435 static ex tan_eval(const ex & x)
 436 {
 437         // tan(n/d*Pi) -> { all known non-nested radicals }
 438         const ex SixtyExOverPi = _ex60*x/Pi;
 439         ex sign = _ex1;
 440         if (SixtyExOverPi.info(info_flags::integer)) {
 441                 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num60_p);
 442                 if (z>=*_num60_p) {
 443                         // wrap to interval [0, Pi)
 444                         z -= *_num60_p;
 445                 }
 446                 if (z>=*_num30_p) {
 447                         // wrap to interval [0, Pi/2)
 448                         z = *_num60_p-z;
 449                         sign = _ex_1;
 450                 }
 451                 if (z.is_equal(*_num0_p))  // tan(0)       -> 0
 452                         return _ex0;
 453                 if (z.is_equal(*_num5_p))  // tan(Pi/12)   -> 2-sqrt(3)
 454                         return sign*(_ex2-sqrt(_ex3));
 455                 if (z.is_equal(*_num10_p)) // tan(Pi/6)    -> sqrt(3)/3
 456                         return sign*_ex1_3*sqrt(_ex3);
 457                 if (z.is_equal(*_num15_p)) // tan(Pi/4)    -> 1
 458                         return sign;
 459                 if (z.is_equal(*_num20_p)) // tan(Pi/3)    -> sqrt(3)
 460                         return sign*sqrt(_ex3);
 461                 if (z.is_equal(*_num25_p)) // tan(5/12*Pi) -> 2+sqrt(3)
 462                         return sign*(sqrt(_ex3)+_ex2);
 463                 if (z.is_equal(*_num30_p)) // tan(Pi/2)    -> infinity
 464                         throw (pole_error("tan_eval(): simple pole",1));
 465         }
 466
 467         if (is_exactly_a<function>(x)) {
 468                 const ex &t = x.op(0);
 469
 470                 // tan(atan(x)) -> x
 471                 if (is_ex_the_function(x, atan))
 472                         return t;
 473
 474                 // tan(asin(x)) -> x/sqrt(1+x^2)
 475                 if (is_ex_the_function(x, asin))
 476                         return t*power(_ex1-power(t,_ex2),_ex_1_2);
 477
 478                 // tan(acos(x)) -> sqrt(1-x^2)/x
 479                 if (is_ex_the_function(x, acos))
 480                         return power(t,_ex_1)*sqrt(_ex1-power(t,_ex2));
 481         }
 482
 483         // tan(float) -> float
 484         if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
 485                 return tan(ex_to<numeric>(x));
 486         }
 487
 488         // tan() is odd
 489         if (x.info(info_flags::negative))
 490                 return -tan(-x);
 491
 492         return tan(x).hold();
 493 }
 494
 495 static ex tan_deriv(const ex & x, unsigned deriv_param)
 496 {
 497         GINAC_ASSERT(deriv_param==0);
 498
 499         // d/dx tan(x) -> 1+tan(x)^2;
 500         return (_ex1+power(tan(x),_ex2));
 501 }
 502
 503 static ex tan_series(const ex &x,
 504                      const relational &rel,
 505                      int order,
 506                      unsigned options)
 507 {
 508         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
 509         // method:
 510         // Taylor series where there is no pole falls back to tan_deriv.
 511         // On a pole simply expand sin(x)/cos(x).
 512         const ex x_pt = x.subs(rel, subs_options::no_pattern);
 513         if (!(2*x_pt/Pi).info(info_flags::odd))
 514                 throw do_taylor();  // caught by function::series()
 515         // if we got here we have to care for a simple pole
 516         return (sin(x)/cos(x)).series(rel, order, options);
 517 }
 518
 519 REGISTER_FUNCTION(tan, eval_func(tan_eval).
 520                        evalf_func(tan_evalf).
 521                        derivative_func(tan_deriv).
 522                        series_func(tan_series).
 523                        latex_name("\\tan"));
 524
 525 //////////
 526 // inverse sine (arc sine)
 527 //////////
 528
 529 static ex asin_evalf(const ex & x)
 530 {
 531         if (is_exactly_a<numeric>(x))
 532                 return asin(ex_to<numeric>(x));
 533
 534         return asin(x).hold();
 535 }
 536
 537 static ex asin_eval(const ex & x)
 538 {
 539         if (x.info(info_flags::numeric)) {
 540
 541                 // asin(0) -> 0
 542                 if (x.is_zero())
 543                         return x;
 544
 545                 // asin(1/2) -> Pi/6
 546                 if (x.is_equal(_ex1_2))
 547                         return numeric(1,6)*Pi;
 548
 549                 // asin(1) -> Pi/2
 550                 if (x.is_equal(_ex1))
 551                         return _ex1_2*Pi;
 552
 553                 // asin(-1/2) -> -Pi/6
 554                 if (x.is_equal(_ex_1_2))
 555                         return numeric(-1,6)*Pi;
 556
 557                 // asin(-1) -> -Pi/2
 558                 if (x.is_equal(_ex_1))
 559                         return _ex_1_2*Pi;
 560
 561                 // asin(float) -> float
 562                 if (!x.info(info_flags::crational))
 563                         return asin(ex_to<numeric>(x));
 564
 565                 // asin() is odd
 566                 if (x.info(info_flags::negative))
 567                         return -asin(-x);
 568         }
 569
 570         return asin(x).hold();
 571 }
 572
 573 static ex asin_deriv(const ex & x, unsigned deriv_param)
 574 {
 575         GINAC_ASSERT(deriv_param==0);
 576
 577         // d/dx asin(x) -> 1/sqrt(1-x^2)
 578         return power(1-power(x,_ex2),_ex_1_2);
 579 }
 580
 581 REGISTER_FUNCTION(asin, eval_func(asin_eval).
 582                         evalf_func(asin_evalf).
 583                         derivative_func(asin_deriv).
 584                         latex_name("\\arcsin"));
 585
 586 //////////
 587 // inverse cosine (arc cosine)
 588 //////////
 589
 590 static ex acos_evalf(const ex & x)
 591 {
 592         if (is_exactly_a<numeric>(x))
 593                 return acos(ex_to<numeric>(x));
 594
 595         return acos(x).hold();
 596 }
 597
 598 static ex acos_eval(const ex & x)
 599 {
 600         if (x.info(info_flags::numeric)) {
 601
 602                 // acos(1) -> 0
 603                 if (x.is_equal(_ex1))
 604                         return _ex0;
 605
 606                 // acos(1/2) -> Pi/3
 607                 if (x.is_equal(_ex1_2))
 608                         return _ex1_3*Pi;
 609
 610                 // acos(0) -> Pi/2
 611                 if (x.is_zero())
 612                         return _ex1_2*Pi;
 613
 614                 // acos(-1/2) -> 2/3*Pi
 615                 if (x.is_equal(_ex_1_2))
 616                         return numeric(2,3)*Pi;
 617
 618                 // acos(-1) -> Pi
 619                 if (x.is_equal(_ex_1))
 620                         return Pi;
 621
 622                 // acos(float) -> float
 623                 if (!x.info(info_flags::crational))
 624                         return acos(ex_to<numeric>(x));
 625
 626                 // acos(-x) -> Pi-acos(x)
 627                 if (x.info(info_flags::negative))
 628                         return Pi-acos(-x);
 629         }
 630
 631         return acos(x).hold();
 632 }
 633
 634 static ex acos_deriv(const ex & x, unsigned deriv_param)
 635 {
 636         GINAC_ASSERT(deriv_param==0);
 637
 638         // d/dx acos(x) -> -1/sqrt(1-x^2)
 639         return -power(1-power(x,_ex2),_ex_1_2);
 640 }
 641
 642 REGISTER_FUNCTION(acos, eval_func(acos_eval).
 643                         evalf_func(acos_evalf).
 644                         derivative_func(acos_deriv).
 645                         latex_name("\\arccos"));
 646
 647 //////////
 648 // inverse tangent (arc tangent)
 649 //////////
 650
 651 static ex atan_evalf(const ex & x)
 652 {
 653         if (is_exactly_a<numeric>(x))
 654                 return atan(ex_to<numeric>(x));
 655
 656         return atan(x).hold();
 657 }
 658
 659 static ex atan_eval(const ex & x)
 660 {
 661         if (x.info(info_flags::numeric)) {
 662
 663                 // atan(0) -> 0
 664                 if (x.is_zero())
 665                         return _ex0;
 666
 667                 // atan(1) -> Pi/4
 668                 if (x.is_equal(_ex1))
 669                         return _ex1_4*Pi;
 670
 671                 // atan(-1) -> -Pi/4
 672                 if (x.is_equal(_ex_1))
 673                         return _ex_1_4*Pi;
 674
 675                 if (x.is_equal(I) || x.is_equal(-I))
 676                         throw (pole_error("atan_eval(): logarithmic pole",0));
 677
 678                 // atan(float) -> float
 679                 if (!x.info(info_flags::crational))
 680                         return atan(ex_to<numeric>(x));
 681
 682                 // atan() is odd
 683                 if (x.info(info_flags::negative))
 684                         return -atan(-x);
 685         }
 686
 687         return atan(x).hold();
 688 }
 689
 690 static ex atan_deriv(const ex & x, unsigned deriv_param)
 691 {
 692         GINAC_ASSERT(deriv_param==0);
 693
 694         // d/dx atan(x) -> 1/(1+x^2)
 695         return power(_ex1+power(x,_ex2), _ex_1);
 696 }
 697
 698 static ex atan_series(const ex &arg,
 699                       const relational &rel,
 700                       int order,
 701                       unsigned options)
 702 {
 703         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
 704         // method:
 705         // Taylor series where there is no pole or cut falls back to atan_deriv.
 706         // There are two branch cuts, one runnig from I up the imaginary axis and
 707         // one running from -I down the imaginary axis.  The points I and -I are
 708         // poles.
 709         // On the branch cuts and the poles series expand
 710         //     (log(1+I*x)-log(1-I*x))/(2*I)
 711         // instead.
 712         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
 713         if (!(I*arg_pt).info(info_flags::real))
 714                 throw do_taylor();     // Re(x) != 0
 715         if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1)
 716                 throw do_taylor();     // Re(x) == 0, but abs(x)<1
 717         // care for the poles, using the defining formula for atan()...
 718         if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
 719                 return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
 720         if (!(options & series_options::suppress_branchcut)) {
 721                 // method:
 722                 // This is the branch cut: assemble the primitive series manually and
 723                 // then add the corresponding complex step function.
 724                 const symbol &s = ex_to<symbol>(rel.lhs());
 725                 const ex &point = rel.rhs();
 726                 const symbol foo;
 727                 const ex replarg = series(atan(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
 728                 ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2;
 729                 if ((I*arg_pt)<_ex0)
 730                         Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2;
 731                 else
 732                         Order0correction += log((I*arg_pt+_ex1)/(I*arg_pt+_ex_1))*I*_ex1_2;
 733                 epvector seq;
 734                 seq.push_back(expair(Order0correction, _ex0));
 735                 seq.push_back(expair(Order(_ex1), order));
 736                 return series(replarg - pseries(rel, seq), rel, order);
 737         }
 738         throw do_taylor();
 739 }
 740
 741 REGISTER_FUNCTION(atan, eval_func(atan_eval).
 742                         evalf_func(atan_evalf).
 743                         derivative_func(atan_deriv).
 744                         series_func(atan_series).
 745                         latex_name("\\arctan"));
 746
 747 //////////
 748 // inverse tangent (atan2(y,x))
 749 //////////
 750
 751 static ex atan2_evalf(const ex &y, const ex &x)
 752 {
 753         if (is_exactly_a<numeric>(y) && is_exactly_a<numeric>(x))
 754                 return atan(ex_to<numeric>(y), ex_to<numeric>(x));
 755
 756         return atan2(y, x).hold();
 757 }
 758
 759 static ex atan2_eval(const ex & y, const ex & x)
 760 {
 761         if (y.info(info_flags::numeric) && x.info(info_flags::numeric)) {
 762
 763                 if (y.is_zero()) {
 764
 765                         // atan(0, 0) -> 0
 766                         if (x.is_zero())
 767                                 return _ex0;
 768
 769                         // atan(0, x), x real and positive -> 0
 770                         if (x.info(info_flags::positive))
 771                                 return _ex0;
 772
 773                         // atan(0, x), x real and negative -> -Pi
 774                         if (x.info(info_flags::negative))
 775                                 return _ex_1*Pi;
 776                 }
 777
 778                 if (x.is_zero()) {
 779
 780                         // atan(y, 0), y real and positive -> Pi/2
 781                         if (y.info(info_flags::positive))
 782                                 return _ex1_2*Pi;
 783
 784                         // atan(y, 0), y real and negative -> -Pi/2
 785                         if (y.info(info_flags::negative))
 786                                 return _ex_1_2*Pi;
 787                 }
 788
 789                 if (y.is_equal(x)) {
 790
 791                         // atan(y, y), y real and positive -> Pi/4
 792                         if (y.info(info_flags::positive))
 793                                 return _ex1_4*Pi;
 794
 795                         // atan(y, y), y real and negative -> -3/4*Pi
 796                         if (y.info(info_flags::negative))
 797                                 return numeric(-3, 4)*Pi;
 798                 }
 799
 800                 if (y.is_equal(-x)) {
 801
 802                         // atan(y, -y), y real and positive -> 3*Pi/4
 803                         if (y.info(info_flags::positive))
 804                                 return numeric(3, 4)*Pi;
 805
 806                         // atan(y, -y), y real and negative -> -Pi/4
 807                         if (y.info(info_flags::negative))
 808                                 return _ex_1_4*Pi;
 809                 }
 810
 811                 // atan(float, float) -> float
 812                 if (!y.info(info_flags::crational) && !x.info(info_flags::crational))
 813                         return atan(ex_to<numeric>(y), ex_to<numeric>(x));
 814
 815                 // atan(real, real) -> atan(y/x) +/- Pi
 816                 if (y.info(info_flags::real) && x.info(info_flags::real)) {
 817                         if (x.info(info_flags::positive))
 818                                 return atan(y/x);
 819                         else if(y.info(info_flags::positive))
 820                                 return atan(y/x)+Pi;
 821                         else
 822                                 return atan(y/x)-Pi;
 823                 }
 824         }
 825
 826         return atan2(y, x).hold();
 827 }
 828
 829 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
 830 {
 831         GINAC_ASSERT(deriv_param<2);
 832
 833         if (deriv_param==0) {
 834                 // d/dy atan(y,x)
 835                 return x*power(power(x,_ex2)+power(y,_ex2),_ex_1);
 836         }
 837         // d/dx atan(y,x)
 838         return -y*power(power(x,_ex2)+power(y,_ex2),_ex_1);
 839 }
 840
 841 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
 842                          evalf_func(atan2_evalf).
 843                          derivative_func(atan2_deriv));
 844
 845 //////////
 846 // hyperbolic sine (trigonometric function)
 847 //////////
 848
 849 static ex sinh_evalf(const ex & x)
 850 {
 851         if (is_exactly_a<numeric>(x))
 852                 return sinh(ex_to<numeric>(x));
 853
 854         return sinh(x).hold();
 855 }
 856
 857 static ex sinh_eval(const ex & x)
 858 {
 859         if (x.info(info_flags::numeric)) {
 860
 861                 // sinh(0) -> 0
 862                 if (x.is_zero())
 863                         return _ex0;
 864
 865                 // sinh(float) -> float
 866                 if (!x.info(info_flags::crational))
 867                         return sinh(ex_to<numeric>(x));
 868
 869                 // sinh() is odd
 870                 if (x.info(info_flags::negative))
 871                         return -sinh(-x);
 872         }
 873
 874         if ((x/Pi).info(info_flags::numeric) &&
 875                 ex_to<numeric>(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
 876                 return I*sin(x/I);
 877
 878         if (is_exactly_a<function>(x)) {
 879                 const ex &t = x.op(0);
 880
 881                 // sinh(asinh(x)) -> x
 882                 if (is_ex_the_function(x, asinh))
 883                         return t;
 884
 885                 // sinh(acosh(x)) -> sqrt(x-1) * sqrt(x+1)
 886                 if (is_ex_the_function(x, acosh))
 887                         return sqrt(t-_ex1)*sqrt(t+_ex1);
 888
 889                 // sinh(atanh(x)) -> x/sqrt(1-x^2)
 890                 if (is_ex_the_function(x, atanh))
 891                         return t*power(_ex1-power(t,_ex2),_ex_1_2);
 892         }
 893
 894         return sinh(x).hold();
 895 }
 896
 897 static ex sinh_deriv(const ex & x, unsigned deriv_param)
 898 {
 899         GINAC_ASSERT(deriv_param==0);
 900
 901         // d/dx sinh(x) -> cosh(x)
 902         return cosh(x);
 903 }
 904
 905 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
 906                         evalf_func(sinh_evalf).
 907                         derivative_func(sinh_deriv).
 908                         latex_name("\\sinh"));
 909
 910 //////////
 911 // hyperbolic cosine (trigonometric function)
 912 //////////
 913
 914 static ex cosh_evalf(const ex & x)
 915 {
 916         if (is_exactly_a<numeric>(x))
 917                 return cosh(ex_to<numeric>(x));
 918
 919         return cosh(x).hold();
 920 }
 921
 922 static ex cosh_eval(const ex & x)
 923 {
 924         if (x.info(info_flags::numeric)) {
 925
 926                 // cosh(0) -> 1
 927                 if (x.is_zero())
 928                         return _ex1;
 929
 930                 // cosh(float) -> float
 931                 if (!x.info(info_flags::crational))
 932                         return cosh(ex_to<numeric>(x));
 933
 934                 // cosh() is even
 935                 if (x.info(info_flags::negative))
 936                         return cosh(-x);
 937         }
 938
 939         if ((x/Pi).info(info_flags::numeric) &&
 940                 ex_to<numeric>(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
 941                 return cos(x/I);
 942
 943         if (is_exactly_a<function>(x)) {
 944                 const ex &t = x.op(0);
 945
 946                 // cosh(acosh(x)) -> x
 947                 if (is_ex_the_function(x, acosh))
 948                         return t;
 949
 950                 // cosh(asinh(x)) -> sqrt(1+x^2)
 951                 if (is_ex_the_function(x, asinh))
 952                         return sqrt(_ex1+power(t,_ex2));
 953
 954                 // cosh(atanh(x)) -> 1/sqrt(1-x^2)
 955                 if (is_ex_the_function(x, atanh))
 956                         return power(_ex1-power(t,_ex2),_ex_1_2);
 957         }
 958
 959         return cosh(x).hold();
 960 }
 961
 962 static ex cosh_deriv(const ex & x, unsigned deriv_param)
 963 {
 964         GINAC_ASSERT(deriv_param==0);
 965
 966         // d/dx cosh(x) -> sinh(x)
 967         return sinh(x);
 968 }
 969
 970 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
 971                         evalf_func(cosh_evalf).
 972                         derivative_func(cosh_deriv).
 973                         latex_name("\\cosh"));
 974
 975 //////////
 976 // hyperbolic tangent (trigonometric function)
 977 //////////
 978
 979 static ex tanh_evalf(const ex & x)
 980 {
 981         if (is_exactly_a<numeric>(x))
 982                 return tanh(ex_to<numeric>(x));
 983
 984         return tanh(x).hold();
 985 }
 986
 987 static ex tanh_eval(const ex & x)
 988 {
 989         if (x.info(info_flags::numeric)) {
 990
 991                 // tanh(0) -> 0
 992                 if (x.is_zero())
 993                         return _ex0;
 994
 995                 // tanh(float) -> float
 996                 if (!x.info(info_flags::crational))
 997                         return tanh(ex_to<numeric>(x));
 998
 999                 // tanh() is odd
1000                 if (x.info(info_flags::negative))
1001                         return -tanh(-x);
1002         }
1003
1004         if ((x/Pi).info(info_flags::numeric) &&
1005                 ex_to<numeric>(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
1006                 return I*tan(x/I);
1007
1008         if (is_exactly_a<function>(x)) {
1009                 const ex &t = x.op(0);
1010
1011                 // tanh(atanh(x)) -> x
1012                 if (is_ex_the_function(x, atanh))
1013                         return t;
1014
1015                 // tanh(asinh(x)) -> x/sqrt(1+x^2)
1016                 if (is_ex_the_function(x, asinh))
1017                         return t*power(_ex1+power(t,_ex2),_ex_1_2);
1018
1019                 // tanh(acosh(x)) -> sqrt(x-1)*sqrt(x+1)/x
1020                 if (is_ex_the_function(x, acosh))
1021                         return sqrt(t-_ex1)*sqrt(t+_ex1)*power(t,_ex_1);
1022         }
1023
1024         return tanh(x).hold();
1025 }
1026
1027 static ex tanh_deriv(const ex & x, unsigned deriv_param)
1028 {
1029         GINAC_ASSERT(deriv_param==0);
1030
1031         // d/dx tanh(x) -> 1-tanh(x)^2
1032         return _ex1-power(tanh(x),_ex2);
1033 }
1034
1035 static ex tanh_series(const ex &x,
1036                       const relational &rel,
1037                       int order,
1038                       unsigned options)
1039 {
1040         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
1041         // method:
1042         // Taylor series where there is no pole falls back to tanh_deriv.
1043         // On a pole simply expand sinh(x)/cosh(x).
1044         const ex x_pt = x.subs(rel, subs_options::no_pattern);
1045         if (!(2*I*x_pt/Pi).info(info_flags::odd))
1046                 throw do_taylor();  // caught by function::series()
1047         // if we got here we have to care for a simple pole
1048         return (sinh(x)/cosh(x)).series(rel, order, options);
1049 }
1050
1051 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
1052                         evalf_func(tanh_evalf).
1053                         derivative_func(tanh_deriv).
1054                         series_func(tanh_series).
1055                         latex_name("\\tanh"));
1056
1057 //////////
1058 // inverse hyperbolic sine (trigonometric function)
1059 //////////
1060
1061 static ex asinh_evalf(const ex & x)
1062 {
1063         if (is_exactly_a<numeric>(x))
1064                 return asinh(ex_to<numeric>(x));
1065
1066         return asinh(x).hold();
1067 }
1068
1069 static ex asinh_eval(const ex & x)
1070 {
1071         if (x.info(info_flags::numeric)) {
1072
1073                 // asinh(0) -> 0
1074                 if (x.is_zero())
1075                         return _ex0;
1076
1077                 // asinh(float) -> float
1078                 if (!x.info(info_flags::crational))
1079                         return asinh(ex_to<numeric>(x));
1080
1081                 // asinh() is odd
1082                 if (x.info(info_flags::negative))
1083                         return -asinh(-x);
1084         }
1085
1086         return asinh(x).hold();
1087 }
1088
1089 static ex asinh_deriv(const ex & x, unsigned deriv_param)
1090 {
1091         GINAC_ASSERT(deriv_param==0);
1092
1093         // d/dx asinh(x) -> 1/sqrt(1+x^2)
1094         return power(_ex1+power(x,_ex2),_ex_1_2);
1095 }
1096
1097 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
1098                          evalf_func(asinh_evalf).
1099                          derivative_func(asinh_deriv));
1100
1101 //////////
1102 // inverse hyperbolic cosine (trigonometric function)
1103 //////////
1104
1105 static ex acosh_evalf(const ex & x)
1106 {
1107         if (is_exactly_a<numeric>(x))
1108                 return acosh(ex_to<numeric>(x));
1109
1110         return acosh(x).hold();
1111 }
1112
1113 static ex acosh_eval(const ex & x)
1114 {
1115         if (x.info(info_flags::numeric)) {
1116
1117                 // acosh(0) -> Pi*I/2
1118                 if (x.is_zero())
1119                         return Pi*I*numeric(1,2);
1120
1121                 // acosh(1) -> 0
1122                 if (x.is_equal(_ex1))
1123                         return _ex0;
1124
1125                 // acosh(-1) -> Pi*I
1126                 if (x.is_equal(_ex_1))
1127                         return Pi*I;
1128
1129                 // acosh(float) -> float
1130                 if (!x.info(info_flags::crational))
1131                         return acosh(ex_to<numeric>(x));
1132
1133                 // acosh(-x) -> Pi*I-acosh(x)
1134                 if (x.info(info_flags::negative))
1135                         return Pi*I-acosh(-x);
1136         }
1137
1138         return acosh(x).hold();
1139 }
1140
1141 static ex acosh_deriv(const ex & x, unsigned deriv_param)
1142 {
1143         GINAC_ASSERT(deriv_param==0);
1144
1145         // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
1146         return power(x+_ex_1,_ex_1_2)*power(x+_ex1,_ex_1_2);
1147 }
1148
1149 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
1150                          evalf_func(acosh_evalf).
1151                          derivative_func(acosh_deriv));
1152
1153 //////////
1154 // inverse hyperbolic tangent (trigonometric function)
1155 //////////
1156
1157 static ex atanh_evalf(const ex & x)
1158 {
1159         if (is_exactly_a<numeric>(x))
1160                 return atanh(ex_to<numeric>(x));
1161
1162         return atanh(x).hold();
1163 }
1164
1165 static ex atanh_eval(const ex & x)
1166 {
1167         if (x.info(info_flags::numeric)) {
1168
1169                 // atanh(0) -> 0
1170                 if (x.is_zero())
1171                         return _ex0;
1172
1173                 // atanh({+|-}1) -> throw
1174                 if (x.is_equal(_ex1) || x.is_equal(_ex_1))
1175                         throw (pole_error("atanh_eval(): logarithmic pole",0));
1176
1177                 // atanh(float) -> float
1178                 if (!x.info(info_flags::crational))
1179                         return atanh(ex_to<numeric>(x));
1180
1181                 // atanh() is odd
1182                 if (x.info(info_flags::negative))
1183                         return -atanh(-x);
1184         }
1185
1186         return atanh(x).hold();
1187 }
1188
1189 static ex atanh_deriv(const ex & x, unsigned deriv_param)
1190 {
1191         GINAC_ASSERT(deriv_param==0);
1192
1193         // d/dx atanh(x) -> 1/(1-x^2)
1194         return power(_ex1-power(x,_ex2),_ex_1);
1195 }
1196
1197 static ex atanh_series(const ex &arg,
1198                        const relational &rel,
1199                        int order,
1200                        unsigned options)
1201 {
1202         GINAC_ASSERT(is_a<symbol>(rel.lhs()));
1203         // method:
1204         // Taylor series where there is no pole or cut falls back to atanh_deriv.
1205         // There are two branch cuts, one runnig from 1 up the real axis and one
1206         // one running from -1 down the real axis.  The points 1 and -1 are poles
1207         // On the branch cuts and the poles series expand
1208         //     (log(1+x)-log(1-x))/2
1209         // instead.
1210         const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
1211         if (!(arg_pt).info(info_flags::real))
1212                 throw do_taylor();     // Im(x) != 0
1213         if ((arg_pt).info(info_flags::real) && abs(arg_pt)<_ex1)
1214                 throw do_taylor();     // Im(x) == 0, but abs(x)<1
1215         // care for the poles, using the defining formula for atanh()...
1216         if (arg_pt.is_equal(_ex1) || arg_pt.is_equal(_ex_1))
1217                 return ((log(_ex1+arg)-log(_ex1-arg))*_ex1_2).series(rel, order, options);
1218         // ...and the branch cuts (the discontinuity at the cut being just I*Pi)
1219         if (!(options & series_options::suppress_branchcut)) {
1220                 // method:
1221                 // This is the branch cut: assemble the primitive series manually and
1222                 // then add the corresponding complex step function.
1223                 const symbol &s = ex_to<symbol>(rel.lhs());
1224                 const ex &point = rel.rhs();
1225                 const symbol foo;
1226                 const ex replarg = series(atanh(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
1227                 ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2;
1228                 if (arg_pt<_ex0)
1229                         Order0correction += log((arg_pt+_ex_1)/(arg_pt+_ex1))*_ex1_2;
1230                 else
1231                         Order0correction += log((arg_pt+_ex1)/(arg_pt+_ex_1))*_ex_1_2;
1232                 epvector seq;
1233                 seq.push_back(expair(Order0correction, _ex0));
1234                 seq.push_back(expair(Order(_ex1), order));
1235                 return series(replarg - pseries(rel, seq), rel, order);
1236         }
1237         throw do_taylor();
1238 }
1239
1240 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
1241                          evalf_func(atanh_evalf).
1242                          derivative_func(atanh_deriv).
1243                          series_func(atanh_series));
1244
1245
1246 } // namespace GiNaC