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