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