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