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