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