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         // method:
 451         // Taylor series where there is no pole falls back to tan_deriv.
 452         // On a pole simply expand sin(x)/cos(x).
 453         const ex x_pt = x.subs(rel);
 454         if (!(2*x_pt/Pi).info(info_flags::odd))
 455                 throw do_taylor();  // caught by function::series()
 456         // if we got here we have to care for a simple pole
 457         return (sin(x)/cos(x)).series(rel, order+2);
 458 }
 459
 460 REGISTER_FUNCTION(tan, eval_func(tan_eval).
 461                        evalf_func(tan_evalf).
 462                        derivative_func(tan_deriv).
 463                        series_func(tan_series));
 464
 465 //////////
 466 // inverse sine (arc sine)
 467 //////////
 468
 469 static ex asin_evalf(const ex & x)
 470 {
 471         BEGIN_TYPECHECK
 472            TYPECHECK(x,numeric)
 473         END_TYPECHECK(asin(x))
 474
 475         return asin(ex_to_numeric(x)); // -> numeric asin(numeric)
 476 }
 477
 478 static ex asin_eval(const ex & x)
 479 {
 480         if (x.info(info_flags::numeric)) {
 481                 // asin(0) -> 0
 482                 if (x.is_zero())
 483                         return x;
 484                 // asin(1/2) -> Pi/6
 485                 if (x.is_equal(_ex1_2()))
 486                         return numeric(1,6)*Pi;
 487                 // asin(1) -> Pi/2
 488                 if (x.is_equal(_ex1()))
 489                         return _num1_2()*Pi;
 490                 // asin(-1/2) -> -Pi/6
 491                 if (x.is_equal(_ex_1_2()))
 492                         return numeric(-1,6)*Pi;
 493                 // asin(-1) -> -Pi/2
 494                 if (x.is_equal(_ex_1()))
 495                         return _num_1_2()*Pi;
 496                 // asin(float) -> float
 497                 if (!x.info(info_flags::crational))
 498                         return asin_evalf(x);
 499         }
 500
 501         return asin(x).hold();
 502 }
 503
 504 static ex asin_deriv(const ex & x, unsigned deriv_param)
 505 {
 506         GINAC_ASSERT(deriv_param==0);
 507
 508         // d/dx asin(x) -> 1/sqrt(1-x^2)
 509         return power(1-power(x,_ex2()),_ex_1_2());
 510 }
 511
 512 REGISTER_FUNCTION(asin, eval_func(asin_eval).
 513                         evalf_func(asin_evalf).
 514                         derivative_func(asin_deriv));
 515
 516 //////////
 517 // inverse cosine (arc cosine)
 518 //////////
 519
 520 static ex acos_evalf(const ex & x)
 521 {
 522         BEGIN_TYPECHECK
 523            TYPECHECK(x,numeric)
 524         END_TYPECHECK(acos(x))
 525
 526         return acos(ex_to_numeric(x)); // -> numeric acos(numeric)
 527 }
 528
 529 static ex acos_eval(const ex & x)
 530 {
 531         if (x.info(info_flags::numeric)) {
 532                 // acos(1) -> 0
 533                 if (x.is_equal(_ex1()))
 534                         return _ex0();
 535                 // acos(1/2) -> Pi/3
 536                 if (x.is_equal(_ex1_2()))
 537                         return _ex1_3()*Pi;
 538                 // acos(0) -> Pi/2
 539                 if (x.is_zero())
 540                         return _ex1_2()*Pi;
 541                 // acos(-1/2) -> 2/3*Pi
 542                 if (x.is_equal(_ex_1_2()))
 543                         return numeric(2,3)*Pi;
 544                 // acos(-1) -> Pi
 545                 if (x.is_equal(_ex_1()))
 546                         return Pi;
 547                 // acos(float) -> float
 548                 if (!x.info(info_flags::crational))
 549                         return acos_evalf(x);
 550         }
 551
 552         return acos(x).hold();
 553 }
 554
 555 static ex acos_deriv(const ex & x, unsigned deriv_param)
 556 {
 557         GINAC_ASSERT(deriv_param==0);
 558
 559         // d/dx acos(x) -> -1/sqrt(1-x^2)
 560         return _ex_1()*power(1-power(x,_ex2()),_ex_1_2());
 561 }
 562
 563 REGISTER_FUNCTION(acos, eval_func(acos_eval).
 564                         evalf_func(acos_evalf).
 565                         derivative_func(acos_deriv));
 566
 567 //////////
 568 // inverse tangent (arc tangent)
 569 //////////
 570
 571 static ex atan_evalf(const ex & x)
 572 {
 573         BEGIN_TYPECHECK
 574                 TYPECHECK(x,numeric)
 575         END_TYPECHECK(atan(x))
 576
 577         return atan(ex_to_numeric(x)); // -> numeric atan(numeric)
 578 }
 579
 580 static ex atan_eval(const ex & x)
 581 {
 582         if (x.info(info_flags::numeric)) {
 583                 // atan(0) -> 0
 584                 if (x.is_equal(_ex0()))
 585                         return _ex0();
 586                 // atan(1) -> Pi/4
 587                 if (x.is_equal(_ex1()))
 588                         return _ex1_4()*Pi;
 589                 // atan(-1) -> -Pi/4
 590                 if (x.is_equal(_ex_1()))
 591                         return _ex_1_4()*Pi;
 592                 if (x.is_equal(I) || x.is_equal(-I))
 593                         throw (pole_error("atan_eval(): logarithmic pole",0));
 594                 // atan(float) -> float
 595                 if (!x.info(info_flags::crational))
 596                         return atan_evalf(x);
 597         }
 598
 599         return atan(x).hold();
 600 }
 601
 602 static ex atan_deriv(const ex & x, unsigned deriv_param)
 603 {
 604         GINAC_ASSERT(deriv_param==0);
 605
 606         // d/dx atan(x) -> 1/(1+x^2)
 607         return power(_ex1()+power(x,_ex2()), _ex_1());
 608 }
 609
 610 REGISTER_FUNCTION(atan, eval_func(atan_eval).
 611                         evalf_func(atan_evalf).
 612                         derivative_func(atan_deriv));
 613
 614 //////////
 615 // inverse tangent (atan2(y,x))
 616 //////////
 617
 618 static ex atan2_evalf(const ex & y, const ex & x)
 619 {
 620         BEGIN_TYPECHECK
 621                 TYPECHECK(y,numeric)
 622                 TYPECHECK(x,numeric)
 623         END_TYPECHECK(atan2(y,x))
 624
 625         return atan(ex_to_numeric(y),ex_to_numeric(x)); // -> numeric atan(numeric)
 626 }
 627
 628 static ex atan2_eval(const ex & y, const ex & x)
 629 {
 630         if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
 631                 x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
 632                 return atan2_evalf(y,x);
 633         }
 634
 635         return atan2(y,x).hold();
 636 }
 637
 638 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
 639 {
 640         GINAC_ASSERT(deriv_param<2);
 641
 642         if (deriv_param==0) {
 643                 // d/dy atan(y,x)
 644                 return x*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
 645         }
 646         // d/dx atan(y,x)
 647         return -y*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
 648 }
 649
 650 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
 651                          evalf_func(atan2_evalf).
 652                          derivative_func(atan2_deriv));
 653
 654 //////////
 655 // hyperbolic sine (trigonometric function)
 656 //////////
 657
 658 static ex sinh_evalf(const ex & x)
 659 {
 660         BEGIN_TYPECHECK
 661            TYPECHECK(x,numeric)
 662         END_TYPECHECK(sinh(x))
 663
 664         return sinh(ex_to_numeric(x)); // -> numeric sinh(numeric)
 665 }
 666
 667 static ex sinh_eval(const ex & x)
 668 {
 669         if (x.info(info_flags::numeric)) {
 670                 if (x.is_zero())  // sinh(0) -> 0
 671                         return _ex0();
 672                 if (!x.info(info_flags::crational))  // sinh(float) -> float
 673                         return sinh_evalf(x);
 674         }
 675
 676         if ((x/Pi).info(info_flags::numeric) &&
 677                 ex_to_numeric(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
 678                 return I*sin(x/I);
 679
 680         if (is_ex_exactly_of_type(x, function)) {
 681                 ex t = x.op(0);
 682                 // sinh(asinh(x)) -> x
 683                 if (is_ex_the_function(x, asinh))
 684                         return t;
 685                 // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2)
 686                 if (is_ex_the_function(x, acosh))
 687                         return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2());
 688                 // sinh(atanh(x)) -> x*(1-x^2)^(-1/2)
 689                 if (is_ex_the_function(x, atanh))
 690                         return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
 691         }
 692
 693         return sinh(x).hold();
 694 }
 695
 696 static ex sinh_deriv(const ex & x, unsigned deriv_param)
 697 {
 698         GINAC_ASSERT(deriv_param==0);
 699
 700         // d/dx sinh(x) -> cosh(x)
 701         return cosh(x);
 702 }
 703
 704 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
 705                         evalf_func(sinh_evalf).
 706                         derivative_func(sinh_deriv));
 707
 708 //////////
 709 // hyperbolic cosine (trigonometric function)
 710 //////////
 711
 712 static ex cosh_evalf(const ex & x)
 713 {
 714         BEGIN_TYPECHECK
 715            TYPECHECK(x,numeric)
 716         END_TYPECHECK(cosh(x))
 717
 718         return cosh(ex_to_numeric(x)); // -> numeric cosh(numeric)
 719 }
 720
 721 static ex cosh_eval(const ex & x)
 722 {
 723         if (x.info(info_flags::numeric)) {
 724                 if (x.is_zero())  // cosh(0) -> 1
 725                         return _ex1();
 726                 if (!x.info(info_flags::crational))  // cosh(float) -> float
 727                         return cosh_evalf(x);
 728         }
 729
 730         if ((x/Pi).info(info_flags::numeric) &&
 731                 ex_to_numeric(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
 732                 return cos(x/I);
 733
 734         if (is_ex_exactly_of_type(x, function)) {
 735                 ex t = x.op(0);
 736                 // cosh(acosh(x)) -> x
 737                 if (is_ex_the_function(x, acosh))
 738                         return t;
 739                 // cosh(asinh(x)) -> (1+x^2)^(1/2)
 740                 if (is_ex_the_function(x, asinh))
 741                         return power(_ex1()+power(t,_ex2()),_ex1_2());
 742                 // cosh(atanh(x)) -> (1-x^2)^(-1/2)
 743                 if (is_ex_the_function(x, atanh))
 744                         return power(_ex1()-power(t,_ex2()),_ex_1_2());
 745         }
 746
 747         return cosh(x).hold();
 748 }
 749
 750 static ex cosh_deriv(const ex & x, unsigned deriv_param)
 751 {
 752         GINAC_ASSERT(deriv_param==0);
 753
 754         // d/dx cosh(x) -> sinh(x)
 755         return sinh(x);
 756 }
 757
 758 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
 759                         evalf_func(cosh_evalf).
 760                         derivative_func(cosh_deriv));
 761
 762
 763 //////////
 764 // hyperbolic tangent (trigonometric function)
 765 //////////
 766
 767 static ex tanh_evalf(const ex & x)
 768 {
 769         BEGIN_TYPECHECK
 770            TYPECHECK(x,numeric)
 771         END_TYPECHECK(tanh(x))
 772
 773         return tanh(ex_to_numeric(x)); // -> numeric tanh(numeric)
 774 }
 775
 776 static ex tanh_eval(const ex & x)
 777 {
 778         if (x.info(info_flags::numeric)) {
 779                 if (x.is_zero())  // tanh(0) -> 0
 780                         return _ex0();
 781                 if (!x.info(info_flags::crational))  // tanh(float) -> float
 782                         return tanh_evalf(x);
 783         }
 784
 785         if ((x/Pi).info(info_flags::numeric) &&
 786                 ex_to_numeric(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
 787                 return I*tan(x/I);
 788
 789         if (is_ex_exactly_of_type(x, function)) {
 790                 ex t = x.op(0);
 791                 // tanh(atanh(x)) -> x
 792                 if (is_ex_the_function(x, atanh))
 793                         return t;
 794                 // tanh(asinh(x)) -> x*(1+x^2)^(-1/2)
 795                 if (is_ex_the_function(x, asinh))
 796                         return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
 797                 // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x
 798                 if (is_ex_the_function(x, acosh))
 799                         return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2())*power(t,_ex_1());
 800         }
 801
 802         return tanh(x).hold();
 803 }
 804
 805 static ex tanh_deriv(const ex & x, unsigned deriv_param)
 806 {
 807         GINAC_ASSERT(deriv_param==0);
 808
 809         // d/dx tanh(x) -> 1-tanh(x)^2
 810         return _ex1()-power(tanh(x),_ex2());
 811 }
 812
 813 static ex tanh_series(const ex &x,
 814                       const relational &rel,
 815                       int order,
 816                       unsigned options)
 817 {
 818         // method:
 819         // Taylor series where there is no pole falls back to tanh_deriv.
 820         // On a pole simply expand sinh(x)/cosh(x).
 821         const ex x_pt = x.subs(rel);
 822         if (!(2*I*x_pt/Pi).info(info_flags::odd))
 823                 throw do_taylor();  // caught by function::series()
 824         // if we got here we have to care for a simple pole
 825         return (sinh(x)/cosh(x)).series(rel, order+2);
 826 }
 827
 828 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
 829                         evalf_func(tanh_evalf).
 830                         derivative_func(tanh_deriv).
 831                         series_func(tanh_series));
 832
 833 //////////
 834 // inverse hyperbolic sine (trigonometric function)
 835 //////////
 836
 837 static ex asinh_evalf(const ex & x)
 838 {
 839         BEGIN_TYPECHECK
 840            TYPECHECK(x,numeric)
 841         END_TYPECHECK(asinh(x))
 842
 843         return asinh(ex_to_numeric(x)); // -> numeric asinh(numeric)
 844 }
 845
 846 static ex asinh_eval(const ex & x)
 847 {
 848         if (x.info(info_flags::numeric)) {
 849                 // asinh(0) -> 0
 850                 if (x.is_zero())
 851                         return _ex0();
 852                 // asinh(float) -> float
 853                 if (!x.info(info_flags::crational))
 854                         return asinh_evalf(x);
 855         }
 856
 857         return asinh(x).hold();
 858 }
 859
 860 static ex asinh_deriv(const ex & x, unsigned deriv_param)
 861 {
 862         GINAC_ASSERT(deriv_param==0);
 863
 864         // d/dx asinh(x) -> 1/sqrt(1+x^2)
 865         return power(_ex1()+power(x,_ex2()),_ex_1_2());
 866 }
 867
 868 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
 869                          evalf_func(asinh_evalf).
 870                          derivative_func(asinh_deriv));
 871
 872 //////////
 873 // inverse hyperbolic cosine (trigonometric function)
 874 //////////
 875
 876 static ex acosh_evalf(const ex & x)
 877 {
 878         BEGIN_TYPECHECK
 879            TYPECHECK(x,numeric)
 880         END_TYPECHECK(acosh(x))
 881
 882         return acosh(ex_to_numeric(x)); // -> numeric acosh(numeric)
 883 }
 884
 885 static ex acosh_eval(const ex & x)
 886 {
 887         if (x.info(info_flags::numeric)) {
 888                 // acosh(0) -> Pi*I/2
 889                 if (x.is_zero())
 890                         return Pi*I*numeric(1,2);
 891                 // acosh(1) -> 0
 892                 if (x.is_equal(_ex1()))
 893                         return _ex0();
 894                 // acosh(-1) -> Pi*I
 895                 if (x.is_equal(_ex_1()))
 896                         return Pi*I;
 897                 // acosh(float) -> float
 898                 if (!x.info(info_flags::crational))
 899                         return acosh_evalf(x);
 900         }
 901
 902         return acosh(x).hold();
 903 }
 904
 905 static ex acosh_deriv(const ex & x, unsigned deriv_param)
 906 {
 907         GINAC_ASSERT(deriv_param==0);
 908
 909         // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
 910         return power(x+_ex_1(),_ex_1_2())*power(x+_ex1(),_ex_1_2());
 911 }
 912
 913 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
 914                          evalf_func(acosh_evalf).
 915                          derivative_func(acosh_deriv));
 916
 917 //////////
 918 // inverse hyperbolic tangent (trigonometric function)
 919 //////////
 920
 921 static ex atanh_evalf(const ex & x)
 922 {
 923         BEGIN_TYPECHECK
 924            TYPECHECK(x,numeric)
 925         END_TYPECHECK(atanh(x))
 926
 927         return atanh(ex_to_numeric(x)); // -> numeric atanh(numeric)
 928 }
 929
 930 static ex atanh_eval(const ex & x)
 931 {
 932         if (x.info(info_flags::numeric)) {
 933                 // atanh(0) -> 0
 934                 if (x.is_zero())
 935                         return _ex0();
 936                 // atanh({+|-}1) -> throw
 937                 if (x.is_equal(_ex1()) || x.is_equal(_ex_1()))
 938                         throw (pole_error("atanh_eval(): logarithmic pole",0));
 939                 // atanh(float) -> float
 940                 if (!x.info(info_flags::crational))
 941                         return atanh_evalf(x);
 942         }
 943
 944         return atanh(x).hold();
 945 }
 946
 947 static ex atanh_deriv(const ex & x, unsigned deriv_param)
 948 {
 949         GINAC_ASSERT(deriv_param==0);
 950
 951         // d/dx atanh(x) -> 1/(1-x^2)
 952         return power(_ex1()-power(x,_ex2()),_ex_1());
 953 }
 954
 955 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
 956                          evalf_func(atanh_evalf).
 957                          derivative_func(atanh_deriv));
 958
 959 #ifndef NO_NAMESPACE_GINAC
 960 } // namespace GiNaC
 961 #endif // ndef NO_NAMESPACE_GINAC