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