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