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-2000 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                 // atan(float) -> float
 593                 if (!x.info(info_flags::crational))
 594                         return atan_evalf(x);
 595         }
 596
 597         return atan(x).hold();
 598 }
 599
 600 static ex atan_deriv(const ex & x, unsigned deriv_param)
 601 {
 602         GINAC_ASSERT(deriv_param==0);
 603
 604         // d/dx atan(x) -> 1/(1+x^2)
 605         return power(_ex1()+power(x,_ex2()), _ex_1());
 606 }
 607
 608 REGISTER_FUNCTION(atan, eval_func(atan_eval).
 609                         evalf_func(atan_evalf).
 610                         derivative_func(atan_deriv));
 611
 612 //////////
 613 // inverse tangent (atan2(y,x))
 614 //////////
 615
 616 static ex atan2_evalf(const ex & y, const ex & x)
 617 {
 618         BEGIN_TYPECHECK
 619                 TYPECHECK(y,numeric)
 620                 TYPECHECK(x,numeric)
 621         END_TYPECHECK(atan2(y,x))
 622
 623         return atan(ex_to_numeric(y),ex_to_numeric(x)); // -> numeric atan(numeric)
 624 }
 625
 626 static ex atan2_eval(const ex & y, const ex & x)
 627 {
 628         if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
 629                 x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
 630                 return atan2_evalf(y,x);
 631         }
 632
 633         return atan2(y,x).hold();
 634 }
 635
 636 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param)
 637 {
 638         GINAC_ASSERT(deriv_param<2);
 639
 640         if (deriv_param==0) {
 641                 // d/dy atan(y,x)
 642                 return x*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
 643         }
 644         // d/dx atan(y,x)
 645         return -y*power(power(x,_ex2())+power(y,_ex2()),_ex_1());
 646 }
 647
 648 REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
 649                          evalf_func(atan2_evalf).
 650                          derivative_func(atan2_deriv));
 651
 652 //////////
 653 // hyperbolic sine (trigonometric function)
 654 //////////
 655
 656 static ex sinh_evalf(const ex & x)
 657 {
 658         BEGIN_TYPECHECK
 659            TYPECHECK(x,numeric)
 660         END_TYPECHECK(sinh(x))
 661
 662         return sinh(ex_to_numeric(x)); // -> numeric sinh(numeric)
 663 }
 664
 665 static ex sinh_eval(const ex & x)
 666 {
 667         if (x.info(info_flags::numeric)) {
 668                 if (x.is_zero())  // sinh(0) -> 0
 669                         return _ex0();
 670                 if (!x.info(info_flags::crational))  // sinh(float) -> float
 671                         return sinh_evalf(x);
 672         }
 673
 674         if ((x/Pi).info(info_flags::numeric) &&
 675                 ex_to_numeric(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
 676                 return I*sin(x/I);
 677
 678         if (is_ex_exactly_of_type(x, function)) {
 679                 ex t = x.op(0);
 680                 // sinh(asinh(x)) -> x
 681                 if (is_ex_the_function(x, asinh))
 682                         return t;
 683                 // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2)
 684                 if (is_ex_the_function(x, acosh))
 685                         return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2());
 686                 // sinh(atanh(x)) -> x*(1-x^2)^(-1/2)
 687                 if (is_ex_the_function(x, atanh))
 688                         return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
 689         }
 690
 691         return sinh(x).hold();
 692 }
 693
 694 static ex sinh_deriv(const ex & x, unsigned deriv_param)
 695 {
 696         GINAC_ASSERT(deriv_param==0);
 697
 698         // d/dx sinh(x) -> cosh(x)
 699         return cosh(x);
 700 }
 701
 702 REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
 703                         evalf_func(sinh_evalf).
 704                         derivative_func(sinh_deriv));
 705
 706 //////////
 707 // hyperbolic cosine (trigonometric function)
 708 //////////
 709
 710 static ex cosh_evalf(const ex & x)
 711 {
 712         BEGIN_TYPECHECK
 713            TYPECHECK(x,numeric)
 714         END_TYPECHECK(cosh(x))
 715
 716         return cosh(ex_to_numeric(x)); // -> numeric cosh(numeric)
 717 }
 718
 719 static ex cosh_eval(const ex & x)
 720 {
 721         if (x.info(info_flags::numeric)) {
 722                 if (x.is_zero())  // cosh(0) -> 1
 723                         return _ex1();
 724                 if (!x.info(info_flags::crational))  // cosh(float) -> float
 725                         return cosh_evalf(x);
 726         }
 727
 728         if ((x/Pi).info(info_flags::numeric) &&
 729                 ex_to_numeric(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
 730                 return cos(x/I);
 731
 732         if (is_ex_exactly_of_type(x, function)) {
 733                 ex t = x.op(0);
 734                 // cosh(acosh(x)) -> x
 735                 if (is_ex_the_function(x, acosh))
 736                         return t;
 737                 // cosh(asinh(x)) -> (1+x^2)^(1/2)
 738                 if (is_ex_the_function(x, asinh))
 739                         return power(_ex1()+power(t,_ex2()),_ex1_2());
 740                 // cosh(atanh(x)) -> (1-x^2)^(-1/2)
 741                 if (is_ex_the_function(x, atanh))
 742                         return power(_ex1()-power(t,_ex2()),_ex_1_2());
 743         }
 744
 745         return cosh(x).hold();
 746 }
 747
 748 static ex cosh_deriv(const ex & x, unsigned deriv_param)
 749 {
 750         GINAC_ASSERT(deriv_param==0);
 751
 752         // d/dx cosh(x) -> sinh(x)
 753         return sinh(x);
 754 }
 755
 756 REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
 757                         evalf_func(cosh_evalf).
 758                         derivative_func(cosh_deriv));
 759
 760
 761 //////////
 762 // hyperbolic tangent (trigonometric function)
 763 //////////
 764
 765 static ex tanh_evalf(const ex & x)
 766 {
 767         BEGIN_TYPECHECK
 768            TYPECHECK(x,numeric)
 769         END_TYPECHECK(tanh(x))
 770
 771         return tanh(ex_to_numeric(x)); // -> numeric tanh(numeric)
 772 }
 773
 774 static ex tanh_eval(const ex & x)
 775 {
 776         if (x.info(info_flags::numeric)) {
 777                 if (x.is_zero())  // tanh(0) -> 0
 778                         return _ex0();
 779                 if (!x.info(info_flags::crational))  // tanh(float) -> float
 780                         return tanh_evalf(x);
 781         }
 782
 783         if ((x/Pi).info(info_flags::numeric) &&
 784                 ex_to_numeric(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
 785                 return I*tan(x/I);
 786
 787         if (is_ex_exactly_of_type(x, function)) {
 788                 ex t = x.op(0);
 789                 // tanh(atanh(x)) -> x
 790                 if (is_ex_the_function(x, atanh))
 791                         return t;
 792                 // tanh(asinh(x)) -> x*(1+x^2)^(-1/2)
 793                 if (is_ex_the_function(x, asinh))
 794                         return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
 795                 // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x
 796                 if (is_ex_the_function(x, acosh))
 797                         return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2())*power(t,_ex_1());
 798         }
 799
 800         return tanh(x).hold();
 801 }
 802
 803 static ex tanh_deriv(const ex & x, unsigned deriv_param)
 804 {
 805         GINAC_ASSERT(deriv_param==0);
 806
 807         // d/dx tanh(x) -> 1-tanh(x)^2
 808         return _ex1()-power(tanh(x),_ex2());
 809 }
 810
 811 static ex tanh_series(const ex &x,
 812                       const relational &rel,
 813                       int order,
 814                       unsigned options)
 815 {
 816         // method:
 817         // Taylor series where there is no pole falls back to tanh_deriv.
 818         // On a pole simply expand sinh(x)/cosh(x).
 819         const ex x_pt = x.subs(rel);
 820         if (!(2*I*x_pt/Pi).info(info_flags::odd))
 821                 throw do_taylor();  // caught by function::series()
 822         // if we got here we have to care for a simple pole
 823         return (sinh(x)/cosh(x)).series(rel, order+2);
 824 }
 825
 826 REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
 827                         evalf_func(tanh_evalf).
 828                         derivative_func(tanh_deriv).
 829                         series_func(tanh_series));
 830
 831 //////////
 832 // inverse hyperbolic sine (trigonometric function)
 833 //////////
 834
 835 static ex asinh_evalf(const ex & x)
 836 {
 837         BEGIN_TYPECHECK
 838            TYPECHECK(x,numeric)
 839         END_TYPECHECK(asinh(x))
 840
 841         return asinh(ex_to_numeric(x)); // -> numeric asinh(numeric)
 842 }
 843
 844 static ex asinh_eval(const ex & x)
 845 {
 846         if (x.info(info_flags::numeric)) {
 847                 // asinh(0) -> 0
 848                 if (x.is_zero())
 849                         return _ex0();
 850                 // asinh(float) -> float
 851                 if (!x.info(info_flags::crational))
 852                         return asinh_evalf(x);
 853         }
 854
 855         return asinh(x).hold();
 856 }
 857
 858 static ex asinh_deriv(const ex & x, unsigned deriv_param)
 859 {
 860         GINAC_ASSERT(deriv_param==0);
 861
 862         // d/dx asinh(x) -> 1/sqrt(1+x^2)
 863         return power(_ex1()+power(x,_ex2()),_ex_1_2());
 864 }
 865
 866 REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
 867                          evalf_func(asinh_evalf).
 868                          derivative_func(asinh_deriv));
 869
 870 //////////
 871 // inverse hyperbolic cosine (trigonometric function)
 872 //////////
 873
 874 static ex acosh_evalf(const ex & x)
 875 {
 876         BEGIN_TYPECHECK
 877            TYPECHECK(x,numeric)
 878         END_TYPECHECK(acosh(x))
 879
 880         return acosh(ex_to_numeric(x)); // -> numeric acosh(numeric)
 881 }
 882
 883 static ex acosh_eval(const ex & x)
 884 {
 885         if (x.info(info_flags::numeric)) {
 886                 // acosh(0) -> Pi*I/2
 887                 if (x.is_zero())
 888                         return Pi*I*numeric(1,2);
 889                 // acosh(1) -> 0
 890                 if (x.is_equal(_ex1()))
 891                         return _ex0();
 892                 // acosh(-1) -> Pi*I
 893                 if (x.is_equal(_ex_1()))
 894                         return Pi*I;
 895                 // acosh(float) -> float
 896                 if (!x.info(info_flags::crational))
 897                         return acosh_evalf(x);
 898         }
 899
 900         return acosh(x).hold();
 901 }
 902
 903 static ex acosh_deriv(const ex & x, unsigned deriv_param)
 904 {
 905         GINAC_ASSERT(deriv_param==0);
 906
 907         // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
 908         return power(x+_ex_1(),_ex_1_2())*power(x+_ex1(),_ex_1_2());
 909 }
 910
 911 REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
 912                          evalf_func(acosh_evalf).
 913                          derivative_func(acosh_deriv));
 914
 915 //////////
 916 // inverse hyperbolic tangent (trigonometric function)
 917 //////////
 918
 919 static ex atanh_evalf(const ex & x)
 920 {
 921         BEGIN_TYPECHECK
 922            TYPECHECK(x,numeric)
 923         END_TYPECHECK(atanh(x))
 924
 925         return atanh(ex_to_numeric(x)); // -> numeric atanh(numeric)
 926 }
 927
 928 static ex atanh_eval(const ex & x)
 929 {
 930         if (x.info(info_flags::numeric)) {
 931                 // atanh(0) -> 0
 932                 if (x.is_zero())
 933                         return _ex0();
 934                 // atanh({+|-}1) -> throw
 935                 if (x.is_equal(_ex1()) || x.is_equal(_ex_1()))
 936                         throw (pole_error("atanh_eval(): logarithmic pole",0));
 937                 // atanh(float) -> float
 938                 if (!x.info(info_flags::crational))
 939                         return atanh_evalf(x);
 940         }
 941
 942         return atanh(x).hold();
 943 }
 944
 945 static ex atanh_deriv(const ex & x, unsigned deriv_param)
 946 {
 947         GINAC_ASSERT(deriv_param==0);
 948
 949         // d/dx atanh(x) -> 1/(1-x^2)
 950         return power(_ex1()-power(x,_ex2()),_ex_1());
 951 }
 952
 953 REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
 954                          evalf_func(atanh_evalf).
 955                          derivative_func(atanh_deriv));
 956
 957 #ifndef NO_NAMESPACE_GINAC
 958 } // namespace GiNaC
 959 #endif // ndef NO_NAMESPACE_GINAC