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 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 "utils.h"
  35
  36 #ifndef NO_GINAC_NAMESPACE
  37 namespace GiNaC {
  38 #endif // ndef NO_GINAC_NAMESPACE
  39
  40 //////////
  41 // exponential function
  42 //////////
  43
  44 static ex exp_evalf(ex const & x)
  45 {
  46     BEGIN_TYPECHECK
  47         TYPECHECK(x,numeric)
  48     END_TYPECHECK(exp(x))
  49
  50     return exp(ex_to_numeric(x)); // -> numeric exp(numeric)
  51 }
  52
  53 static ex exp_eval(ex const & x)
  54 {
  55     // exp(0) -> 1
  56     if (x.is_zero()) {
  57         return _ex1();
  58     }
  59     // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
  60     ex TwoExOverPiI=(_ex2()*x)/(Pi*I);
  61     if (TwoExOverPiI.info(info_flags::integer)) {
  62         numeric z=mod(ex_to_numeric(TwoExOverPiI),_num4());
  63         if (z.is_equal(_num0()))
  64             return _ex1();
  65         if (z.is_equal(_num1()))
  66             return ex(I);
  67         if (z.is_equal(_num2()))
  68             return _ex_1();
  69         if (z.is_equal(_num3()))
  70             return ex(-I);
  71     }
  72     // exp(log(x)) -> x
  73     if (is_ex_the_function(x, log))
  74         return x.op(0);
  75
  76     // exp(float)
  77     if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
  78         return exp_evalf(x);
  79
  80     return exp(x).hold();
  81 }
  82
  83 static ex exp_diff(ex const & x, unsigned diff_param)
  84 {
  85     GINAC_ASSERT(diff_param==0);
  86
  87     // d/dx exp(x) -> exp(x)
  88     return exp(x);
  89 }
  90
  91 REGISTER_FUNCTION(exp, exp_eval, exp_evalf, exp_diff, NULL);
  92
  93 //////////
  94 // natural logarithm
  95 //////////
  96
  97 static ex log_evalf(ex const & x)
  98 {
  99     BEGIN_TYPECHECK
 100         TYPECHECK(x,numeric)
 101     END_TYPECHECK(log(x))
 102
 103     return log(ex_to_numeric(x)); // -> numeric log(numeric)
 104 }
 105
 106 static ex log_eval(ex const & x)
 107 {
 108     if (x.info(info_flags::numeric)) {
 109         if (x.is_equal(_ex1()))  // log(1) -> 0
 110             return _ex0();
 111         if (x.is_equal(_ex_1())) // log(-1) -> I*Pi
 112             return (I*Pi);
 113         if (x.is_equal(I))       // log(I) -> Pi*I/2
 114             return (Pi*I*_num1_2());
 115         if (x.is_equal(-I))      // log(-I) -> -Pi*I/2
 116             return (Pi*I*_num_1_2());
 117         if (x.is_equal(_ex0()))  // log(0) -> infinity
 118             throw(std::domain_error("log_eval(): log(0)"));
 119         // log(float)
 120         if (!x.info(info_flags::crational))
 121             return log_evalf(x);
 122     }
 123     // log(exp(t)) -> t (for real-valued t):
 124     if (is_ex_the_function(x, exp)) {
 125         ex t=x.op(0);
 126         if (t.info(info_flags::real))
 127             return t;
 128     }
 129
 130     return log(x).hold();
 131 }
 132
 133 static ex log_diff(ex const & x, unsigned diff_param)
 134 {
 135     GINAC_ASSERT(diff_param==0);
 136
 137     // d/dx log(x) -> 1/x
 138     return power(x, -1);
 139 }
 140
 141 REGISTER_FUNCTION(log, log_eval, log_evalf, log_diff, NULL);
 142
 143 //////////
 144 // sine (trigonometric function)
 145 //////////
 146
 147 static ex sin_evalf(ex const & x)
 148 {
 149     BEGIN_TYPECHECK
 150        TYPECHECK(x,numeric)
 151     END_TYPECHECK(sin(x))
 152
 153     return sin(ex_to_numeric(x)); // -> numeric sin(numeric)
 154 }
 155
 156 static ex sin_eval(ex const & x)
 157 {
 158     // sin(n*Pi/6) -> { 0 | +/-1/2 | +/-sqrt(3)/2 | +/-1 }
 159     ex SixExOverPi = _ex6()*x/Pi;
 160     if (SixExOverPi.info(info_flags::integer)) {
 161         numeric z = smod(ex_to_numeric(SixExOverPi),_num12());
 162         if (z.is_equal(_num_5()))  // sin(7*Pi/6)  -> -1/2
 163             return _ex_1_2();
 164         if (z.is_equal(_num_4()))  // sin(8*Pi/6)  -> -sqrt(3)/2
 165             return _ex_1_2()*power(_ex3(),_ex1_2());
 166         if (z.is_equal(_num_3()))  // sin(9*Pi/6)  -> -1
 167             return _ex_1();
 168         if (z.is_equal(_num_2()))  // sin(10*Pi/6) -> -sqrt(3)/2
 169             return _ex_1_2()*power(_ex3(),_ex1_2());
 170         if (z.is_equal(_num_1()))  // sin(11*Pi/6) -> -1/2
 171             return _ex_1_2();
 172         if (z.is_equal(_num0()))   // sin(0)       ->  0
 173             return _ex0();
 174         if (z.is_equal(_num1()))   // sin(1*Pi/6)  ->  1/2
 175             return _ex1_2();
 176         if (z.is_equal(_num2()))   // sin(2*Pi/6)  ->  sqrt(3)/2
 177             return _ex1_2()*power(_ex3(),_ex1_2());
 178         if (z.is_equal(_num3()))   // sin(3*Pi/6)  ->  1
 179             return _ex1();
 180         if (z.is_equal(_num4()))   // sin(4*Pi/6)  ->  sqrt(3)/2
 181             return _ex1_2()*power(_ex3(),_ex1_2());
 182         if (z.is_equal(_num5()))   // sin(5*Pi/6)  ->  1/2
 183             return _ex1_2();
 184         if (z.is_equal(_num6()))   // sin(6*Pi/6)  ->  0
 185             return _ex0();
 186     }
 187
 188     if (is_ex_exactly_of_type(x, function)) {
 189         ex t=x.op(0);
 190         // sin(asin(x)) -> x
 191         if (is_ex_the_function(x, asin))
 192             return t;
 193         // sin(acos(x)) -> sqrt(1-x^2)
 194         if (is_ex_the_function(x, acos))
 195             return power(_ex1()-power(t,_ex2()),_ex1_2());
 196         // sin(atan(x)) -> x*(1+x^2)^(-1/2)
 197         if (is_ex_the_function(x, atan))
 198             return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
 199     }
 200
 201     // sin(float) -> float
 202     if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
 203         return sin_evalf(x);
 204
 205     return sin(x).hold();
 206 }
 207
 208 static ex sin_diff(ex const & x, unsigned diff_param)
 209 {
 210     GINAC_ASSERT(diff_param==0);
 211
 212     // d/dx sin(x) -> cos(x)
 213     return cos(x);
 214 }
 215
 216 REGISTER_FUNCTION(sin, sin_eval, sin_evalf, sin_diff, NULL);
 217
 218 //////////
 219 // cosine (trigonometric function)
 220 //////////
 221
 222 static ex cos_evalf(ex const & x)
 223 {
 224     BEGIN_TYPECHECK
 225         TYPECHECK(x,numeric)
 226     END_TYPECHECK(cos(x))
 227
 228     return cos(ex_to_numeric(x)); // -> numeric cos(numeric)
 229 }
 230
 231 static ex cos_eval(ex const & x)
 232 {
 233     // cos(n*Pi/6) -> { 0 | +/-1/2 | +/-sqrt(3)/2 | +/-1 }
 234     ex SixExOverPi = _ex6()*x/Pi;
 235     if (SixExOverPi.info(info_flags::integer)) {
 236         numeric z = smod(ex_to_numeric(SixExOverPi),_num12());
 237         if (z.is_equal(_num_5()))  // cos(7*Pi/6)  -> -sqrt(3)/2
 238             return _ex_1_2()*power(_ex3(),_ex1_2());
 239         if (z.is_equal(_num_4()))  // cos(8*Pi/6)  -> -1/2
 240             return _ex_1_2();
 241         if (z.is_equal(_num_3()))  // cos(9*Pi/6)  ->  0
 242             return _ex0();
 243         if (z.is_equal(_num_2()))  // cos(10*Pi/6) ->  1/2
 244             return _ex1_2();
 245         if (z.is_equal(_num_1()))  // cos(11*Pi/6) ->  sqrt(3)/2
 246             return _ex1_2()*power(_ex3(),_ex1_2());
 247         if (z.is_equal(_num0()))   // cos(0)       ->  1
 248             return _ex1();
 249         if (z.is_equal(_num1()))   // cos(1*Pi/6)  ->  sqrt(3)/2
 250             return _ex1_2()*power(_ex3(),_ex1_2());
 251         if (z.is_equal(_num2()))   // cos(2*Pi/6)  ->  1/2
 252             return _ex1_2();
 253         if (z.is_equal(_num3()))   // cos(3*Pi/6)  ->  0
 254             return _ex0();
 255         if (z.is_equal(_num4()))   // cos(4*Pi/6)  -> -1/2
 256             return _ex_1_2();
 257         if (z.is_equal(_num5()))   // cos(5*Pi/6)  -> -sqrt(3)/2
 258             return _ex_1_2()*power(_ex3(),_ex1_2());
 259         if (z.is_equal(_num6()))   // cos(6*Pi/6)  -> -1
 260             return _ex_1();
 261     }
 262
 263     if (is_ex_exactly_of_type(x, function)) {
 264         ex t=x.op(0);
 265         // cos(acos(x)) -> x
 266         if (is_ex_the_function(x, acos))
 267             return t;
 268         // cos(asin(x)) -> (1-x^2)^(1/2)
 269         if (is_ex_the_function(x, asin))
 270             return power(_ex1()-power(t,_ex2()),_ex1_2());
 271         // cos(atan(x)) -> (1+x^2)^(-1/2)
 272         if (is_ex_the_function(x, atan))
 273             return power(_ex1()+power(t,_ex2()),_ex_1_2());
 274     }
 275
 276     // cos(float) -> float
 277     if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
 278         return cos_evalf(x);
 279
 280     return cos(x).hold();
 281 }
 282
 283 static ex cos_diff(ex const & x, unsigned diff_param)
 284 {
 285     GINAC_ASSERT(diff_param==0);
 286
 287     // d/dx cos(x) -> -sin(x)
 288     return _ex_1()*sin(x);
 289 }
 290
 291 REGISTER_FUNCTION(cos, cos_eval, cos_evalf, cos_diff, NULL);
 292
 293 //////////
 294 // tangent (trigonometric function)
 295 //////////
 296
 297 static ex tan_evalf(ex const & x)
 298 {
 299     BEGIN_TYPECHECK
 300        TYPECHECK(x,numeric)
 301     END_TYPECHECK(tan(x)) // -> numeric tan(numeric)
 302
 303     return tan(ex_to_numeric(x));
 304 }
 305
 306 static ex tan_eval(ex const & x)
 307 {
 308     // tan(n*Pi/6) -> { 0 | +/-sqrt(3) | +/-sqrt(3)/2 }
 309     ex SixExOverPi = _ex6()*x/Pi;
 310     if (SixExOverPi.info(info_flags::integer)) {
 311         numeric z = smod(ex_to_numeric(SixExOverPi),_num6());
 312         if (z.is_equal(_num_2()))  // tan(4*Pi/6) -> -sqrt(3)
 313             return _ex_1()*power(_ex3(),_ex1_2());
 314         if (z.is_equal(_num_1()))  // tan(5*Pi/6) -> -sqrt(3)/3
 315             return _ex_1_3()*power(_ex3(),_ex1_2());
 316         if (z.is_equal(_num0()))   // tan(0)      ->  0
 317             return _ex0();
 318         if (z.is_equal(_num1()))   // tan(1*Pi/6) ->  sqrt(3)/3
 319             return _ex1_3()*power(_ex3(),_ex1_2());
 320         if (z.is_equal(_num2()))   // tan(2*Pi/6) ->  sqrt(3)
 321             return power(_ex3(),_ex1_2());
 322         if (z.is_equal(_num3()))   // tan(3*Pi/6) ->  infinity
 323             throw (std::domain_error("tan_eval(): infinity"));
 324     }
 325
 326     if (is_ex_exactly_of_type(x, function)) {
 327         ex t=x.op(0);
 328         // tan(atan(x)) -> x
 329         if (is_ex_the_function(x, atan))
 330             return t;
 331         // tan(asin(x)) -> x*(1+x^2)^(-1/2)
 332         if (is_ex_the_function(x, asin))
 333             return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
 334         // tan(acos(x)) -> (1-x^2)^(1/2)/x
 335         if (is_ex_the_function(x, acos))
 336             return power(t,_ex_1())*power(_ex1()-power(t,_ex2()),_ex1_2());
 337     }
 338
 339     // tan(float) -> float
 340     if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
 341         return tan_evalf(x);
 342     }
 343
 344     return tan(x).hold();
 345 }
 346
 347 static ex tan_diff(ex const & x, unsigned diff_param)
 348 {
 349     GINAC_ASSERT(diff_param==0);
 350
 351     // d/dx tan(x) -> 1+tan(x)^2;
 352     return (1+power(tan(x),_ex2()));
 353 }
 354
 355 static ex tan_series(ex const & x, symbol const & s, ex const & point, int order)
 356 {
 357     // method:
 358     // Taylor series where there is no pole falls back to tan_diff.
 359     // On a pole simply expand sin(x)/cos(x).
 360     ex xpoint = x.subs(s==point);
 361     if (!(2*xpoint/Pi).info(info_flags::odd))
 362         throw do_taylor();  // caught by function::series()
 363     // if we got here we have to care for a simple pole
 364     return (sin(x)/cos(x)).series(s, point, order+2);
 365 }
 366
 367 REGISTER_FUNCTION(tan, tan_eval, tan_evalf, tan_diff, tan_series);
 368
 369 //////////
 370 // inverse sine (arc sine)
 371 //////////
 372
 373 static ex asin_evalf(ex const & x)
 374 {
 375     BEGIN_TYPECHECK
 376        TYPECHECK(x,numeric)
 377     END_TYPECHECK(asin(x))
 378
 379     return asin(ex_to_numeric(x)); // -> numeric asin(numeric)
 380 }
 381
 382 static ex asin_eval(ex const & x)
 383 {
 384     if (x.info(info_flags::numeric)) {
 385         // asin(0) -> 0
 386         if (x.is_zero())
 387             return x;
 388         // asin(1/2) -> Pi/6
 389         if (x.is_equal(_ex1_2()))
 390             return numeric(1,6)*Pi;
 391         // asin(1) -> Pi/2
 392         if (x.is_equal(_ex1()))
 393             return _num1_2()*Pi;
 394         // asin(-1/2) -> -Pi/6
 395         if (x.is_equal(_ex_1_2()))
 396             return numeric(-1,6)*Pi;
 397         // asin(-1) -> -Pi/2
 398         if (x.is_equal(_ex_1()))
 399             return _num_1_2()*Pi;
 400         // asin(float) -> float
 401         if (!x.info(info_flags::crational))
 402             return asin_evalf(x);
 403     }
 404
 405     return asin(x).hold();
 406 }
 407
 408 static ex asin_diff(ex const & x, unsigned diff_param)
 409 {
 410     GINAC_ASSERT(diff_param==0);
 411
 412     // d/dx asin(x) -> 1/sqrt(1-x^2)
 413     return power(1-power(x,_ex2()),_ex_1_2());
 414 }
 415
 416 REGISTER_FUNCTION(asin, asin_eval, asin_evalf, asin_diff, NULL);
 417
 418 //////////
 419 // inverse cosine (arc cosine)
 420 //////////
 421
 422 static ex acos_evalf(ex const & x)
 423 {
 424     BEGIN_TYPECHECK
 425        TYPECHECK(x,numeric)
 426     END_TYPECHECK(acos(x))
 427
 428     return acos(ex_to_numeric(x)); // -> numeric acos(numeric)
 429 }
 430
 431 static ex acos_eval(ex const & x)
 432 {
 433     if (x.info(info_flags::numeric)) {
 434         // acos(1) -> 0
 435         if (x.is_equal(_ex1()))
 436             return _ex0();
 437         // acos(1/2) -> Pi/3
 438         if (x.is_equal(_ex1_2()))
 439             return numeric(1,3)*Pi;
 440         // acos(0) -> Pi/2
 441         if (x.is_zero())
 442             return numeric(1,2)*Pi;
 443         // acos(-1/2) -> 2/3*Pi
 444         if (x.is_equal(_ex_1_2()))
 445             return numeric(2,3)*Pi;
 446         // acos(-1) -> Pi
 447         if (x.is_equal(_ex_1()))
 448             return Pi;
 449         // acos(float) -> float
 450         if (!x.info(info_flags::crational))
 451             return acos_evalf(x);
 452     }
 453
 454     return acos(x).hold();
 455 }
 456
 457 static ex acos_diff(ex const & x, unsigned diff_param)
 458 {
 459     GINAC_ASSERT(diff_param==0);
 460
 461     // d/dx acos(x) -> -1/sqrt(1-x^2)
 462     return _ex_1()*power(1-power(x,_ex2()),_ex_1_2());
 463 }
 464
 465 REGISTER_FUNCTION(acos, acos_eval, acos_evalf, acos_diff, NULL);
 466
 467 //////////
 468 // inverse tangent (arc tangent)
 469 //////////
 470
 471 static ex atan_evalf(ex const & x)
 472 {
 473     BEGIN_TYPECHECK
 474         TYPECHECK(x,numeric)
 475     END_TYPECHECK(atan(x))
 476
 477     return atan(ex_to_numeric(x)); // -> numeric atan(numeric)
 478 }
 479
 480 static ex atan_eval(ex const & x)
 481 {
 482     if (x.info(info_flags::numeric)) {
 483         // atan(0) -> 0
 484         if (x.is_equal(_ex0()))
 485             return _ex0();
 486         // atan(float) -> float
 487         if (!x.info(info_flags::crational))
 488             return atan_evalf(x);
 489     }
 490
 491     return atan(x).hold();
 492 }
 493
 494 static ex atan_diff(ex const & x, unsigned diff_param)
 495 {
 496     GINAC_ASSERT(diff_param==0);
 497
 498     return power(1+x*x, -1);
 499 }
 500
 501 REGISTER_FUNCTION(atan, atan_eval, atan_evalf, atan_diff, NULL);
 502
 503 //////////
 504 // inverse tangent (atan2(y,x))
 505 //////////
 506
 507 static ex atan2_evalf(ex const & y, ex const & x)
 508 {
 509     BEGIN_TYPECHECK
 510         TYPECHECK(y,numeric)
 511         TYPECHECK(x,numeric)
 512     END_TYPECHECK(atan2(y,x))
 513
 514     return atan(ex_to_numeric(y),ex_to_numeric(x)); // -> numeric atan(numeric)
 515 }
 516
 517 static ex atan2_eval(ex const & y, ex const & x)
 518 {
 519     if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
 520         x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
 521         return atan2_evalf(y,x);
 522     }
 523
 524     return atan2(y,x).hold();
 525 }
 526
 527 static ex atan2_diff(ex const & y, ex const & x, unsigned diff_param)
 528 {
 529     GINAC_ASSERT(diff_param<2);
 530
 531     if (diff_param==0) {
 532         // d/dy atan(y,x)
 533         return x*pow(pow(x,2)+pow(y,2),-1);
 534     }
 535     // d/dx atan(y,x)
 536     return -y*pow(pow(x,2)+pow(y,2),-1);
 537 }
 538
 539 REGISTER_FUNCTION(atan2, atan2_eval, atan2_evalf, atan2_diff, NULL);
 540
 541 //////////
 542 // hyperbolic sine (trigonometric function)
 543 //////////
 544
 545 static ex sinh_evalf(ex const & x)
 546 {
 547     BEGIN_TYPECHECK
 548        TYPECHECK(x,numeric)
 549     END_TYPECHECK(sinh(x))
 550
 551     return sinh(ex_to_numeric(x)); // -> numeric sinh(numeric)
 552 }
 553
 554 static ex sinh_eval(ex const & x)
 555 {
 556     if (x.info(info_flags::numeric)) {
 557         if (x.is_zero())  // sinh(0) -> 0
 558             return _ex0();
 559         if (!x.info(info_flags::crational))  // sinh(float) -> float
 560             return sinh_evalf(x);
 561     }
 562
 563     if ((x/Pi).info(info_flags::numeric) &&
 564         ex_to_numeric(x/Pi).real().is_zero())  // sinh(I*x) -> I*sin(x)
 565         return I*sin(x/I);
 566
 567     if (is_ex_exactly_of_type(x, function)) {
 568         ex t=x.op(0);
 569         // sinh(asinh(x)) -> x
 570         if (is_ex_the_function(x, asinh))
 571             return t;
 572         // sinh(acosh(x)) -> (x-1)^(1/2) * (x+1)^(1/2)
 573         if (is_ex_the_function(x, acosh))
 574             return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2());
 575         // sinh(atanh(x)) -> x*(1-x^2)^(-1/2)
 576         if (is_ex_the_function(x, atanh))
 577             return t*power(_ex1()-power(t,_ex2()),_ex_1_2());
 578     }
 579
 580     return sinh(x).hold();
 581 }
 582
 583 static ex sinh_diff(ex const & x, unsigned diff_param)
 584 {
 585     GINAC_ASSERT(diff_param==0);
 586
 587     // d/dx sinh(x) -> cosh(x)
 588     return cosh(x);
 589 }
 590
 591 REGISTER_FUNCTION(sinh, sinh_eval, sinh_evalf, sinh_diff, NULL);
 592
 593 //////////
 594 // hyperbolic cosine (trigonometric function)
 595 //////////
 596
 597 static ex cosh_evalf(ex const & x)
 598 {
 599     BEGIN_TYPECHECK
 600        TYPECHECK(x,numeric)
 601     END_TYPECHECK(cosh(x))
 602
 603     return cosh(ex_to_numeric(x)); // -> numeric cosh(numeric)
 604 }
 605
 606 static ex cosh_eval(ex const & x)
 607 {
 608     if (x.info(info_flags::numeric)) {
 609         if (x.is_zero())  // cosh(0) -> 1
 610             return _ex1();
 611         if (!x.info(info_flags::crational))  // cosh(float) -> float
 612             return cosh_evalf(x);
 613     }
 614
 615     if ((x/Pi).info(info_flags::numeric) &&
 616         ex_to_numeric(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
 617         return cos(x/I);
 618
 619     if (is_ex_exactly_of_type(x, function)) {
 620         ex t=x.op(0);
 621         // cosh(acosh(x)) -> x
 622         if (is_ex_the_function(x, acosh))
 623             return t;
 624         // cosh(asinh(x)) -> (1+x^2)^(1/2)
 625         if (is_ex_the_function(x, asinh))
 626             return power(_ex1()+power(t,_ex2()),_ex1_2());
 627         // cosh(atanh(x)) -> (1-x^2)^(-1/2)
 628         if (is_ex_the_function(x, atanh))
 629             return power(_ex1()-power(t,_ex2()),_ex_1_2());
 630     }
 631
 632     return cosh(x).hold();
 633 }
 634
 635 static ex cosh_diff(ex const & x, unsigned diff_param)
 636 {
 637     GINAC_ASSERT(diff_param==0);
 638
 639     // d/dx cosh(x) -> sinh(x)
 640     return sinh(x);
 641 }
 642
 643 REGISTER_FUNCTION(cosh, cosh_eval, cosh_evalf, cosh_diff, NULL);
 644
 645 //////////
 646 // hyperbolic tangent (trigonometric function)
 647 //////////
 648
 649 static ex tanh_evalf(ex const & x)
 650 {
 651     BEGIN_TYPECHECK
 652        TYPECHECK(x,numeric)
 653     END_TYPECHECK(tanh(x))
 654
 655     return tanh(ex_to_numeric(x)); // -> numeric tanh(numeric)
 656 }
 657
 658 static ex tanh_eval(ex const & x)
 659 {
 660     if (x.info(info_flags::numeric)) {
 661         if (x.is_zero())  // tanh(0) -> 0
 662             return _ex0();
 663         if (!x.info(info_flags::crational))  // tanh(float) -> float
 664             return tanh_evalf(x);
 665     }
 666
 667     if ((x/Pi).info(info_flags::numeric) &&
 668         ex_to_numeric(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
 669         return I*tan(x/I);
 670
 671     if (is_ex_exactly_of_type(x, function)) {
 672         ex t=x.op(0);
 673         // tanh(atanh(x)) -> x
 674         if (is_ex_the_function(x, atanh))
 675             return t;
 676         // tanh(asinh(x)) -> x*(1+x^2)^(-1/2)
 677         if (is_ex_the_function(x, asinh))
 678             return t*power(_ex1()+power(t,_ex2()),_ex_1_2());
 679         // tanh(acosh(x)) -> (x-1)^(1/2)*(x+1)^(1/2)/x
 680         if (is_ex_the_function(x, acosh))
 681             return power(t-_ex1(),_ex1_2())*power(t+_ex1(),_ex1_2())*power(t,_ex_1());
 682     }
 683
 684     return tanh(x).hold();
 685 }
 686
 687 static ex tanh_diff(ex const & x, unsigned diff_param)
 688 {
 689     GINAC_ASSERT(diff_param==0);
 690
 691     // d/dx tanh(x) -> 1-tanh(x)^2
 692     return _ex1()-power(tanh(x),_ex2());
 693 }
 694
 695 static ex tanh_series(ex const & x, symbol const & s, ex const & point, int order)
 696 {
 697     // method:
 698     // Taylor series where there is no pole falls back to tanh_diff.
 699     // On a pole simply expand sinh(x)/cosh(x).
 700     ex xpoint = x.subs(s==point);
 701     if (!(2*I*xpoint/Pi).info(info_flags::odd))
 702         throw do_taylor();  // caught by function::series()
 703     // if we got here we have to care for a simple pole
 704     return (sinh(x)/cosh(x)).series(s, point, order+2);
 705 }
 706
 707 REGISTER_FUNCTION(tanh, tanh_eval, tanh_evalf, tanh_diff, tanh_series);
 708
 709 //////////
 710 // inverse hyperbolic sine (trigonometric function)
 711 //////////
 712
 713 static ex asinh_evalf(ex const & x)
 714 {
 715     BEGIN_TYPECHECK
 716        TYPECHECK(x,numeric)
 717     END_TYPECHECK(asinh(x))
 718
 719     return asinh(ex_to_numeric(x)); // -> numeric asinh(numeric)
 720 }
 721
 722 static ex asinh_eval(ex const & x)
 723 {
 724     if (x.info(info_flags::numeric)) {
 725         // asinh(0) -> 0
 726         if (x.is_zero())
 727             return _ex0();
 728         // asinh(float) -> float
 729         if (!x.info(info_flags::crational))
 730             return asinh_evalf(x);
 731     }
 732
 733     return asinh(x).hold();
 734 }
 735
 736 static ex asinh_diff(ex const & x, unsigned diff_param)
 737 {
 738     GINAC_ASSERT(diff_param==0);
 739
 740     // d/dx asinh(x) -> 1/sqrt(1+x^2)
 741     return power(1+power(x,_ex2()),_ex_1_2());
 742 }
 743
 744 REGISTER_FUNCTION(asinh, asinh_eval, asinh_evalf, asinh_diff, NULL);
 745
 746 //////////
 747 // inverse hyperbolic cosine (trigonometric function)
 748 //////////
 749
 750 static ex acosh_evalf(ex const & x)
 751 {
 752     BEGIN_TYPECHECK
 753        TYPECHECK(x,numeric)
 754     END_TYPECHECK(acosh(x))
 755
 756     return acosh(ex_to_numeric(x)); // -> numeric acosh(numeric)
 757 }
 758
 759 static ex acosh_eval(ex const & x)
 760 {
 761     if (x.info(info_flags::numeric)) {
 762         // acosh(0) -> Pi*I/2
 763         if (x.is_zero())
 764             return Pi*I*numeric(1,2);
 765         // acosh(1) -> 0
 766         if (x.is_equal(_ex1()))
 767             return _ex0();
 768         // acosh(-1) -> Pi*I
 769         if (x.is_equal(_ex_1()))
 770             return Pi*I;
 771         // acosh(float) -> float
 772         if (!x.info(info_flags::crational))
 773             return acosh_evalf(x);
 774     }
 775
 776     return acosh(x).hold();
 777 }
 778
 779 static ex acosh_diff(ex const & x, unsigned diff_param)
 780 {
 781     GINAC_ASSERT(diff_param==0);
 782
 783     // d/dx acosh(x) -> 1/(sqrt(x-1)*sqrt(x+1))
 784     return power(x+_ex_1(),_ex_1_2())*power(x+_ex1(),_ex_1_2());
 785 }
 786
 787 REGISTER_FUNCTION(acosh, acosh_eval, acosh_evalf, acosh_diff, NULL);
 788
 789 //////////
 790 // inverse hyperbolic tangent (trigonometric function)
 791 //////////
 792
 793 static ex atanh_evalf(ex const & x)
 794 {
 795     BEGIN_TYPECHECK
 796        TYPECHECK(x,numeric)
 797     END_TYPECHECK(atanh(x))
 798
 799     return atanh(ex_to_numeric(x)); // -> numeric atanh(numeric)
 800 }
 801
 802 static ex atanh_eval(ex const & x)
 803 {
 804     if (x.info(info_flags::numeric)) {
 805         // atanh(0) -> 0
 806         if (x.is_zero())
 807             return _ex0();
 808         // atanh({+|-}1) -> throw
 809         if (x.is_equal(_ex1()) || x.is_equal(_ex1()))
 810             throw (std::domain_error("atanh_eval(): infinity"));
 811         // atanh(float) -> float
 812         if (!x.info(info_flags::crational))
 813             return atanh_evalf(x);
 814     }
 815
 816     return atanh(x).hold();
 817 }
 818
 819 static ex atanh_diff(ex const & x, unsigned diff_param)
 820 {
 821     GINAC_ASSERT(diff_param==0);
 822
 823     // d/dx atanh(x) -> 1/(1-x^2)
 824     return power(_ex1()-power(x,_ex2()),_ex_1());
 825 }
 826
 827 REGISTER_FUNCTION(atanh, atanh_eval, atanh_evalf, atanh_diff, NULL);
 828
 829 #ifndef NO_GINAC_NAMESPACE
 830 } // namespace GiNaC
 831 #endif // ndef NO_GINAC_NAMESPACE