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