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