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