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