1 /** @file inifcns_trans.cpp 2 * 3 * Implementation of trigonometric functions. */ 4 5 /* 6 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany 7 * 8 * This program is free software; you can redistribute it and/or modify 9 * it under the terms of the GNU General Public License as published by 10 * the Free Software Foundation; either version 2 of the License, or 11 * (at your option) any later version. 12 * 13 * This program is distributed in the hope that it will be useful, 14 * but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 * GNU General Public License for more details. 17 * 18 * You should have received a copy of the GNU General Public License 19 * along with this program; if not, write to the Free Software 20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 21 */ 22 23 #include "inifcns.h" 24 #include "ex.h" 25 #include "ex_utils.h" 26 #include "constant.h" 27 #include "infinity.h" 28 #include "symbol.h" 29 #include "mul.h" 30 #include "power.h" 31 #include "operators.h" 32 #include "pseries.h" 33 #include "utils.h" 34 #include "add.h" 35 36 namespace GiNaC { 37 38 39 /* In Sage all the arc trig functions are printed with "arc" instead 40 of "a" at the beginning. This is for consistency with other 41 computer algebra systems. These print methods are registered 42 below with each of the corresponding inverse trig function. */ 43 44 static bool has_pi(const ex & the_ex) { 45 if (is_exactly_a<constant>(the_ex) and ex_to<constant>(the_ex) == Pi) 46 return true; 47 for (size_t i=0; i<the_ex.nops(); ++i) 48 if (has_pi(the_ex.op(i))) 49 return true; 50 return false; 51 } 52 53 // helper function: returns whether the expression is a multiple of I 54 static bool is_multiple_of_I(const ex & the_ex) 55 { 56 57 if (is_exactly_a<numeric>(the_ex) 58 and the_ex.real_part().is_zero()) 59 return true; 60 61 if (is_exactly_a<mul>(the_ex)) { 62 for (size_t i=0; i < the_ex.nops(); ++i) 63 if (is_multiple_of_I(the_ex.op(i))) 64 return true; 65 } 66 67 if (is_exactly_a<add>(the_ex)) { 68 for (size_t i=0; i < the_ex.nops(); ++i) 69 if (!is_multiple_of_I(the_ex.op(i))) 70 return false; 71 return true; 72 } 73 return false; 74 }; 75 76 77 ////////// 78 // sine (trigonometric function) 79 ////////// 80 81 static ex sin_eval(const ex & x) 82 { 83 // sin(oo) -> error 84 // This should be before the tests below, since multiplying infinity 85 // with other values raises runtime_errors 86 if (x.info(info_flags::infinity)) { 87 throw (std::runtime_error("sin_eval(): sin(infinity) encountered")); 88 } 89 90 if (is_exactly_a<numeric>(x) and ex_to<numeric>(x).is_zero()) 91 return _ex0; 92 93 // sin() is odd 94 if (x.info(info_flags::negative)) 95 return -sin(-x); 96 97 ex x_red; 98 if (has_pi(x)) { 99 100 ex coef_pi = x.coeff(Pi,_ex1).expand(); 101 ex rem = _ex0; 102 if (is_exactly_a<add>(coef_pi)) { 103 for (size_t i=0; i < coef_pi.nops(); i++) { 104 if ((coef_pi.op(i) / _ex2).is_integer()) 105 rem += Pi * coef_pi.op(i); 106 } 107 } 108 else if ((coef_pi / _ex2).is_integer()) 109 rem = Pi * coef_pi; 110 x_red = (x - rem).expand(); 111 112 // sin(n/d*Pi) -> { all known radicals with nesting depth 2 } 113 const ex SixtyExOverPi = _ex60*x_red/Pi; 114 ex sign = _ex1; 115 if (is_exactly_a<numeric>(SixtyExOverPi) 116 and SixtyExOverPi.is_integer()) { 117 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p); 118 if (z>=*_num60_p) { 119 // wrap to interval [0, Pi) 120 z -= *_num60_p; 121 sign = _ex_1; 122 } 123 if (z>*_num30_p) { 124 // wrap to interval [0, Pi/2) 125 z = *_num60_p-z; 126 } 127 // Not included were n*Pi/15 which has 3-level-deep roots 128 if (z.is_equal(*_num0_p)) // sin(0) -> 0 129 return _ex0; 130 if (z.is_equal(*_num2_p)) // sin(Pi/30) -> -1/8*sqrt(5) + 1/2*sqrt(-3/8*sqrt(5) + 15/8) - 1/8 131 return sign*(_ex_1*(sqrt(_ex5)+_ex1)/_ex8 + 132 _ex1_4*sqrt(_ex1_2*(_ex15-_ex3*sqrt(_ex5)))); 133 if (z.is_equal(*_num5_p)) // sin(Pi/12) -> (sqrt(6)-sqrt(2))/4 134 return sign*_ex1_4*(sqrt(_ex6)-sqrt(_ex2)); 135 if (z.is_equal(*_num6_p)) // sin(Pi/10) -> sqrt(5)/4-1/4 136 return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4); 137 if (z.is_equal(*_num10_p)) // sin(Pi/6) -> 1/2 138 return sign*_ex1_2; 139 if (z.is_equal(*_num12_p)) // sin(Pi/5) -> 1/4*sqrt(10-2*sqrt(5)) 140 return sign*_ex1_4*sqrt(_ex10-_ex2*sqrt(_ex5)); 141 if (z.is_equal(*_num14_p)) // sin(7*Pi/30) -> -1/8*sqrt(5) + 1/2*sqrt(3/8*sqrt(5) + 15/8) + 1/8 142 return sign*((_ex1-sqrt(_ex5))/_ex8 + 143 _ex1_4*sqrt(_ex1_2*(_ex15+_ex3*sqrt(_ex5)))); 144 if (z.is_equal(*_num15_p)) // sin(Pi/4) -> sqrt(2)/2 145 return sign*_ex1_2*sqrt(_ex2); 146 if (z.is_equal(*_num18_p)) // sin(3/10*Pi) -> sqrt(5)/4+1/4 147 return sign*(_ex1_4*sqrt(_ex5)+_ex1_4); 148 if (z.is_equal(*_num20_p)) // sin(Pi/3) -> sqrt(3)/2 149 return sign*_ex1_2*sqrt(_ex3); 150 if (z.is_equal(*_num22_p)) // sin(11*Pi/30) -> 1/8*sqrt(5) + 1/2*sqrt(-3/8*sqrt(5) + 15/8) + 1/8 151 return sign*((sqrt(_ex5)+_ex1)/_ex8 + 152 _ex1_4*sqrt(_ex1_2*(_ex15-_ex3*sqrt(_ex5)))); 153 if (z.is_equal(*_num24_p)) // sin(2*Pi/5) -> 1/4*sqrt(10+2*sqrt(5)) 154 return sign*_ex1_4*sqrt(_ex10+_ex2*sqrt(_ex5)); 155 if (z.is_equal(*_num25_p)) // sin(5/12*Pi) -> (sqrt(6)+sqrt(2))/4 156 return sign*_ex1_4*(sqrt(_ex6)+sqrt(_ex2)); 157 if (z.is_equal(*_num26_p)) // sin(13*Pi/30) -> 1/8*sqrt(5) + 1/2*sqrt(3/8*sqrt(5) + 15/8) - 1/8 158 return sign*((sqrt(_ex5)-_ex1)/_ex8 + 159 _ex1_4*sqrt(_ex1_2*(_ex15+_ex3*sqrt(_ex5)))); 160 if (z.is_equal(*_num30_p)) // sin(Pi/2) -> 1 161 return sign; 162 } 163 164 const ex TwentyforExOverPi = _ex24*x_red/Pi; 165 sign = _ex1; 166 if (is_exactly_a<numeric>(TwentyforExOverPi) 167 and TwentyforExOverPi.is_integer()) { 168 numeric z = mod(ex_to<numeric>(TwentyforExOverPi),*_num48_p); 169 if (z>=*_num24_p) { 170 // wrap to interval [0, Pi) 171 z -= *_num24_p; 172 sign = _ex_1; 173 } 174 if (z>*_num12_p) { 175 // wrap to interval [0, Pi/2) 176 z = *_num24_p-z; 177 } 178 if (z.is_equal(*_num1_p)) // sin(Pi/24) -> 1/4*sqrt(8-2*sqrt(6)-2*sqrt(2)) 179 return sign*(_ex1_4*sqrt(_ex8 - _ex2*sqrt(_ex6) - _ex2*sqrt(_ex2))); 180 if (z.is_equal(*_num3_p)) // sin(Pi/8) -> 1/2*sqrt(-sqrt(2) + 2) 181 return sign*(_ex1_2*sqrt(_ex2-sqrt(_ex2))); 182 if (z.is_equal(*_num5_p)) // sin(5*Pi/24) -> 1/4*sqrt(8-2*sqrt(6)+2*sqrt(2)) 183 return sign*(_ex1_4*sqrt(_ex8 - _ex2*sqrt(_ex6) + _ex2*sqrt(_ex2))); 184 if (z.is_equal(*_num7_p)) // sin(7*Pi/24) -> 1/4*sqrt(8+2*sqrt(6)-2*sqrt(2)) 185 return sign *(_ex1_4*sqrt(_ex8 + _ex2*sqrt(_ex6) - _ex2*sqrt(_ex2))); 186 if (z.is_equal(*_num9_p)) // sin(3*Pi/8) -> 1/2*sqrt(-sqrt(2) + 2) 187 return sign*(_ex1_2*sqrt(_ex2+sqrt(_ex2))); 188 if (z.is_equal(*_num11_p)) // sin(11*Pi/24) -> 1/4*sqrt(8+2*sqrt(6)+2*sqrt(2)) 189 return sign *(_ex1_4*sqrt(_ex8 + _ex2*sqrt(_ex6) + _ex2*sqrt(_ex2))); 190 } 191 192 // Reflection at Pi/2 193 const ex ExOverPi = x_red/Pi; 194 if (ExOverPi.is_integer()) 195 return _ex0; 196 if (is_exactly_a<numeric>(ExOverPi)) { 197 const numeric& c = ex_to<numeric>(ExOverPi); 198 if (c.is_rational()) { 199 numeric den = c.denom(); 200 numeric num = c.numer().mod(den * *_num2_p); 201 numeric fac = *_num1_p; 202 if (num > den) { 203 num -= den; 204 fac = *_num_1_p; 205 } 206 if (num*(*_num2_p) > den) 207 return fac*sin_eval((den-num)*Pi/den); 208 return fac*sin((num*Pi)/den).hold(); 209 } 210 } 211 } 212 else 213 x_red = x; 214 215 // simplify sin(I*x) --> I*sinh(x) 216 if (is_multiple_of_I(x_red.expand())) 217 return I*sinh(x_red/I); 218 219 if (is_exactly_a<function>(x_red)) { 220 const ex &t = x_red.op(0); 221 222 // sin(asin(x)) -> x 223 if (is_ex_the_function(x_red, asin)) 224 return t; 225 226 // sin(acos(x)) -> sqrt(1-x^2) 227 if (is_ex_the_function(x_red, acos)) 228 return sqrt(_ex1-power(t,_ex2)); 229 230 // sin(atan(x)) -> x/sqrt(1+x^2) 231 if (is_ex_the_function(x_red, atan)) 232 return t*power(_ex1+power(t,_ex2),_ex_1_2); 233 234 // sin(acot(x)) -> 1/sqrt(1+x^2) 235 if (is_ex_the_function(x_red, acot)) 236 return power(_ex1+power(t,_ex2),_ex_1_2); 237 238 // sin(acosec(x)) -> 1/x 239 if (is_ex_the_function(x_red, acsc)) 240 return _ex1/t; 241 242 // sin(asec(x)) -> sqrt(x^2-1)/x 243 if (is_ex_the_function(x_red, asec)) 244 return sqrt(power(t,_ex2)-_ex1)/t; 245 246 // sin(atan2(y,x)) -> y/sqrt(x^2+y^2) 247 if (is_ex_the_function(x_red, atan2)) { 248 const ex &t1 = x_red.op(1); 249 return t*power(power(t,_ex2)+power(t1,_ex2),_ex_1_2); 250 } 251 } 252 253 // sin(float) -> float 254 if (is_exactly_a<numeric>(x_red) && x_red.info(info_flags::inexact)) 255 return sin(ex_to<numeric>(x_red)); 256 257 return sin(x_red).hold(); 258 } 259 260 static ex sin_deriv(const ex & x, unsigned deriv_param) 261 { 262 GINAC_ASSERT(deriv_param==0); 263 264 // d/dx sin(x) -> cos(x) 265 return cos(x); 266 } 267 268 static ex sin_real_part(const ex & x) 269 { 270 return cosh(GiNaC::imag_part(x))*sin(GiNaC::real_part(x)); 271 } 272 273 static ex sin_imag_part(const ex & x) 274 { 275 return sinh(GiNaC::imag_part(x))*cos(GiNaC::real_part(x)); 276 } 277 278 static ex sin_conjugate(const ex & x) 279 { 280 // conjugate(sin(x))==sin(conjugate(x)) 281 return sin(x.conjugate()); 282 } 283 284 REGISTER_FUNCTION(sin, eval_func(sin_eval). 285 derivative_func(sin_deriv). 286 real_part_func(sin_real_part). 287 imag_part_func(sin_imag_part). 288 conjugate_func(sin_conjugate). 289 latex_name("\\sin")); 290 291 ////////// 292 // cosine (trigonometric function) 293 ////////// 294 295 static ex cos_eval(const ex & x) 296 { 297 // cos(oo) -> error 298 // This should be before the tests below, since multiplying infinity 299 // with other values raises runtime_errors 300 if (x.info(info_flags::infinity)) { 301 throw (std::runtime_error("cos_eval(): cos(infinity) encountered")); 302 } 303 304 // cos() is even 305 if (x.info(info_flags::negative)) 306 return cos(-x); 307 308 if (is_exactly_a<numeric>(x) and ex_to<numeric>(x).is_zero()) 309 return _ex1; 310 311 ex x_red; 312 if (has_pi(x)) { 313 314 ex coef_pi = x.coeff(Pi,_ex1).expand(); 315 ex rem = _ex0; 316 if (is_exactly_a<add>(coef_pi)) { 317 for (size_t i=0; i < coef_pi.nops(); i++) { 318 if ((coef_pi.op(i) / _ex2).is_integer()) 319 rem += Pi * coef_pi.op(i); 320 } 321 } 322 else if ((coef_pi / _ex2).is_integer()) 323 rem = Pi * coef_pi; 324 x_red = (x - rem).expand(); 325 326 // cos(n/d*Pi) -> { all known radicals with nesting depth 2 } 327 const ex SixtyExOverPi = _ex60*x_red/Pi; 328 ex sign = _ex1; 329 if (is_exactly_a<numeric>(SixtyExOverPi) 330 and SixtyExOverPi.is_integer()) { 331 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num120_p); 332 if (z>=*_num60_p) { 333 // wrap to interval [0, Pi) 334 z = *_num120_p-z; 335 } 336 if (z>=*_num30_p) { 337 // wrap to interval [0, Pi/2) 338 z = *_num60_p-z; 339 sign = _ex_1; 340 } 341 if (z.is_equal(*_num0_p)) // cos(0) -> 1 342 return sign; 343 if (z.is_equal(*_num4_p)) // cos(Pi/15) -> 1/8*sqrt(5) + 1/2*sqrt(3/8*sqrt(5) + 15/8) - 1/8 344 return sign*((sqrt(_ex5)-_ex1)/_ex8 + _ex1_4*sqrt((_ex15+_ex3*sqrt(_ex5))/_ex2)); 345 if (z.is_equal(*_num5_p)) // cos(Pi/12) -> (sqrt(6)+sqrt(2))/4 346 return sign*_ex1_4*(sqrt(_ex6)+sqrt(_ex2)); 347 if (z.is_equal(*_num6_p)) // cos(Pi/10) -> 1/4*sqrt(2*sqrt(5) + 10) 348 return sign*_ex1_4*sqrt(_ex10+_ex2*sqrt(_ex5)); 349 if (z.is_equal(*_num8_p)) // cos(2*Pi/15) -> 1/8*sqrt(5) + 1/2*sqrt(-3/8*sqrt(5) + 15/8) + 1/8 350 return sign*((sqrt(_ex5)+_ex1)/_ex8 + _ex1_4*sqrt((_ex15-_ex3*sqrt(_ex5))/_ex2)); 351 if (z.is_equal(*_num10_p)) // cos(Pi/6) -> sqrt(3)/2 352 return sign*_ex1_2*sqrt(_ex3); 353 if (z.is_equal(*_num12_p)) // cos(Pi/5) -> sqrt(5)/4+1/4 354 return sign*(_ex1_4*sqrt(_ex5)+_ex1_4); 355 if (z.is_equal(*_num15_p)) // cos(Pi/4) -> sqrt(2)/2 356 return sign*_ex1_2*sqrt(_ex2); 357 if (z.is_equal(*_num16_p)) // cos(4*Pi/15) -> -1/8*sqrt(5) + 1/2*sqrt(3/8*sqrt(5) + 15/8) + 1/8 358 return sign*((_ex1-sqrt(_ex5))/_ex8 + _ex1_4*sqrt((_ex15+_ex3*sqrt(_ex5))/_ex2)); 359 if (z.is_equal(*_num18_p)) // cos(3*Pi/10) -> 1/4*sqrt(-2*sqrt(5) + 10) 360 return sign*_ex1_4*sqrt(_ex10-_ex2*sqrt(_ex5)); 361 if (z.is_equal(*_num20_p)) // cos(Pi/3) -> 1/2 362 return sign*_ex1_2; 363 if (z.is_equal(*_num24_p)) // cos(2/5*Pi) -> sqrt(5)/4-1/4x 364 return sign*(_ex1_4*sqrt(_ex5)+_ex_1_4); 365 if (z.is_equal(*_num25_p)) // cos(5/12*Pi) -> (sqrt(6)-sqrt(2))/4 366 return sign*_ex1_4*(sqrt(_ex6)-sqrt(_ex2)); 367 if (z.is_equal(*_num28_p)) // cos(7*Pi/15) -> -1/8*sqrt(5) + 1/2*sqrt(-3/8*sqrt(5) + 15/8) - 1/8 368 return sign*(_ex_1*(_ex1+sqrt(_ex5))/_ex8 + _ex1_4*sqrt((_ex15-_ex3*sqrt(_ex5))/_ex2)); 369 if (z.is_equal(*_num30_p)) // cos(Pi/2) -> 0 370 return _ex0; 371 } 372 373 const ex TwentyforExOverPi = _ex24*x_red/Pi; 374 sign = _ex1; 375 if (is_exactly_a<numeric>(TwentyforExOverPi) 376 and TwentyforExOverPi.is_integer()) { 377 numeric z = mod(ex_to<numeric>(TwentyforExOverPi),*_num48_p); 378 if (z>=*_num24_p) { 379 // wrap to interval [0, Pi) 380 z -= *_num48_p; 381 } 382 if (z>*_num12_p) { 383 // wrap to interval [0, Pi/2) 384 z = *_num24_p-z; 385 sign = _ex_1; 386 } 387 if (z.is_equal(*_num1_p)) // cos(Pi/24) -> 1/4*sqrt(8+2*sqrt(6)+2*sqrt(2)) 388 return sign*(_ex1_4*sqrt(_ex8 + _ex2*sqrt(_ex6) + _ex2*sqrt(_ex2))); 389 if (z.is_equal(*_num3_p)) // cos(Pi/8) -> 1/2*sqrt(sqrt(2) + 2) 390 return sign*(_ex1_2*sqrt(_ex2+sqrt(_ex2))); 391 if (z.is_equal(*_num5_p)) // cos(5*Pi/24) -> 1/4*sqrt(8+2*sqrt(6)-2*sqrt(2)) 392 return sign*(_ex1_4*sqrt(_ex8 + _ex2*sqrt(_ex6) - _ex2*sqrt(_ex2))); 393 if (z.is_equal(*_num7_p)) // cos(7*Pi/24) -> 1/4*sqrt(8-2*sqrt(6)+2*sqrt(2)) 394 return sign*(_ex1_4*sqrt(_ex8 - _ex2*sqrt(_ex6) + _ex2*sqrt(_ex2))); 395 if (z.is_equal(*_num9_p)) // cos(3*Pi/8) -> 1/2*sqrt(-sqrt(2) + 2) 396 return sign*(_ex1_2*sqrt(_ex2-sqrt(_ex2))); 397 if (z.is_equal(*_num11_p)) // cos(11*Pi/24) -> 1/4*sqrt(8-2*sqrt(6)-2*sqrt(2)) 398 return sign*(_ex1_4*sqrt(_ex8 - _ex2*sqrt(_ex6) - _ex2*sqrt(_ex2))); 399 } 400 401 402 // Reflection at Pi/2 403 const ex ExOverPi = x_red/Pi; 404 if (is_exactly_a<numeric>(ExOverPi)) { 405 const numeric& c = ex_to<numeric>(ExOverPi); 406 // cos(integer*pi) --> (-1)^integer 407 if (c.is_integer()) 408 return pow(*_num_1_p, c); 409 if (c.is_rational()) { 410 numeric den = c.denom(); 411 numeric num = c.numer().mod(den * *_num2_p); 412 if (num > den) 413 num = den * *_num2_p - num; 414 if (num*(*_num2_p) > den) 415 return mul(_ex_1, cos_eval((den-num)*Pi/den)); 416 417 return cos((num*Pi)/den).hold(); 418 } 419 } 420 } 421 else 422 x_red = x; 423 424 // simplify cos(I*x) --> cosh(x) 425 if (is_multiple_of_I(x_red.expand())) 426 return cosh(x_red/I); 427 428 if (is_exactly_a<function>(x_red)) { 429 const ex &t = x_red.op(0); 430 431 // cos(acos(x)) -> x 432 if (is_ex_the_function(x_red, acos)) 433 return t; 434 435 // cos(asin(x)) -> sqrt(1-x^2) 436 if (is_ex_the_function(x_red, asin)) 437 return sqrt(_ex1-power(t,_ex2)); 438 439 // cos(atan(x)) -> 1/sqrt(1+x^2) 440 if (is_ex_the_function(x_red, atan)) 441 return power(_ex1+power(t,_ex2),_ex_1_2); 442 443 // cos(acot(x)) -> x/sqrt(1+x^2) 444 if (is_ex_the_function(x_red, acot)) 445 return t*power(_ex1+power(t,_ex2),_ex_1_2); 446 447 // cos(acsc(x)) -> sqrt(x^2-1)/x 448 if (is_ex_the_function(x_red, acsc)) 449 return sqrt(power(t,_ex2)-_ex1)/t; 450 451 // cos(asec(x)) -> 1/x 452 if (is_ex_the_function(x_red, asec)) 453 return _ex1/t; 454 455 // cos(atan2(y,x)) -> x/sqrt(x^2+y^2) 456 if (is_ex_the_function(x_red, atan2)) { 457 const ex &t1 = x_red.op(1); 458 return t1*power(power(t,_ex2)+power(t1,_ex2),_ex_1_2); 459 } 460 461 } 462 463 // cos(float) -> float 464 if (is_exactly_a<numeric>(x_red) && x_red.info(info_flags::inexact)) 465 return cos(ex_to<numeric>(x_red)); 466 467 return cos(x_red).hold(); 468 } 469 470 static ex cos_deriv(const ex & x, unsigned deriv_param) 471 { 472 GINAC_ASSERT(deriv_param==0); 473 474 // d/dx cos(x) -> -sin(x) 475 return -sin(x); 476 } 477 478 static ex cos_real_part(const ex & x) 479 { 480 return cosh(GiNaC::imag_part(x))*cos(GiNaC::real_part(x)); 481 } 482 483 static ex cos_imag_part(const ex & x) 484 { 485 return -sinh(GiNaC::imag_part(x))*sin(GiNaC::real_part(x)); 486 } 487 488 static ex cos_conjugate(const ex & x) 489 { 490 // conjugate(cos(x))==cos(conjugate(x)) 491 return cos(x.conjugate()); 492 } 493 494 REGISTER_FUNCTION(cos, eval_func(cos_eval). 495 derivative_func(cos_deriv). 496 real_part_func(cos_real_part). 497 imag_part_func(cos_imag_part). 498 conjugate_func(cos_conjugate). 499 latex_name("\\cos")); 500 501 ////////// 502 // tangent (trigonometric function) 503 ////////// 504 505 static ex tan_eval(const ex & x) 506 { 507 // tan(oo) -> error 508 // This should be before the tests below, since multiplying infinity 509 // with other values raises runtime_errors 510 if (x.info(info_flags::infinity)) { 511 throw (std::runtime_error("tan_eval(): tan(infinity) encountered")); 512 } 513 514 if (is_exactly_a<numeric>(x) and ex_to<numeric>(x).is_zero()) 515 return _ex0; 516 517 // tan() is odd 518 if (x.info(info_flags::negative)) 519 return -tan(-x); 520 521 ex x_red; 522 if (has_pi(x)) { 523 524 ex coef_pi = x.coeff(Pi,_ex1).expand(); 525 ex rem = _ex0; 526 if (is_exactly_a<add>(coef_pi)) { 527 for (size_t i=0; i < coef_pi.nops(); i++) { 528 if (coef_pi.op(i).is_integer()) 529 rem += Pi * coef_pi.op(i); 530 } 531 } 532 else if (coef_pi.is_integer()) 533 rem = Pi * coef_pi; 534 x_red = (x - rem).expand(); 535 536 537 // tan(n/d*Pi) -> { all known non-nested radicals } 538 const ex SixtyExOverPi = _ex60*x_red/Pi; 539 ex sign = _ex1; 540 if (is_exactly_a<numeric>(SixtyExOverPi) 541 and SixtyExOverPi.is_integer()) { 542 numeric z = mod(ex_to<numeric>(SixtyExOverPi),*_num60_p); 543 if (z>=*_num60_p) { 544 // wrap to interval [0, Pi) 545 z -= *_num60_p; 546 } 547 if (z>=*_num30_p) { 548 // wrap to interval [0, Pi/2) 549 z = *_num60_p-z; 550 sign = _ex_1; 551 } 552 if (z.is_equal(*_num0_p)) // tan(0) -> 0 553 return _ex0; 554 if (z.is_equal(*_num3_p)) // tan(Pi/20) -> sqrt(5) - sqrt(2*sqrt(5) + 5) + 1 555 return sign*(sqrt(_ex5)+_ex1-sqrt(_ex5+_ex2*sqrt(_ex5))); 556 if (z.is_equal(*_num5_p)) // tan(Pi/12) -> 2-sqrt(3) 557 return sign*(_ex2-sqrt(_ex3)); 558 if (z.is_equal(*_num6_p)) // tan(Pi/10) -> 1/5*sqrt(25-10*sqrt(5)) 559 return sign*sqrt(_ex25-_ex10*sqrt(_ex5))/_ex5; 560 if (z.is_equal(*_num9_p)) // tan(3*Pi/20) -> sqrt(5) - sqrt(-2*sqrt(5) + 5) - 1 561 return sign*(sqrt(_ex5)-_ex1-sqrt(_ex5-_ex2*sqrt(_ex5))); 562 if (z.is_equal(*_num10_p)) // tan(Pi/6) -> sqrt(3)/3 563 return sign*_ex1_3*sqrt(_ex3); 564 if (z.is_equal(*_num12_p)) // tan(Pi/5) -> sqrt(5-2*sqrt(5)) 565 return sign*sqrt(_ex5-_ex2*sqrt(_ex5)); 566 if (z.is_equal(*_num15_p)) // tan(Pi/4) -> 1 567 return sign; 568 if (z.is_equal(*_num18_p)) // tan(3*Pi/10) -> 1/5*sqrt(25+10*sqrt(5)) 569 return sign*sqrt(_ex25+_ex10*sqrt(_ex5))/_ex5; 570 if (z.is_equal(*_num20_p)) // tan(Pi/3) -> sqrt(3) 571 return sign*sqrt(_ex3); 572 if (z.is_equal(*_num21_p)) // tan(7*Pi/20) -> sqrt(5) + sqrt(-2*sqrt(5) + 5) - 1 573 return sign*(sqrt(_ex5)-_ex1+sqrt(_ex5-_ex2*sqrt(_ex5))); 574 if (z.is_equal(*_num24_p)) // tan(2*Pi/5) -> sqrt(5+2*sqrt(5)) 575 return sign*sqrt(_ex5+_ex2*sqrt(_ex5)); 576 if (z.is_equal(*_num25_p)) // tan(5/12*Pi) -> 2+sqrt(3) 577 return sign*(sqrt(_ex3)+_ex2); 578 if (z.is_equal(*_num27_p)) // tan(9*Pi/20) -> sqrt(5) + sqrt(2*sqrt(5) + 5) + 1 579 return sign*(sqrt(_ex5)+_ex1+sqrt(_ex5+_ex2*sqrt(_ex5))); 580 if (z.is_equal(*_num30_p)) // tan(Pi/2) -> infinity 581 //throw (pole_error("tan_eval(): simple pole",1)); 582 return UnsignedInfinity; 583 } 584 585 const ex FortyeightExOverPi = _ex48*x_red/Pi; 586 sign = _ex1; 587 if (is_exactly_a<numeric>(FortyeightExOverPi) 588 and FortyeightExOverPi.is_integer()) { 589 numeric z = mod(ex_to<numeric>(FortyeightExOverPi),*_num48_p); 590 if (z>=*_num48_p) { 591 // wrap to interval [0, Pi) 592 z -= *_num48_p; 593 } 594 if (z>*_num24_p) { 595 // wrap to interval [0, Pi/2) 596 z = *_num48_p-z; 597 sign = _ex_1; 598 } 599 if (z.is_equal(*_num2_p)) // tan(Pi/24) -> sqrt(6) - sqrt(-2*sqrt(6) + 5) - 2 600 return sign *(_ex_2+sqrt(_ex6)+sqrt(_ex2)-sqrt(_ex3)); 601 if (z.is_equal(*_num3_p)) // tan(Pi/16) -> -sqrt(2) + sqrt(2*sqrt(2) + 4) - 1 602 return sign*(_ex_1-sqrt(_ex2)+sqrt(_ex2*sqrt(_ex2)+_ex4)); 603 if (z.is_equal(*_num6_p)) // tan(Pi/8) -> sqrt(2)-1 604 return sign*(sqrt(_ex2)-_ex1); 605 if (z.is_equal(*_num9_p)) // tan(3*Pi/16) -> -sqrt(2) + sqrt(-2*sqrt(2) + 4) + 1 606 return sign*(_ex1-sqrt(_ex2)+sqrt(_ex4-_ex2*sqrt(_ex2))); 607 if (z.is_equal(*_num10_p)) // tan(5*Pi/24) -> sqrt(6) + sqrt(-2*sqrt(6) + 5) - 2 608 return sign*(_ex_2+sqrt(_ex6)-sqrt(_ex2)+sqrt(_ex3)); 609 if (z.is_equal(*_num14_p)) // tan(7*Pi/24) -> sqrt(6) - sqrt(2*sqrt(6) + 5) + 2 610 return sign*(_ex2+sqrt(_ex6)-sqrt(_ex2)-sqrt(_ex3)); 611 if (z.is_equal(*_num15_p)) // tan(5*Pi/16) -> sqrt(2) + sqrt(-2*sqrt(2) + 4) + 1 612 return sign*(_ex_1+sqrt(_ex2)+sqrt(_ex4-_ex2*sqrt(_ex2))); 613 if (z.is_equal(*_num18_p)) // tan(3*Pi/8) -> sqrt(2)+1 614 return sign*(sqrt(_ex2)+_ex1); 615 if (z.is_equal(*_num21_p)) // tan(7*Pi/16) -> sqrt(2) + sqrt(2*sqrt(2) + 4) + 1 616 return sign*(_ex1+sqrt(_ex2)+sqrt(_ex2*sqrt(_ex2)+_ex4)); 617 if (z.is_equal(*_num22_p)) // tan(11*Pi/24) -> sqrt(6) + sqrt(2*sqrt(6) + 5) + 2 618 return sign*(_ex2+sqrt(_ex6)+sqrt(_ex2)+sqrt(_ex3)); 619 } 620 621 // Reflection at Pi/2 622 const ex ExOverPi = x_red/Pi; 623 if (is_exactly_a<numeric>(ExOverPi)) { 624 const numeric& c = ex_to<numeric>(ExOverPi); 625 if (c.is_rational()) { 626 numeric den = c.denom(); 627 numeric num = c.numer().mod(den); 628 if (num*(*_num2_p) > den) 629 return mul(_ex_1, tan((den-num)*Pi/den).hold()); 630 return tan((num*Pi)/den).hold(); 631 } 632 } 633 } 634 else 635 x_red = x; 636 637 if (is_multiple_of_I(x_red.expand())) 638 return I*tanh(x_red/I); 639 640 if (is_exactly_a<function>(x_red)) { 641 const ex &t = x_red.op(0); 642 643 // tan(atan(x)) -> x 644 if (is_ex_the_function(x_red, atan)) 645 return t; 646 647 // tan(asin(x)) -> x/sqrt(1-x^2) 648 if (is_ex_the_function(x_red, asin)) 649 return t*power(_ex1-power(t,_ex2),_ex_1_2); 650 651 // tan(acos(x)) -> sqrt(1-x^2)/x 652 if (is_ex_the_function(x_red, acos)) 653 return power(t,_ex_1)*sqrt(_ex1-power(t,_ex2)); 654 655 // tan(acot(x)) -> 1/x 656 if (is_ex_the_function(x_red, acot)) 657 return _ex1/t; 658 659 // tan(acsc(x)) -> 1/sqrt(x^2-1) 660 if (is_ex_the_function(x_red, acsc)) 661 return power(power(t,_ex2)-_ex1,_ex_1_2); 662 663 // tan(asec(x)) -> sqrt(x^2-1) 664 if (is_ex_the_function(x_red, asec)) 665 return sqrt(power(t,_ex2)-_ex1); 666 667 // tan(atan2(y,x)) -> y/x 668 if (is_ex_the_function(x_red, atan2)) { 669 const ex &t1 = x_red.op(1); 670 return t/t1; 671 } 672 673 } 674 675 // tan(float) -> float 676 if (is_exactly_a<numeric>(x_red) && x_red.info(info_flags::inexact)) { 677 return tan(ex_to<numeric>(x_red)); 678 } 679 680 return tan(x_red).hold(); 681 } 682 683 static ex tan_deriv(const ex & x, unsigned deriv_param) 684 { 685 GINAC_ASSERT(deriv_param==0); 686 687 // d/dx tan(x) -> 1+tan(x)^2; 688 return (_ex1+power(tan(x),_ex2)); 689 } 690 691 static ex tan_real_part(const ex & x) 692 { 693 ex a = GiNaC::real_part(mul(x, _ex2)); 694 ex b = GiNaC::imag_part(mul(x, _ex2)); 695 return sin(a)/(cos(a) + cosh(b)); 696 } 697 698 static ex tan_imag_part(const ex & x) 699 { 700 ex a = GiNaC::real_part(mul(x, _ex2)); 701 ex b = GiNaC::imag_part(mul(x, _ex2)); 702 return sinh(b)/(cos(a) + cosh(b)); 703 } 704 705 static ex tan_series(const ex &x, 706 const relational &rel, 707 int order, 708 unsigned options) 709 { 710 GINAC_ASSERT(is_a<symbol>(rel.lhs())); 711 // method: 712 // Taylor series where there is no pole falls back to tan_deriv. 713 // On a pole simply expand sin(x)/cos(x). 714 const ex x_pt = x.subs(rel, subs_options::no_pattern); 715 if (!(2*x_pt/Pi).info(info_flags::odd)) 716 throw do_taylor(); // caught by function::series() 717 // if we got here we have to care for a simple pole 718 return (sin(x)/cos(x)).series(rel, order, options); 719 } 720 721 static ex tan_conjugate(const ex & x) 722 { 723 // conjugate(tan(x))==tan(conjugate(x)) 724 return tan(x.conjugate()); 725 } 726 727 REGISTER_FUNCTION(tan, eval_func(tan_eval). 728 derivative_func(tan_deriv). 729 series_func(tan_series). 730 real_part_func(tan_real_part). 731 imag_part_func(tan_imag_part). 732 conjugate_func(tan_conjugate). 733 latex_name("\\tan")); 734 735 ////////// 736 // cotangent (trigonometric function) 737 ////////// 738 739 static ex cot_eval(const ex & x) 740 { 741 // This should be before the tests below, since multiplying infinity 742 // with other values raises runtime_errors 743 if (x.is_zero()) 744 return UnsignedInfinity; 745 746 if (is_multiple_of_I(x.expand())) 747 return -I*coth(x/I); 748 749 if (is_exactly_a<function>(x)) { 750 const ex &t = x.op(0); 751 752 // cot(acot(x)) -> x 753 if (is_ex_the_function(x, acot)) 754 return t; 755 756 // cot(asin(x)) -> sqrt(1-x^2)/x 757 if (is_ex_the_function(x, asin)) 758 return sqrt(_ex1-power(t,_ex2))/t; 759 760 // cot(acos(x)) -> x/sqrt(1-x^2) 761 if (is_ex_the_function(x, acos)) 762 return t*power(_ex1-power(t,_ex2),_ex_1_2); 763 764 // cot(atan(x)) -> 1/x; 765 if (is_ex_the_function(x, atan)) 766 return _ex1/t; 767 768 // cot(acsc(x)) -> sqrt(x^2-1) 769 if (is_ex_the_function(x, acsc)) 770 return sqrt(power(t,_ex2)-_ex1); 771 772 // cot(asec(x)) -> 1/sqrt(x^2-1) 773 if (is_ex_the_function(x, asec)) 774 return power(power(t,_ex2)-_ex1,_ex_1_2); 775 776 // cot(atan2(y,x)) -> x/y; 777 if (is_ex_the_function(x, atan2)) { 778 const ex &t1 = x.op(1); 779 return t1/t; 780 } 781 } 782 783 // cot(float) -> float 784 if (is_exactly_a<numeric>(x) && x.info(info_flags::inexact)) { 785 if (ex_to<numeric>(x).is_zero()) 786 return UnsignedInfinity; 787 return tan(ex_to<numeric>(x)).inverse(); 788 } 789 790 ex res = tan_eval(x); 791 if (not is_ex_the_function(res, tan) && not is_ex_the_function(_ex_1*res, tan)) { 792 if (not res.is_zero()) { 793 return tan_eval(Pi/2-x); 794 } 795 796 return UnsignedInfinity; 797 } 798 // Reflection at Pi/2 799 const ex ExOverPi = x/Pi; 800 if(is_exactly_a<numeric>(ExOverPi)) { 801 ex coef_pi = x.coeff(Pi,_ex1).expand(); 802 if (is_exactly_a<numeric>(coef_pi)) { 803 const numeric& c = ex_to<numeric>(coef_pi); 804 if (c.is_rational()) { 805 const numeric num = c.numer(); 806 const numeric den = c.denom(); 807 const numeric rm = num.mod(den); 808 809 if (rm.mul(2) == den) { 810 return _ex0; 811 } 812 if (rm.mul(2) > den) { 813 return _ex_1*cot(den.sub(rm)*Pi/den).hold(); 814 } 815 return cot(rm*Pi/den).hold(); 816 } 817 } 818 } 819 return cot(x).hold(); 820 } 821 822 static ex cot_deriv(const ex & x, unsigned deriv_param) 823 { 824 GINAC_ASSERT(deriv_param==0); 825 826 // d/dx cot(x) -> -1-cot(x)^2; 827 return (_ex_1-power(cot(x),_ex2)); 828 } 829 830 // Ref.: http://dlmf.nist.gov/4.21.E40 831 static ex cot_real_part(const ex & x) 832 { 833 ex a = GiNaC::real_part(mul(x, _ex2)); 834 ex b = GiNaC::imag_part(mul(x, _ex2)); 835 return sin(a)/(cosh(b) - cos(a)); 836 } 837 838 static ex cot_imag_part(const ex & x) 839 { 840 ex a = GiNaC::real_part(mul(x, _ex2)); 841 ex b = GiNaC::imag_part(mul(x, _ex2)); 842 return -sinh(b)/(cosh(b) - cos(a)); 843 } 844 845 static ex cot_series(const ex &x, 846 const relational &rel, 847 int order, 848 unsigned options) 849 { 850 GINAC_ASSERT(is_a<symbol>(rel.lhs())); 851 // method: 852 // Taylor series where there is no pole falls back to tan_deriv. 853 // On a pole simply expand cos(x)/sin(x). 854 const ex x_pt = x.subs(rel, subs_options::no_pattern); 855 if (!(2*x_pt/Pi).info(info_flags::even)) 856 throw do_taylor(); // caught by function::series() 857 // if we got here we have to care for a simple pole 858 return (cos(x)/sin(x)).series(rel, order, options); 859 } 860 861 static ex cot_conjugate(const ex & x) 862 { 863 // conjugate(tan(x))==1/tan(conjugate(x)) 864 return power(tan(x.conjugate()), _ex_1); 865 } 866 867 REGISTER_FUNCTION(cot, eval_func(cot_eval). 868 derivative_func(cot_deriv). 869 series_func(cot_series). 870 real_part_func(cot_real_part). 871 imag_part_func(cot_imag_part). 872 conjugate_func(cot_conjugate). 873 latex_name("\\cot")); 874 875 ////////// 876 // secant (trigonometric function) 877 ////////// 878 879 static ex sec_eval(const ex & x) 880 { 881 if (is_multiple_of_I(x.expand())) 882 return sech(x/I); 883 884 if (is_exactly_a<function>(x)) { 885 const ex &t = x.op(0); 886 887 // sec(asec(x)) -> x 888 if (is_ex_the_function(x, asec)) 889 return t; 890 891 // sec(asin(x)) -> 1/sqrt(1-x^2) 892 if (is_ex_the_function(x, asin)) 893 return power(_ex1-power(t,_ex2),_ex_1_2); 894 895 // sec(acos(x)) -> 1/x 896 if (is_ex_the_function(x, acos)) 897 return _ex1/t; 898 899 // sec(atan(x)) -> sqrt(x^2+1) 900 if (is_ex_the_function(x, atan)) 901 return sqrt(power(t,_ex2)+_ex1); 902 903 // sec(acot(x)) -> sqrt(x^2+1)/x 904 if (is_ex_the_function(x, acot)) 905 return sqrt(power(t,_ex2)+_ex1)/t; 906 907 // sec(acsch(x)) -> x/sqrt(x^2-1) 908 if (is_ex_the_function(x, acsc)) 909 return t*power(power(t,_ex2)-_ex1,_ex_1_2); 910 911 // sec(atan2(y,x)) -> sqrt(x^2+y^2)/x 912 if (is_ex_the_function(x, atan2)) { 913 const ex &t1 = x.op(1); 914 return sqrt(power(t1,_ex1)+power(t,_ex2))/t1; 915 } 916 917 } 918 919 // sec() is even 920 if (x.info(info_flags::negative)) 921 return sec(-x); 922 923 // sec(float) -> float 924 if (is_exactly_a<numeric>(x) && x.info(info_flags::inexact)) { 925 return cos(ex_to<numeric>(x)).inverse(); 926 } 927 928 // Handle simplification via cos 929 ex res = cos_eval(x); 930 if (not is_ex_the_function(res, cos) && not is_ex_the_function(_ex_1*res, cos)) { 931 if (res.is_zero()) 932 return UnsignedInfinity; 933 return power(res, _ex_1); 934 } 935 // cos has reflected also the argument so take it 936 if (is_ex_the_function(res, cos)) 937 return sec(res.op(0)).hold(); 938 939 return -sec((-res).op(0)).hold(); 940 } 941 942 static ex sec_deriv(const ex & x, unsigned deriv_param) 943 { 944 GINAC_ASSERT(deriv_param==0); 945 946 // d/dx sec(x) -> sec(x)*tan(x); 947 return sec(x)*tan(x); 948 } 949 950 static ex sec_real_part(const ex & x) 951 { 952 ex a = GiNaC::real_part(x); 953 ex b = GiNaC::imag_part(x); 954 return cos(a)*cosh(b)/(sin(a)*sin(a)*sinh(b)*sinh(b) \ 955 + cos(a)*cos(a)*cosh(b)*cosh(b)); 956 } 957 958 static ex sec_imag_part(const ex & x) 959 { 960 ex a = GiNaC::real_part(x); 961 ex b = GiNaC::imag_part(x); 962 return sin(a)*sinh(b)/(sin(a)*sin(a)*sinh(b)*sinh(b) \ 963 + cos(a)*cos(a)*cosh(b)*cosh(b)); 964 } 965 966 static ex sec_series(const ex &x, 967 const relational &rel, 968 int order, 969 unsigned options) 970 { 971 GINAC_ASSERT(is_a<symbol>(rel.lhs())); 972 return (_ex1/cos(x)).series(rel, order, options); 973 } 974 975 static ex sec_conjugate(const ex & x) 976 { 977 // conjugate(tan(x))==1/tan(conjugate(x)) 978 return power(cos(x.conjugate()), _ex_1); 979 } 980 981 REGISTER_FUNCTION(sec, eval_func(sec_eval). 982 derivative_func(sec_deriv). 983 series_func(sec_series). 984 real_part_func(sec_real_part). 985 imag_part_func(sec_imag_part). 986 conjugate_func(sec_conjugate). 987 latex_name("\\sec")); 988 989 ////////// 990 // cosecant (trigonometric function) 991 ////////// 992 993 static ex csc_eval(const ex & x) 994 { 995 996 if (is_multiple_of_I(x.expand())) 997 return -I*csch(x/I); 998 999 if (is_exactly_a<function>(x)) { 1000 const ex &t = x.op(0); 1001 1002 // csc(acsc(x)) -> x 1003 if (is_ex_the_function(x, acsc)) 1004 return t; 1005 1006 // csc(asin(x)) -> 1/x 1007 if (is_ex_the_function(x, asin)) 1008 return _ex1/t; 1009 1010 // csc(acos(x)) -> 1/sqrt(1-x^2) 1011 if (is_ex_the_function(x, acos)) 1012 return power(_ex1-power(t,_ex2),_ex_1_2); 1013 1014 // csc(atan(x)) -> sqrt(x^2+1)/x; 1015 if (is_ex_the_function(x, atan)) 1016 return sqrt(power(t,_ex2)+_ex1)/t; 1017 1018 // csc(acot(x)) -> sqrt(x^2+1); 1019 if (is_ex_the_function(x, acot)) 1020 return sqrt(power(t,_ex2)+_ex1); 1021 1022 // csc(asec(x)) -> x/sqrt(x^2-1) 1023 if (is_ex_the_function(x, asec)) 1024 return t*power(power(t,_ex2)-_ex1,_ex_1_2); 1025 1026 // csc(atan2(y,x)) -> sqrt(x^2+y^2)/y 1027 if (is_ex_the_function(x, atan2)) { 1028 const ex &t1 = x.op(1); 1029 return sqrt(power(t1,_ex1)+power(t,_ex2))/t; 1030 } 1031 1032 } 1033 1034 // csc(float) -> float 1035 if (is_exactly_a<numeric>(x) && x.info(info_flags::inexact)) { 1036 if (ex_to<numeric>(x).is_zero()) 1037 return UnsignedInfinity; 1038 return sin(ex_to<numeric>(x)).inverse(); 1039 } 1040 1041 // Handle simplification via sin 1042 ex res = sin_eval(x); 1043 if (not is_ex_the_function(res, sin) && not is_ex_the_function(_ex_1*res, sin)) { 1044 if (res.is_zero()) 1045 return UnsignedInfinity; 1046 1047 return power(res, _ex_1); 1048 } 1049 // sin has reflected also the argument so take it 1050 if (is_ex_the_function(res, sin)) 1051 return csc(res.op(0)).hold(); 1052 return -csc((-res).op(0)).hold(); 1053 } 1054 1055 static ex csc_deriv(const ex & x, unsigned deriv_param) 1056 { 1057 GINAC_ASSERT(deriv_param==0); 1058 1059 // d/dx cot(x) -> -1-cot(x)^2; 1060 return -csc(x)*cot(x); 1061 } 1062 1063 static ex csc_real_part(const ex & x) 1064 { 1065 ex a = GiNaC::real_part(x); 1066 ex b = GiNaC::imag_part(x); 1067 return sin(a)*cosh(b)/(sin(a)*sin(a)*cosh(b)*cosh(b) \ 1068 + cos(a)*cos(a)*sinh(b)*sinh(b)); 1069 } 1070 1071 static ex csc_imag_part(const ex & x) 1072 { 1073 ex a = GiNaC::real_part(x); 1074 ex b = GiNaC::imag_part(x); 1075 return -cos(a)*sinh(b)/(sin(a)*sin(a)*cosh(b)*cosh(b) \ 1076 + cos(a)*cos(a)*sinh(b)*sinh(b)); 1077 } 1078 1079 static ex csc_series(const ex &x, 1080 const relational &rel, 1081 int order, 1082 unsigned options) 1083 { 1084 GINAC_ASSERT(is_a<symbol>(rel.lhs())); 1085 return (_ex1/sin(x)).series(rel, order, options); 1086 } 1087 1088 static ex csc_conjugate(const ex & x) 1089 { 1090 // conjugate(tan(x))==1/tan(conjugate(x)) 1091 return power(sin(x.conjugate()), _ex_1); 1092 } 1093 1094 REGISTER_FUNCTION(csc, eval_func(csc_eval). 1095 derivative_func(csc_deriv). 1096 series_func(csc_series). 1097 real_part_func(csc_real_part). 1098 imag_part_func(csc_imag_part). 1099 conjugate_func(csc_conjugate). 1100 latex_name("\\csc")); 1101 1102 ////////// 1103 // inverse sine (arc sine) 1104 ////////// 1105 1106 // Needed because there is no RR member 1107 static ex asin_evalf(const ex & x, PyObject* parent) 1108 { 1109 if (is_exactly_a<numeric>(x)) 1110 return asin(ex_to<numeric>(x), parent); 1111 1112 return asin(x).hold(); 1113 } 1114 1115 static ex asin_eval(const ex & x) 1116 { 1117 // asin() is odd 1118 if (x.info(info_flags::negative)) 1119 return -asin(-x); 1120 1121 if (is_exactly_a<numeric>(x)) { 1122 // asin(0) -> 0 1123 if (x.is_zero()) 1124 return x; 1125 1126 // asin(1/2) -> Pi/6 1127 if (x.is_equal(_ex1_2)) 1128 return numeric(1,6)*Pi; 1129 1130 // asin(1) -> Pi/2 1131 if (x.is_one()) 1132 return _ex1_2*Pi; 1133 1134 if (x.info(info_flags::inexact)) { 1135 const numeric& num = ex_to<numeric>(x); 1136 1137 // asin(float) -> float 1138 return asin(num); 1139 } 1140 1141 } 1142 1143 // asin(oo) -> error 1144 // asin(UnsignedInfinity) -> UnsignedInfinity 1145 if (x.info(info_flags::infinity)) { 1146 if (x.is_equal(UnsignedInfinity)) 1147 return UnsignedInfinity; 1148 throw (std::runtime_error("arcsin_eval(): arcsin(infinity) encountered")); 1149 } 1150 1151 if (x.is_equal(mul(pow(_ex2, _ex1_2), _ex1_2))) 1152 return mul(Pi, _ex1_4); 1153 1154 if (x.is_equal(mul(pow(_ex3, _ex1_2), _ex1_2))) 1155 return mul(Pi, _ex1_3); 1156 1157 return asin(x).hold(); 1158 } 1159 1160 static ex asin_deriv(const ex & x, unsigned deriv_param) 1161 { 1162 GINAC_ASSERT(deriv_param==0); 1163 1164 // d/dx asin(x) -> 1/sqrt(1-x^2) 1165 return power(1-power(x,_ex2),_ex_1_2); 1166 } 1167 1168 static ex asin_conjugate(const ex & x) 1169 { 1170 // conjugate(asin(x))==asin(conjugate(x)) unless on the branch cuts which 1171 // run along the real axis outside the interval [-1, +1]. 1172 if (is_exactly_a<numeric>(x) && 1173 (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) { 1174 return asin(x.conjugate()); 1175 } 1176 return conjugate_function(asin(x)).hold(); 1177 } 1178 1179 REGISTER_FUNCTION(asin, eval_func(asin_eval). 1180 evalf_func(asin_evalf). 1181 derivative_func(asin_deriv). 1182 conjugate_func(asin_conjugate). 1183 set_name("arcsin", "\\arcsin")); 1184 1185 ////////// 1186 // inverse cosine (arc cosine) 1187 ////////// 1188 1189 static ex acos_eval(const ex & x) 1190 { 1191 if (is_exactly_a<numeric>(x)) { 1192 1193 // acos(1) -> 0 1194 if (x.is_one()) 1195 return _ex0; 1196 1197 // acos(1/2) -> Pi/3 1198 if (x.is_equal(_ex1_2)) 1199 return _ex1_3*Pi; 1200 1201 // acos(0) -> Pi/2 1202 if (x.is_zero()) 1203 return _ex1_2*Pi; 1204 1205 // acos(-1/2) -> 2/3*Pi 1206 if (x.is_equal(_ex_1_2)) 1207 return numeric(2,3)*Pi; 1208 1209 // acos(-1) -> Pi 1210 if (x.is_minus_one()) 1211 return Pi; 1212 1213 if (x.info(info_flags::inexact)) { 1214 const numeric& num = ex_to<numeric>(x); 1215 1216 // acos(float) -> float 1217 return acos(num); 1218 } 1219 1220 // acos(-x) -> Pi-acos(x) 1221 if (x.info(info_flags::negative)) 1222 return Pi-acos(-x); 1223 } 1224 1225 // acos(oo) -> error 1226 // acos(UnsignedInfinity) -> UnsignedInfinity 1227 if (x.info(info_flags::infinity)) { 1228 if (x.is_equal(UnsignedInfinity)) 1229 return UnsignedInfinity; 1230 throw (std::runtime_error("arccos_eval(): arccos(infinity) encountered")); 1231 } 1232 return acos(x).hold(); 1233 } 1234 1235 static ex acos_deriv(const ex & x, unsigned deriv_param) 1236 { 1237 GINAC_ASSERT(deriv_param==0); 1238 1239 // d/dx acos(x) -> -1/sqrt(1-x^2) 1240 return -power(1-power(x,_ex2),_ex_1_2); 1241 } 1242 1243 static ex acos_conjugate(const ex & x) 1244 { 1245 // conjugate(acos(x))==acos(conjugate(x)) unless on the branch cuts which 1246 // run along the real axis outside the interval [-1, +1]. 1247 if (is_exactly_a<numeric>(x) && 1248 (!x.imag_part().is_zero() || (x > *_num_1_p && x < *_num1_p))) { 1249 return acos(x.conjugate()); 1250 } 1251 return conjugate_function(acos(x)).hold(); 1252 } 1253 1254 1255 REGISTER_FUNCTION(acos, eval_func(acos_eval). 1256 derivative_func(acos_deriv). 1257 conjugate_func(acos_conjugate). 1258 set_name("arccos", "\\arccos")); 1259 1260 ////////// 1261 // inverse tangent (arc tangent) 1262 ////////// 1263 1264 static ex atan_eval(const ex & x) 1265 { 1266 // atan() is odd 1267 if (x.info(info_flags::negative)) 1268 return -atan(-x); 1269 1270 if (is_exactly_a<numeric>(x)) { 1271 // atan(0) -> 0 1272 if (x.is_zero()) 1273 return _ex0; 1274 1275 // atan(1) -> Pi/4 1276 if (x.is_one()) 1277 return _ex1_4*Pi; 1278 1279 if (x.is_equal(I)) 1280 throw (pole_error("atan_eval(): logarithmic pole",0)); 1281 1282 // atan(float) -> float 1283 if (x.info(info_flags::inexact)) 1284 return atan(ex_to<numeric>(x)); 1285 } 1286 1287 // arctan(oo) -> Pi/2 1288 // arctan(-oo) -> -Pi/2 1289 // arctan(UnsignedInfinity) -> error 1290 if (is_exactly_a<infinity>(x)) { 1291 const infinity& xinf = ex_to<infinity>(x); 1292 if (xinf.is_plus_infinity()) 1293 return _ex1_2*Pi; 1294 if (xinf.is_minus_infinity()) 1295 return _ex_1_2*Pi; 1296 // x is UnsignedInfinity 1297 throw (std::runtime_error("arctan_eval(): arctan(unsigned_infinity) encountered")); 1298 } 1299 1300 if (x.is_equal(pow(_ex3, _ex1_2))) 1301 return mul(Pi, _ex1_3); 1302 1303 if (x.is_equal(mul(pow(_ex3, _ex1_2), _ex1_3))) 1304 return mul(Pi, numeric(1,6)); 1305 1306 return atan(x).hold(); 1307 } 1308 1309 static ex atan_deriv(const ex & x, unsigned deriv_param) 1310 { 1311 GINAC_ASSERT(deriv_param==0); 1312 1313 // d/dx atan(x) -> 1/(1+x^2) 1314 return power(_ex1+power(x,_ex2), _ex_1); 1315 } 1316 1317 static ex atan_series(const ex &arg, 1318 const relational &rel, 1319 int order, 1320 unsigned options) 1321 { 1322 GINAC_ASSERT(is_a<symbol>(rel.lhs())); 1323 // method: 1324 // Taylor series where there is no pole or cut falls back to atan_deriv. 1325 // There are two branch cuts, one running from I up the imaginary axis and 1326 // one running from -I down the imaginary axis. The points I and -I are 1327 // poles. 1328 // On the branch cuts and the poles series expand 1329 // (log(1+I*x)-log(1-I*x))/(2*I) 1330 // instead. 1331 const ex arg_pt = arg.subs(rel, subs_options::no_pattern); 1332 if (not (I*arg_pt).is_real()) 1333 throw do_taylor(); // Re(x) != 0 1334 if ((I*arg_pt).is_real() && abs(I*arg_pt)<_ex1) 1335 throw do_taylor(); // Re(x) == 0, but abs(x)<1 1336 // care for the poles, using the defining formula for atan()... 1337 if (arg_pt.is_equal(I) || arg_pt.is_equal(-I)) 1338 return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options); 1339 if ((options & series_options::suppress_branchcut) == 0u) { 1340 // method: 1341 // This is the branch cut: assemble the primitive series manually and 1342 // then add the corresponding complex step function. 1343 const symbol &s = ex_to<symbol>(rel.lhs()); 1344 const ex &point = rel.rhs(); 1345 const symbol foo; 1346 const ex replarg = series(atan(arg), s==foo, order).subs(foo==point, subs_options::no_pattern); 1347 ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2; 1348 if ((I*arg_pt)<_ex0) 1349 Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2; 1350 else 1351 Order0correction += log((I*arg_pt+_ex1)/(I*arg_pt+_ex_1))*I*_ex1_2; 1352 epvector seq; 1353 seq.emplace_back(Order0correction, _ex0); 1354 seq.emplace_back(Order(_ex1), order); 1355 return series(replarg - pseries(rel, seq), rel, order); 1356 } 1357 throw do_taylor(); 1358 } 1359 1360 static ex atan_conjugate(const ex & x) 1361 { 1362 // conjugate(atan(x))==atan(conjugate(x)) unless on the branch cuts which 1363 // run along the imaginary axis outside the interval [-I, +I]. 1364 if (x.is_real()) 1365 return atan(x); 1366 if (is_exactly_a<numeric>(x)) { 1367 const numeric x_re = ex_to<numeric>(x.real_part()); 1368 const numeric x_im = ex_to<numeric>(x.imag_part()); 1369 if (!x_re.is_zero() || 1370 (x_im > *_num_1_p && x_im < *_num1_p)) 1371 return atan(x.conjugate()); 1372 } 1373 return conjugate_function(atan(x)).hold(); 1374 } 1375 1376 REGISTER_FUNCTION(atan, eval_func(atan_eval). 1377 derivative_func(atan_deriv). 1378 series_func(atan_series). 1379 conjugate_func(atan_conjugate). 1380 set_name("arctan", "\\arctan")); 1381 1382 ////////// 1383 // inverse tangent (atan2(y,x)) 1384 ////////// 1385 1386 static ex atan2_evalf(const ex &y, const ex &x, PyObject* parent) 1387 { 1388 if (is_exactly_a<numeric>(y) and is_exactly_a<numeric>(x)) 1389 return atan(ex_to<numeric>(y), ex_to<numeric>(x)); 1390 1391 return atan2(y, x).hold(); 1392 } 1393 1394 static ex atan2_eval(const ex & y, const ex & x) 1395 { 1396 if (y.is_zero()) { 1397 1398 // atan2(0, 0) -> undefined 1399 if (x.is_zero()) 1400 return NaN; 1401 1402 // atan2(0, x), x real and positive -> 0 1403 if (x.is_positive()) 1404 return _ex0; 1405 1406 // atan2(0, x), x real and negative -> Pi 1407 if (x.info(info_flags::negative)) 1408 return Pi; 1409 } 1410 1411 if (x.is_zero()) { 1412 1413 // atan2(y, 0), y real and positive -> Pi/2 1414 if (y.is_positive()) 1415 return _ex1_2*Pi; 1416 1417 // atan2(y, 0), y real and negative -> -Pi/2 1418 if (y.info(info_flags::negative)) 1419 return _ex_1_2*Pi; 1420 } 1421 1422 if (y.is_equal(x)) { 1423 1424 // atan2(y, y), y real and positive -> Pi/4 1425 if (y.is_positive()) 1426 return _ex1_4*Pi; 1427 1428 // atan2(y, y), y real and negative -> -3/4*Pi 1429 if (y.info(info_flags::negative)) 1430 return numeric(-3, 4)*Pi; 1431 } 1432 1433 if (y.is_equal(-x)) { 1434 1435 // atan2(y, -y), y real and positive -> 3*Pi/4 1436 if (y.is_positive()) 1437 return numeric(3, 4)*Pi; 1438 1439 // atan2(y, -y), y real and negative -> -Pi/4 1440 if (y.info(info_flags::negative)) 1441 return _ex_1_4*Pi; 1442 } 1443 1444 // atan2(float, float) -> float 1445 if (is_exactly_a<numeric>(x) 1446 and is_exactly_a<numeric>(y) 1447 and (x.info(info_flags::inexact) or y.info(info_flags::inexact))) 1448 return atan(ex_to<numeric>(y), ex_to<numeric>(x)); 1449 1450 // handle infinities 1451 if (is_a<infinity>(x) || is_a<infinity>(y)) { 1452 if (is_a<infinity>(x) && ex_to<infinity>(x).is_unsigned_infinity()) 1453 throw (std::runtime_error("arctan2_eval(): arctan2(unsigned_infinity, x) encountered")); 1454 if (is_a<infinity>(y) && ex_to<infinity>(y).is_unsigned_infinity()) 1455 throw (std::runtime_error("arctan2_eval(): arctan2(x, unsigned_infinity) encountered")); 1456 1457 if (is_a<infinity>(x) && is_a<infinity>(y)) 1458 return atan2_eval(ex_to<infinity>(x).get_direction(), 1459 ex_to<infinity>(y).get_direction()); 1460 1461 if (is_a<infinity>(x)) 1462 return atan2_eval(ex_to<infinity>(x).get_direction(), 0); 1463 if (is_a<infinity>(y)) 1464 return atan2_eval(0, ex_to<infinity>(y).get_direction()); 1465 } 1466 1467 // atan2(real, real) -> atan(y/x) +/- Pi 1468 if (y.is_real() && x.is_real()) { 1469 if (x.is_positive()) 1470 return atan(y/x); 1471 1472 if (x.info(info_flags::negative)) { 1473 if (y.is_positive()) 1474 return atan(y/x)+Pi; 1475 if (y.info(info_flags::negative)) 1476 return atan(y/x)-Pi; 1477 } 1478 } 1479 1480 return atan2(y, x).hold(); 1481 } 1482 1483 static ex atan2_deriv(const ex & y, const ex & x, unsigned deriv_param) 1484 { 1485 GINAC_ASSERT(deriv_param<2); 1486 1487 if (deriv_param==0) { 1488 // d/dy atan2(y,x) 1489 return x*power(power(x,_ex2)+power(y,_ex2),_ex_1); 1490 } 1491 // d/dx atan2(y,x) 1492 return -y*power(power(x,_ex2)+power(y,_ex2),_ex_1); 1493 } 1494 1495 REGISTER_FUNCTION(atan2, eval_func(atan2_eval). 1496 evalf_func(atan2_evalf). 1497 derivative_func(atan2_deriv). 1498 set_name("arctan2")); 1499 1500 ////////// 1501 // inverse cotangent (arc cotangent) 1502 ////////// 1503 1504 // Needed because there is no Python RR equivalent 1505 static ex acot_evalf(const ex & x, PyObject* parent) 1506 { 1507 if (is_exactly_a<numeric>(x)) 1508 return atan(ex_to<numeric>(x).inverse()); 1509 1510 return acot(x).hold(); 1511 } 1512 1513 static ex acot_eval(const ex & x) 1514 { 1515 if (is_exactly_a<numeric>(x)) { 1516 1517 if (x.is_zero()) 1518 return _ex1_2*Pi; 1519 1520 if (x.is_one()) 1521 return _ex1_4*Pi; 1522 1523 if (x.is_minus_one()) 1524 return _ex_1_4*Pi; 1525 1526 if (x.is_equal(I) || x.is_equal(-I)) 1527 throw (pole_error("acot_eval(): logarithmic pole",0)); 1528 1529 if (x.info(info_flags::inexact)) 1530 return atan(ex_to<numeric>(x).inverse()); 1531 1532 if (x.info(info_flags::negative)) 1533 return -acot(-x); 1534 } 1535 1536 if (x.info(info_flags::infinity)) { 1537 return _ex0; 1538 } 1539 1540 if (is_exactly_a<function>(x)) { 1541 const ex &t = x.op(0); 1542 if (is_ex_the_function(x, cot)) 1543 return t; 1544 } 1545 return acot(x).hold(); 1546 } 1547 1548 static ex acot_deriv(const ex & x, unsigned deriv_param) 1549 { 1550 GINAC_ASSERT(deriv_param==0); 1551 return mul(power(_ex1+power(x,_ex2), _ex_1), _ex_1); 1552 } 1553 1554 static ex acot_series(const ex &arg, 1555 const relational &rel, 1556 int order, 1557 unsigned options) 1558 { 1559 return _ex1_2*Pi - atan_series(arg, rel, order, options); 1560 } 1561 1562 REGISTER_FUNCTION(acot, eval_func(acot_eval). 1563 evalf_func(acot_evalf). 1564 derivative_func(acot_deriv). 1565 series_func(acot_series). 1566 set_name("arccot", "\\operatorname{arccot}")); 1567 1568 1569 ////////// 1570 // inverse secant (arc secant) 1571 ////////// 1572 1573 // Needed because there is no Python RR equivalent 1574 static ex asec_evalf(const ex & x, PyObject* parent) 1575 { 1576 if (is_exactly_a<numeric>(x)) 1577 { 1578 const numeric& num = ex_to<numeric>(x); 1579 return acos(num.inverse()); 1580 } 1581 1582 return asec(x).hold(); 1583 } 1584 1585 static ex asec_eval(const ex & x) 1586 { 1587 if (is_exactly_a<numeric>(x)) { 1588 const numeric& num = ex_to<numeric>(x); 1589 if (num.is_zero()) 1590 return NaN; 1591 if (num.is_equal(*_num1_p)) 1592 return _ex0; 1593 if (num.is_equal(*_num_1_p)) 1594 return Pi; 1595 if (num.info(info_flags::inexact)) 1596 return asec_evalf(x, nullptr); 1597 } 1598 1599 if (x.info(info_flags::infinity)) { 1600 return Pi/_ex2; 1601 } 1602 1603 if (is_exactly_a<function>(x)) { 1604 const ex &t = x.op(0); 1605 if (is_ex_the_function(x, sec)) 1606 return t; 1607 } 1608 return asec(x).hold(); 1609 } 1610 1611 static ex asec_deriv(const ex & x, unsigned deriv_param) 1612 { 1613 GINAC_ASSERT(deriv_param==0); 1614 return power(mul(x, power(_ex_1 + power(x,_ex2), _ex1_2)), _ex_1); 1615 } 1616 1617 static ex asec_series(const ex &arg, 1618 const relational &rel, 1619 int order, 1620 unsigned options) 1621 { 1622 const ex &point = rel.rhs(); 1623 if (not point.info(info_flags::infinity)) { 1624 throw (pole_error("cannot expand asec(x) around finite value",0)); 1625 } 1626 1627 throw (std::runtime_error("series growing in 1/x not implemented")); 1628 } 1629 1630 REGISTER_FUNCTION(asec, eval_func(asec_eval). 1631 evalf_func(asec_evalf). 1632 derivative_func(asec_deriv). 1633 series_func(asec_series). 1634 set_name("arcsec", "\\operatorname{arcsec}")); 1635 1636 1637 ////////// 1638 // inverse cosecant (arc cosecant) 1639 ////////// 1640 1641 // Needed because there is no Python RR equivalent 1642 static ex acsc_evalf(const ex & x, PyObject* parent) 1643 { 1644 if (is_exactly_a<numeric>(x)) { 1645 const numeric& num = ex_to<numeric>(x); 1646 return asin(num.inverse()); 1647 } 1648 1649 return acsc(x).hold(); 1650 } 1651 1652 static ex acsc_eval(const ex & x) 1653 { 1654 if (is_exactly_a<numeric>(x)) { 1655 const numeric& num = ex_to<numeric>(x); 1656 if (num.is_zero()) 1657 return NaN; 1658 if (num.is_equal(*_num1_p)) 1659 return Pi/_ex2; 1660 if (num.is_equal(*_num_1_p)) 1661 return -Pi/_ex2; 1662 if (num.info(info_flags::inexact)) 1663 return asin(num.inverse()); 1664 } 1665 1666 if (x.info(info_flags::infinity)) { 1667 return _ex0; 1668 } 1669 1670 if (is_exactly_a<function>(x)) { 1671 const ex &t = x.op(0); 1672 if (is_ex_the_function(x, csc)) 1673 return t; 1674 } 1675 return acsc(x).hold(); 1676 } 1677 1678 static ex acsc_deriv(const ex & x, unsigned deriv_param) 1679 { 1680 GINAC_ASSERT(deriv_param==0); 1681 return -power(mul(x, power(_ex_1 + power(x,_ex2), _ex1_2)), _ex_1); 1682 } 1683 1684 static ex acsc_series(const ex &arg, 1685 const relational &rel, 1686 int order, 1687 unsigned options) 1688 { 1689 const ex &point = rel.rhs(); 1690 if (not point.info(info_flags::infinity)) { 1691 throw (pole_error("cannot expand acsc(x) around finite value",0)); 1692 } 1693 1694 throw (std::runtime_error("series growing in 1/x not implemented")); 1695 } 1696 1697 REGISTER_FUNCTION(acsc, eval_func(acsc_eval). 1698 evalf_func(acsc_evalf). 1699 derivative_func(acsc_deriv). 1700 series_func(acsc_series). 1701 set_name("arccsc", "\\operatorname{arccsc}")); 1702 1703 1704 } // namespace GiNaC 1705