1"""Power series evaluation and manipulation using sparse Polynomials 2 3Implementing a new function 4--------------------------- 5 6There are a few things to be kept in mind when adding a new function here:: 7 8 - The implementation should work on all possible input domains/rings. 9 Special cases include the ``EX`` ring and a constant term in the series 10 to be expanded. There can be two types of constant terms in the series: 11 12 + A constant value or symbol. 13 + A term of a multivariate series not involving the generator, with 14 respect to which the series is to expanded. 15 16 Strictly speaking, a generator of a ring should not be considered a 17 constant. However, for series expansion both the cases need similar 18 treatment (as the user doesn't care about inner details), i.e, use an 19 addition formula to separate the constant part and the variable part (see 20 rs_sin for reference). 21 22 - All the algorithms used here are primarily designed to work for Taylor 23 series (number of iterations in the algo equals the required order). 24 Hence, it becomes tricky to get the series of the right order if a 25 Puiseux series is input. Use rs_puiseux? in your function if your 26 algorithm is not designed to handle fractional powers. 27 28Extending rs_series 29------------------- 30 31To make a function work with rs_series you need to do two things:: 32 33 - Many sure it works with a constant term (as explained above). 34 - If the series contains constant terms, you might need to extend its ring. 35 You do so by adding the new terms to the rings as generators. 36 ``PolyRing.compose`` and ``PolyRing.add_gens`` are two functions that do 37 so and need to be called every time you expand a series containing a 38 constant term. 39 40Look at rs_sin and rs_series for further reference. 41 42""" 43 44from sympy.polys.domains import QQ, EX 45from sympy.polys.rings import PolyElement, ring, sring 46from sympy.polys.polyerrors import DomainError 47from sympy.polys.monomials import (monomial_min, monomial_mul, monomial_div, 48 monomial_ldiv) 49from mpmath.libmp.libintmath import ifac 50from sympy.core import PoleError, Function, Expr 51from sympy.core.numbers import Rational, igcd 52from sympy.core.compatibility import as_int 53from sympy.functions import sin, cos, tan, atan, exp, atanh, tanh, log, ceiling 54from mpmath.libmp.libintmath import giant_steps 55import math 56 57 58def _invert_monoms(p1): 59 """ 60 Compute ``x**n * p1(1/x)`` for a univariate polynomial ``p1`` in ``x``. 61 62 Examples 63 ======== 64 65 >>> from sympy.polys.domains import ZZ 66 >>> from sympy.polys.rings import ring 67 >>> from sympy.polys.ring_series import _invert_monoms 68 >>> R, x = ring('x', ZZ) 69 >>> p = x**2 + 2*x + 3 70 >>> _invert_monoms(p) 71 3*x**2 + 2*x + 1 72 73 See Also 74 ======== 75 76 sympy.polys.densebasic.dup_reverse 77 """ 78 terms = list(p1.items()) 79 terms.sort() 80 deg = p1.degree() 81 R = p1.ring 82 p = R.zero 83 cv = p1.listcoeffs() 84 mv = p1.listmonoms() 85 for i in range(len(mv)): 86 p[(deg - mv[i][0],)] = cv[i] 87 return p 88 89def _giant_steps(target): 90 """Return a list of precision steps for the Newton's method""" 91 res = giant_steps(2, target) 92 if res[0] != 2: 93 res = [2] + res 94 return res 95 96def rs_trunc(p1, x, prec): 97 """ 98 Truncate the series in the ``x`` variable with precision ``prec``, 99 that is, modulo ``O(x**prec)`` 100 101 Examples 102 ======== 103 104 >>> from sympy.polys.domains import QQ 105 >>> from sympy.polys.rings import ring 106 >>> from sympy.polys.ring_series import rs_trunc 107 >>> R, x = ring('x', QQ) 108 >>> p = x**10 + x**5 + x + 1 109 >>> rs_trunc(p, x, 12) 110 x**10 + x**5 + x + 1 111 >>> rs_trunc(p, x, 10) 112 x**5 + x + 1 113 """ 114 R = p1.ring 115 p = R.zero 116 i = R.gens.index(x) 117 for exp1 in p1: 118 if exp1[i] >= prec: 119 continue 120 p[exp1] = p1[exp1] 121 return p 122 123def rs_is_puiseux(p, x): 124 """ 125 Test if ``p`` is Puiseux series in ``x``. 126 127 Raise an exception if it has a negative power in ``x``. 128 129 Examples 130 ======== 131 132 >>> from sympy.polys.domains import QQ 133 >>> from sympy.polys.rings import ring 134 >>> from sympy.polys.ring_series import rs_is_puiseux 135 >>> R, x = ring('x', QQ) 136 >>> p = x**QQ(2,5) + x**QQ(2,3) + x 137 >>> rs_is_puiseux(p, x) 138 True 139 """ 140 index = p.ring.gens.index(x) 141 for k in p: 142 if k[index] != int(k[index]): 143 return True 144 if k[index] < 0: 145 raise ValueError('The series is not regular in %s' % x) 146 return False 147 148def rs_puiseux(f, p, x, prec): 149 """ 150 Return the puiseux series for `f(p, x, prec)`. 151 152 To be used when function ``f`` is implemented only for regular series. 153 154 Examples 155 ======== 156 157 >>> from sympy.polys.domains import QQ 158 >>> from sympy.polys.rings import ring 159 >>> from sympy.polys.ring_series import rs_puiseux, rs_exp 160 >>> R, x = ring('x', QQ) 161 >>> p = x**QQ(2,5) + x**QQ(2,3) + x 162 >>> rs_puiseux(rs_exp,p, x, 1) 163 1/2*x**(4/5) + x**(2/3) + x**(2/5) + 1 164 """ 165 index = p.ring.gens.index(x) 166 n = 1 167 for k in p: 168 power = k[index] 169 if isinstance(power, Rational): 170 num, den = power.as_numer_denom() 171 n = int(n*den // igcd(n, den)) 172 elif power != int(power): 173 den = power.denominator 174 n = int(n*den // igcd(n, den)) 175 if n != 1: 176 p1 = pow_xin(p, index, n) 177 r = f(p1, x, prec*n) 178 n1 = QQ(1, n) 179 if isinstance(r, tuple): 180 r = tuple([pow_xin(rx, index, n1) for rx in r]) 181 else: 182 r = pow_xin(r, index, n1) 183 else: 184 r = f(p, x, prec) 185 return r 186 187def rs_puiseux2(f, p, q, x, prec): 188 """ 189 Return the puiseux series for `f(p, q, x, prec)`. 190 191 To be used when function ``f`` is implemented only for regular series. 192 """ 193 index = p.ring.gens.index(x) 194 n = 1 195 for k in p: 196 power = k[index] 197 if isinstance(power, Rational): 198 num, den = power.as_numer_denom() 199 n = n*den // igcd(n, den) 200 elif power != int(power): 201 den = power.denominator 202 n = n*den // igcd(n, den) 203 if n != 1: 204 p1 = pow_xin(p, index, n) 205 r = f(p1, q, x, prec*n) 206 n1 = QQ(1, n) 207 r = pow_xin(r, index, n1) 208 else: 209 r = f(p, q, x, prec) 210 return r 211 212def rs_mul(p1, p2, x, prec): 213 """ 214 Return the product of the given two series, modulo ``O(x**prec)``. 215 216 ``x`` is the series variable or its position in the generators. 217 218 Examples 219 ======== 220 221 >>> from sympy.polys.domains import QQ 222 >>> from sympy.polys.rings import ring 223 >>> from sympy.polys.ring_series import rs_mul 224 >>> R, x = ring('x', QQ) 225 >>> p1 = x**2 + 2*x + 1 226 >>> p2 = x + 1 227 >>> rs_mul(p1, p2, x, 3) 228 3*x**2 + 3*x + 1 229 """ 230 R = p1.ring 231 p = R.zero 232 if R.__class__ != p2.ring.__class__ or R != p2.ring: 233 raise ValueError('p1 and p2 must have the same ring') 234 iv = R.gens.index(x) 235 if not isinstance(p2, PolyElement): 236 raise ValueError('p1 and p2 must have the same ring') 237 if R == p2.ring: 238 get = p.get 239 items2 = list(p2.items()) 240 items2.sort(key=lambda e: e[0][iv]) 241 if R.ngens == 1: 242 for exp1, v1 in p1.items(): 243 for exp2, v2 in items2: 244 exp = exp1[0] + exp2[0] 245 if exp < prec: 246 exp = (exp, ) 247 p[exp] = get(exp, 0) + v1*v2 248 else: 249 break 250 else: 251 monomial_mul = R.monomial_mul 252 for exp1, v1 in p1.items(): 253 for exp2, v2 in items2: 254 if exp1[iv] + exp2[iv] < prec: 255 exp = monomial_mul(exp1, exp2) 256 p[exp] = get(exp, 0) + v1*v2 257 else: 258 break 259 260 p.strip_zero() 261 return p 262 263def rs_square(p1, x, prec): 264 """ 265 Square the series modulo ``O(x**prec)`` 266 267 Examples 268 ======== 269 270 >>> from sympy.polys.domains import QQ 271 >>> from sympy.polys.rings import ring 272 >>> from sympy.polys.ring_series import rs_square 273 >>> R, x = ring('x', QQ) 274 >>> p = x**2 + 2*x + 1 275 >>> rs_square(p, x, 3) 276 6*x**2 + 4*x + 1 277 """ 278 R = p1.ring 279 p = R.zero 280 iv = R.gens.index(x) 281 get = p.get 282 items = list(p1.items()) 283 items.sort(key=lambda e: e[0][iv]) 284 monomial_mul = R.monomial_mul 285 for i in range(len(items)): 286 exp1, v1 = items[i] 287 for j in range(i): 288 exp2, v2 = items[j] 289 if exp1[iv] + exp2[iv] < prec: 290 exp = monomial_mul(exp1, exp2) 291 p[exp] = get(exp, 0) + v1*v2 292 else: 293 break 294 p = p.imul_num(2) 295 get = p.get 296 for expv, v in p1.items(): 297 if 2*expv[iv] < prec: 298 e2 = monomial_mul(expv, expv) 299 p[e2] = get(e2, 0) + v**2 300 p.strip_zero() 301 return p 302 303def rs_pow(p1, n, x, prec): 304 """ 305 Return ``p1**n`` modulo ``O(x**prec)`` 306 307 Examples 308 ======== 309 310 >>> from sympy.polys.domains import QQ 311 >>> from sympy.polys.rings import ring 312 >>> from sympy.polys.ring_series import rs_pow 313 >>> R, x = ring('x', QQ) 314 >>> p = x + 1 315 >>> rs_pow(p, 4, x, 3) 316 6*x**2 + 4*x + 1 317 """ 318 R = p1.ring 319 if isinstance(n, Rational): 320 np = int(n.p) 321 nq = int(n.q) 322 if nq != 1: 323 res = rs_nth_root(p1, nq, x, prec) 324 if np != 1: 325 res = rs_pow(res, np, x, prec) 326 else: 327 res = rs_pow(p1, np, x, prec) 328 return res 329 330 n = as_int(n) 331 if n == 0: 332 if p1: 333 return R(1) 334 else: 335 raise ValueError('0**0 is undefined') 336 if n < 0: 337 p1 = rs_pow(p1, -n, x, prec) 338 return rs_series_inversion(p1, x, prec) 339 if n == 1: 340 return rs_trunc(p1, x, prec) 341 if n == 2: 342 return rs_square(p1, x, prec) 343 if n == 3: 344 p2 = rs_square(p1, x, prec) 345 return rs_mul(p1, p2, x, prec) 346 p = R(1) 347 while 1: 348 if n & 1: 349 p = rs_mul(p1, p, x, prec) 350 n -= 1 351 if not n: 352 break 353 p1 = rs_square(p1, x, prec) 354 n = n // 2 355 return p 356 357def rs_subs(p, rules, x, prec): 358 """ 359 Substitution with truncation according to the mapping in ``rules``. 360 361 Return a series with precision ``prec`` in the generator ``x`` 362 363 Note that substitutions are not done one after the other 364 365 >>> from sympy.polys.domains import QQ 366 >>> from sympy.polys.rings import ring 367 >>> from sympy.polys.ring_series import rs_subs 368 >>> R, x, y = ring('x, y', QQ) 369 >>> p = x**2 + y**2 370 >>> rs_subs(p, {x: x+ y, y: x+ 2*y}, x, 3) 371 2*x**2 + 6*x*y + 5*y**2 372 >>> (x + y)**2 + (x + 2*y)**2 373 2*x**2 + 6*x*y + 5*y**2 374 375 which differs from 376 377 >>> rs_subs(rs_subs(p, {x: x+ y}, x, 3), {y: x+ 2*y}, x, 3) 378 5*x**2 + 12*x*y + 8*y**2 379 380 Parameters 381 ---------- 382 p : :class:`~.PolyElement` Input series. 383 rules : ``dict`` with substitution mappings. 384 x : :class:`~.PolyElement` in which the series truncation is to be done. 385 prec : :class:`~.Integer` order of the series after truncation. 386 387 Examples 388 ======== 389 390 >>> from sympy.polys.domains import QQ 391 >>> from sympy.polys.rings import ring 392 >>> from sympy.polys.ring_series import rs_subs 393 >>> R, x, y = ring('x, y', QQ) 394 >>> rs_subs(x**2+y**2, {y: (x+y)**2}, x, 3) 395 6*x**2*y**2 + x**2 + 4*x*y**3 + y**4 396 """ 397 R = p.ring 398 ngens = R.ngens 399 d = R(0) 400 for i in range(ngens): 401 d[(i, 1)] = R.gens[i] 402 for var in rules: 403 d[(R.index(var), 1)] = rules[var] 404 p1 = R(0) 405 p_keys = sorted(p.keys()) 406 for expv in p_keys: 407 p2 = R(1) 408 for i in range(ngens): 409 power = expv[i] 410 if power == 0: 411 continue 412 if (i, power) not in d: 413 q, r = divmod(power, 2) 414 if r == 0 and (i, q) in d: 415 d[(i, power)] = rs_square(d[(i, q)], x, prec) 416 elif (i, power - 1) in d: 417 d[(i, power)] = rs_mul(d[(i, power - 1)], d[(i, 1)], 418 x, prec) 419 else: 420 d[(i, power)] = rs_pow(d[(i, 1)], power, x, prec) 421 p2 = rs_mul(p2, d[(i, power)], x, prec) 422 p1 += p2*p[expv] 423 return p1 424 425def _has_constant_term(p, x): 426 """ 427 Check if ``p`` has a constant term in ``x`` 428 429 Examples 430 ======== 431 432 >>> from sympy.polys.domains import QQ 433 >>> from sympy.polys.rings import ring 434 >>> from sympy.polys.ring_series import _has_constant_term 435 >>> R, x = ring('x', QQ) 436 >>> p = x**2 + x + 1 437 >>> _has_constant_term(p, x) 438 True 439 """ 440 R = p.ring 441 iv = R.gens.index(x) 442 zm = R.zero_monom 443 a = [0]*R.ngens 444 a[iv] = 1 445 miv = tuple(a) 446 for expv in p: 447 if monomial_min(expv, miv) == zm: 448 return True 449 return False 450 451def _get_constant_term(p, x): 452 """Return constant term in p with respect to x 453 454 Note that it is not simply `p[R.zero_monom]` as there might be multiple 455 generators in the ring R. We want the `x`-free term which can contain other 456 generators. 457 """ 458 R = p.ring 459 i = R.gens.index(x) 460 zm = R.zero_monom 461 a = [0]*R.ngens 462 a[i] = 1 463 miv = tuple(a) 464 c = 0 465 for expv in p: 466 if monomial_min(expv, miv) == zm: 467 c += R({expv: p[expv]}) 468 return c 469 470def _check_series_var(p, x, name): 471 index = p.ring.gens.index(x) 472 m = min(p, key=lambda k: k[index])[index] 473 if m < 0: 474 raise PoleError("Asymptotic expansion of %s around [oo] not " 475 "implemented." % name) 476 return index, m 477 478def _series_inversion1(p, x, prec): 479 """ 480 Univariate series inversion ``1/p`` modulo ``O(x**prec)``. 481 482 The Newton method is used. 483 484 Examples 485 ======== 486 487 >>> from sympy.polys.domains import QQ 488 >>> from sympy.polys.rings import ring 489 >>> from sympy.polys.ring_series import _series_inversion1 490 >>> R, x = ring('x', QQ) 491 >>> p = x + 1 492 >>> _series_inversion1(p, x, 4) 493 -x**3 + x**2 - x + 1 494 """ 495 if rs_is_puiseux(p, x): 496 return rs_puiseux(_series_inversion1, p, x, prec) 497 R = p.ring 498 zm = R.zero_monom 499 c = p[zm] 500 501 # giant_steps does not seem to work with PythonRational numbers with 1 as 502 # denominator. This makes sure such a number is converted to integer. 503 if prec == int(prec): 504 prec = int(prec) 505 506 if zm not in p: 507 raise ValueError("No constant term in series") 508 if _has_constant_term(p - c, x): 509 raise ValueError("p cannot contain a constant term depending on " 510 "parameters") 511 one = R(1) 512 if R.domain is EX: 513 one = 1 514 if c != one: 515 # TODO add check that it is a unit 516 p1 = R(1)/c 517 else: 518 p1 = R(1) 519 for precx in _giant_steps(prec): 520 t = 1 - rs_mul(p1, p, x, precx) 521 p1 = p1 + rs_mul(p1, t, x, precx) 522 return p1 523 524def rs_series_inversion(p, x, prec): 525 """ 526 Multivariate series inversion ``1/p`` modulo ``O(x**prec)``. 527 528 Examples 529 ======== 530 531 >>> from sympy.polys.domains import QQ 532 >>> from sympy.polys.rings import ring 533 >>> from sympy.polys.ring_series import rs_series_inversion 534 >>> R, x, y = ring('x, y', QQ) 535 >>> rs_series_inversion(1 + x*y**2, x, 4) 536 -x**3*y**6 + x**2*y**4 - x*y**2 + 1 537 >>> rs_series_inversion(1 + x*y**2, y, 4) 538 -x*y**2 + 1 539 >>> rs_series_inversion(x + x**2, x, 4) 540 x**3 - x**2 + x - 1 + x**(-1) 541 """ 542 R = p.ring 543 if p == R.zero: 544 raise ZeroDivisionError 545 zm = R.zero_monom 546 index = R.gens.index(x) 547 m = min(p, key=lambda k: k[index])[index] 548 if m: 549 p = mul_xin(p, index, -m) 550 prec = prec + m 551 if zm not in p: 552 raise NotImplementedError("No constant term in series") 553 554 if _has_constant_term(p - p[zm], x): 555 raise NotImplementedError("p - p[0] must not have a constant term in " 556 "the series variables") 557 r = _series_inversion1(p, x, prec) 558 if m != 0: 559 r = mul_xin(r, index, -m) 560 return r 561 562def _coefficient_t(p, t): 563 r"""Coefficient of `x_i**j` in p, where ``t`` = (i, j)""" 564 i, j = t 565 R = p.ring 566 expv1 = [0]*R.ngens 567 expv1[i] = j 568 expv1 = tuple(expv1) 569 p1 = R(0) 570 for expv in p: 571 if expv[i] == j: 572 p1[monomial_div(expv, expv1)] = p[expv] 573 return p1 574 575def rs_series_reversion(p, x, n, y): 576 r""" 577 Reversion of a series. 578 579 ``p`` is a series with ``O(x**n)`` of the form $p = ax + f(x)$ 580 where $a$ is a number different from 0. 581 582 $f(x) = \sum_{k=2}^{n-1} a_kx_k$ 583 584 Parameters 585 ========== 586 587 a_k : Can depend polynomially on other variables, not indicated. 588 x : Variable with name x. 589 y : Variable with name y. 590 591 Returns 592 ======= 593 594 Solve $p = y$, that is, given $ax + f(x) - y = 0$, 595 find the solution $x = r(y)$ up to $O(y^n)$. 596 597 Algorithm 598 ========= 599 600 If $r_i$ is the solution at order $i$, then: 601 $ar_i + f(r_i) - y = O\left(y^{i + 1}\right)$ 602 603 and if $r_{i + 1}$ is the solution at order $i + 1$, then: 604 $ar_{i + 1} + f(r_{i + 1}) - y = O\left(y^{i + 2}\right)$ 605 606 We have, $r_{i + 1} = r_i + e$, such that, 607 $ae + f(r_i) = O\left(y^{i + 2}\right)$ 608 or $e = -f(r_i)/a$ 609 610 So we use the recursion relation: 611 $r_{i + 1} = r_i - f(r_i)/a$ 612 with the boundary condition: $r_1 = y$ 613 614 Examples 615 ======== 616 617 >>> from sympy.polys.domains import QQ 618 >>> from sympy.polys.rings import ring 619 >>> from sympy.polys.ring_series import rs_series_reversion, rs_trunc 620 >>> R, x, y, a, b = ring('x, y, a, b', QQ) 621 >>> p = x - x**2 - 2*b*x**2 + 2*a*b*x**2 622 >>> p1 = rs_series_reversion(p, x, 3, y); p1 623 -2*y**2*a*b + 2*y**2*b + y**2 + y 624 >>> rs_trunc(p.compose(x, p1), y, 3) 625 y 626 """ 627 if rs_is_puiseux(p, x): 628 raise NotImplementedError 629 R = p.ring 630 nx = R.gens.index(x) 631 y = R(y) 632 ny = R.gens.index(y) 633 if _has_constant_term(p, x): 634 raise ValueError("p must not contain a constant term in the series " 635 "variable") 636 a = _coefficient_t(p, (nx, 1)) 637 zm = R.zero_monom 638 assert zm in a and len(a) == 1 639 a = a[zm] 640 r = y/a 641 for i in range(2, n): 642 sp = rs_subs(p, {x: r}, y, i + 1) 643 sp = _coefficient_t(sp, (ny, i))*y**i 644 r -= sp/a 645 return r 646 647def rs_series_from_list(p, c, x, prec, concur=1): 648 """ 649 Return a series `sum c[n]*p**n` modulo `O(x**prec)`. 650 651 It reduces the number of multiplications by summing concurrently. 652 653 `ax = [1, p, p**2, .., p**(J - 1)]` 654 `s = sum(c[i]*ax[i]` for i in `range(r, (r + 1)*J))*p**((K - 1)*J)` 655 with `K >= (n + 1)/J` 656 657 Examples 658 ======== 659 660 >>> from sympy.polys.domains import QQ 661 >>> from sympy.polys.rings import ring 662 >>> from sympy.polys.ring_series import rs_series_from_list, rs_trunc 663 >>> R, x = ring('x', QQ) 664 >>> p = x**2 + x + 1 665 >>> c = [1, 2, 3] 666 >>> rs_series_from_list(p, c, x, 4) 667 6*x**3 + 11*x**2 + 8*x + 6 668 >>> rs_trunc(1 + 2*p + 3*p**2, x, 4) 669 6*x**3 + 11*x**2 + 8*x + 6 670 >>> pc = R.from_list(list(reversed(c))) 671 >>> rs_trunc(pc.compose(x, p), x, 4) 672 6*x**3 + 11*x**2 + 8*x + 6 673 674 """ 675 676 # TODO: Add this when it is documented in Sphinx 677 """ 678 See Also 679 ======== 680 681 sympy.polys.rings.PolyRing.compose 682 683 """ 684 R = p.ring 685 n = len(c) 686 if not concur: 687 q = R(1) 688 s = c[0]*q 689 for i in range(1, n): 690 q = rs_mul(q, p, x, prec) 691 s += c[i]*q 692 return s 693 J = int(math.sqrt(n) + 1) 694 K, r = divmod(n, J) 695 if r: 696 K += 1 697 ax = [R(1)] 698 q = R(1) 699 if len(p) < 20: 700 for i in range(1, J): 701 q = rs_mul(q, p, x, prec) 702 ax.append(q) 703 else: 704 for i in range(1, J): 705 if i % 2 == 0: 706 q = rs_square(ax[i//2], x, prec) 707 else: 708 q = rs_mul(q, p, x, prec) 709 ax.append(q) 710 # optimize using rs_square 711 pj = rs_mul(ax[-1], p, x, prec) 712 b = R(1) 713 s = R(0) 714 for k in range(K - 1): 715 r = J*k 716 s1 = c[r] 717 for j in range(1, J): 718 s1 += c[r + j]*ax[j] 719 s1 = rs_mul(s1, b, x, prec) 720 s += s1 721 b = rs_mul(b, pj, x, prec) 722 if not b: 723 break 724 k = K - 1 725 r = J*k 726 if r < n: 727 s1 = c[r]*R(1) 728 for j in range(1, J): 729 if r + j >= n: 730 break 731 s1 += c[r + j]*ax[j] 732 s1 = rs_mul(s1, b, x, prec) 733 s += s1 734 return s 735 736def rs_diff(p, x): 737 """ 738 Return partial derivative of ``p`` with respect to ``x``. 739 740 Parameters 741 ========== 742 743 x : :class:`~.PolyElement` with respect to which ``p`` is differentiated. 744 745 Examples 746 ======== 747 748 >>> from sympy.polys.domains import QQ 749 >>> from sympy.polys.rings import ring 750 >>> from sympy.polys.ring_series import rs_diff 751 >>> R, x, y = ring('x, y', QQ) 752 >>> p = x + x**2*y**3 753 >>> rs_diff(p, x) 754 2*x*y**3 + 1 755 """ 756 R = p.ring 757 n = R.gens.index(x) 758 p1 = R.zero 759 mn = [0]*R.ngens 760 mn[n] = 1 761 mn = tuple(mn) 762 for expv in p: 763 if expv[n]: 764 e = monomial_ldiv(expv, mn) 765 p1[e] = R.domain_new(p[expv]*expv[n]) 766 return p1 767 768def rs_integrate(p, x): 769 """ 770 Integrate ``p`` with respect to ``x``. 771 772 Parameters 773 ========== 774 775 x : :class:`~.PolyElement` with respect to which ``p`` is integrated. 776 777 Examples 778 ======== 779 780 >>> from sympy.polys.domains import QQ 781 >>> from sympy.polys.rings import ring 782 >>> from sympy.polys.ring_series import rs_integrate 783 >>> R, x, y = ring('x, y', QQ) 784 >>> p = x + x**2*y**3 785 >>> rs_integrate(p, x) 786 1/3*x**3*y**3 + 1/2*x**2 787 """ 788 R = p.ring 789 p1 = R.zero 790 n = R.gens.index(x) 791 mn = [0]*R.ngens 792 mn[n] = 1 793 mn = tuple(mn) 794 795 for expv in p: 796 e = monomial_mul(expv, mn) 797 p1[e] = R.domain_new(p[expv]/(expv[n] + 1)) 798 return p1 799 800def rs_fun(p, f, *args): 801 r""" 802 Function of a multivariate series computed by substitution. 803 804 The case with f method name is used to compute `rs\_tan` and `rs\_nth\_root` 805 of a multivariate series: 806 807 `rs\_fun(p, tan, iv, prec)` 808 809 tan series is first computed for a dummy variable _x, 810 i.e, `rs\_tan(\_x, iv, prec)`. Then we substitute _x with p to get the 811 desired series 812 813 Parameters 814 ========== 815 816 p : :class:`~.PolyElement` The multivariate series to be expanded. 817 f : `ring\_series` function to be applied on `p`. 818 args[-2] : :class:`~.PolyElement` with respect to which, the series is to be expanded. 819 args[-1] : Required order of the expanded series. 820 821 Examples 822 ======== 823 824 >>> from sympy.polys.domains import QQ 825 >>> from sympy.polys.rings import ring 826 >>> from sympy.polys.ring_series import rs_fun, _tan1 827 >>> R, x, y = ring('x, y', QQ) 828 >>> p = x + x*y + x**2*y + x**3*y**2 829 >>> rs_fun(p, _tan1, x, 4) 830 1/3*x**3*y**3 + 2*x**3*y**2 + x**3*y + 1/3*x**3 + x**2*y + x*y + x 831 """ 832 _R = p.ring 833 R1, _x = ring('_x', _R.domain) 834 h = int(args[-1]) 835 args1 = args[:-2] + (_x, h) 836 zm = _R.zero_monom 837 # separate the constant term of the series 838 # compute the univariate series f(_x, .., 'x', sum(nv)) 839 if zm in p: 840 x1 = _x + p[zm] 841 p1 = p - p[zm] 842 else: 843 x1 = _x 844 p1 = p 845 if isinstance(f, str): 846 q = getattr(x1, f)(*args1) 847 else: 848 q = f(x1, *args1) 849 a = sorted(q.items()) 850 c = [0]*h 851 for x in a: 852 c[x[0][0]] = x[1] 853 p1 = rs_series_from_list(p1, c, args[-2], args[-1]) 854 return p1 855 856def mul_xin(p, i, n): 857 r""" 858 Return `p*x_i**n`. 859 860 `x\_i` is the ith variable in ``p``. 861 """ 862 R = p.ring 863 q = R(0) 864 for k, v in p.items(): 865 k1 = list(k) 866 k1[i] += n 867 q[tuple(k1)] = v 868 return q 869 870def pow_xin(p, i, n): 871 """ 872 >>> from sympy.polys.domains import QQ 873 >>> from sympy.polys.rings import ring 874 >>> from sympy.polys.ring_series import pow_xin 875 >>> R, x, y = ring('x, y', QQ) 876 >>> p = x**QQ(2,5) + x + x**QQ(2,3) 877 >>> index = p.ring.gens.index(x) 878 >>> pow_xin(p, index, 15) 879 x**15 + x**10 + x**6 880 """ 881 R = p.ring 882 q = R(0) 883 for k, v in p.items(): 884 k1 = list(k) 885 k1[i] *= n 886 q[tuple(k1)] = v 887 return q 888 889def _nth_root1(p, n, x, prec): 890 """ 891 Univariate series expansion of the nth root of ``p``. 892 893 The Newton method is used. 894 """ 895 if rs_is_puiseux(p, x): 896 return rs_puiseux2(_nth_root1, p, n, x, prec) 897 R = p.ring 898 zm = R.zero_monom 899 if zm not in p: 900 raise NotImplementedError('No constant term in series') 901 n = as_int(n) 902 assert p[zm] == 1 903 p1 = R(1) 904 if p == 1: 905 return p 906 if n == 0: 907 return R(1) 908 if n == 1: 909 return p 910 if n < 0: 911 n = -n 912 sign = 1 913 else: 914 sign = 0 915 for precx in _giant_steps(prec): 916 tmp = rs_pow(p1, n + 1, x, precx) 917 tmp = rs_mul(tmp, p, x, precx) 918 p1 += p1/n - tmp/n 919 if sign: 920 return p1 921 else: 922 return _series_inversion1(p1, x, prec) 923 924def rs_nth_root(p, n, x, prec): 925 """ 926 Multivariate series expansion of the nth root of ``p``. 927 928 Parameters 929 ========== 930 931 p : Expr 932 The polynomial to computer the root of. 933 n : integer 934 The order of the root to be computed. 935 x : :class:`~.PolyElement` 936 prec : integer 937 Order of the expanded series. 938 939 Notes 940 ===== 941 942 The result of this function is dependent on the ring over which the 943 polynomial has been defined. If the answer involves a root of a constant, 944 make sure that the polynomial is over a real field. It can not yet handle 945 roots of symbols. 946 947 Examples 948 ======== 949 950 >>> from sympy.polys.domains import QQ, RR 951 >>> from sympy.polys.rings import ring 952 >>> from sympy.polys.ring_series import rs_nth_root 953 >>> R, x, y = ring('x, y', QQ) 954 >>> rs_nth_root(1 + x + x*y, -3, x, 3) 955 2/9*x**2*y**2 + 4/9*x**2*y + 2/9*x**2 - 1/3*x*y - 1/3*x + 1 956 >>> R, x, y = ring('x, y', RR) 957 >>> rs_nth_root(3 + x + x*y, 3, x, 2) 958 0.160249952256379*x*y + 0.160249952256379*x + 1.44224957030741 959 """ 960 if n == 0: 961 if p == 0: 962 raise ValueError('0**0 expression') 963 else: 964 return p.ring(1) 965 if n == 1: 966 return rs_trunc(p, x, prec) 967 R = p.ring 968 index = R.gens.index(x) 969 m = min(p, key=lambda k: k[index])[index] 970 p = mul_xin(p, index, -m) 971 prec -= m 972 973 if _has_constant_term(p - 1, x): 974 zm = R.zero_monom 975 c = p[zm] 976 if R.domain is EX: 977 c_expr = c.as_expr() 978 const = c_expr**QQ(1, n) 979 elif isinstance(c, PolyElement): 980 try: 981 c_expr = c.as_expr() 982 const = R(c_expr**(QQ(1, n))) 983 except ValueError: 984 raise DomainError("The given series can't be expanded in " 985 "this domain.") 986 else: 987 try: # RealElement doesn't support 988 const = R(c**Rational(1, n)) # exponentiation with mpq object 989 except ValueError: # as exponent 990 raise DomainError("The given series can't be expanded in " 991 "this domain.") 992 res = rs_nth_root(p/c, n, x, prec)*const 993 else: 994 res = _nth_root1(p, n, x, prec) 995 if m: 996 m = QQ(m, n) 997 res = mul_xin(res, index, m) 998 return res 999 1000def rs_log(p, x, prec): 1001 """ 1002 The Logarithm of ``p`` modulo ``O(x**prec)``. 1003 1004 Notes 1005 ===== 1006 1007 Truncation of ``integral dx p**-1*d p/dx`` is used. 1008 1009 Examples 1010 ======== 1011 1012 >>> from sympy.polys.domains import QQ 1013 >>> from sympy.polys.rings import ring 1014 >>> from sympy.polys.ring_series import rs_log 1015 >>> R, x = ring('x', QQ) 1016 >>> rs_log(1 + x, x, 8) 1017 1/7*x**7 - 1/6*x**6 + 1/5*x**5 - 1/4*x**4 + 1/3*x**3 - 1/2*x**2 + x 1018 >>> rs_log(x**QQ(3, 2) + 1, x, 5) 1019 1/3*x**(9/2) - 1/2*x**3 + x**(3/2) 1020 """ 1021 if rs_is_puiseux(p, x): 1022 return rs_puiseux(rs_log, p, x, prec) 1023 R = p.ring 1024 if p == 1: 1025 return R.zero 1026 c = _get_constant_term(p, x) 1027 if c: 1028 const = 0 1029 if c == 1: 1030 pass 1031 else: 1032 c_expr = c.as_expr() 1033 if R.domain is EX: 1034 const = log(c_expr) 1035 elif isinstance(c, PolyElement): 1036 try: 1037 const = R(log(c_expr)) 1038 except ValueError: 1039 R = R.add_gens([log(c_expr)]) 1040 p = p.set_ring(R) 1041 x = x.set_ring(R) 1042 c = c.set_ring(R) 1043 const = R(log(c_expr)) 1044 else: 1045 try: 1046 const = R(log(c)) 1047 except ValueError: 1048 raise DomainError("The given series can't be expanded in " 1049 "this domain.") 1050 1051 dlog = p.diff(x) 1052 dlog = rs_mul(dlog, _series_inversion1(p, x, prec), x, prec - 1) 1053 return rs_integrate(dlog, x) + const 1054 else: 1055 raise NotImplementedError 1056 1057def rs_LambertW(p, x, prec): 1058 """ 1059 Calculate the series expansion of the principal branch of the Lambert W 1060 function. 1061 1062 Examples 1063 ======== 1064 1065 >>> from sympy.polys.domains import QQ 1066 >>> from sympy.polys.rings import ring 1067 >>> from sympy.polys.ring_series import rs_LambertW 1068 >>> R, x, y = ring('x, y', QQ) 1069 >>> rs_LambertW(x + x*y, x, 3) 1070 -x**2*y**2 - 2*x**2*y - x**2 + x*y + x 1071 1072 See Also 1073 ======== 1074 1075 LambertW 1076 """ 1077 if rs_is_puiseux(p, x): 1078 return rs_puiseux(rs_LambertW, p, x, prec) 1079 R = p.ring 1080 p1 = R(0) 1081 if _has_constant_term(p, x): 1082 raise NotImplementedError("Polynomial must not have constant term in " 1083 "the series variables") 1084 if x in R.gens: 1085 for precx in _giant_steps(prec): 1086 e = rs_exp(p1, x, precx) 1087 p2 = rs_mul(e, p1, x, precx) - p 1088 p3 = rs_mul(e, p1 + 1, x, precx) 1089 p3 = rs_series_inversion(p3, x, precx) 1090 tmp = rs_mul(p2, p3, x, precx) 1091 p1 -= tmp 1092 return p1 1093 else: 1094 raise NotImplementedError 1095 1096def _exp1(p, x, prec): 1097 r"""Helper function for `rs\_exp`. """ 1098 R = p.ring 1099 p1 = R(1) 1100 for precx in _giant_steps(prec): 1101 pt = p - rs_log(p1, x, precx) 1102 tmp = rs_mul(pt, p1, x, precx) 1103 p1 += tmp 1104 return p1 1105 1106def rs_exp(p, x, prec): 1107 """ 1108 Exponentiation of a series modulo ``O(x**prec)`` 1109 1110 Examples 1111 ======== 1112 1113 >>> from sympy.polys.domains import QQ 1114 >>> from sympy.polys.rings import ring 1115 >>> from sympy.polys.ring_series import rs_exp 1116 >>> R, x = ring('x', QQ) 1117 >>> rs_exp(x**2, x, 7) 1118 1/6*x**6 + 1/2*x**4 + x**2 + 1 1119 """ 1120 if rs_is_puiseux(p, x): 1121 return rs_puiseux(rs_exp, p, x, prec) 1122 R = p.ring 1123 c = _get_constant_term(p, x) 1124 if c: 1125 if R.domain is EX: 1126 c_expr = c.as_expr() 1127 const = exp(c_expr) 1128 elif isinstance(c, PolyElement): 1129 try: 1130 c_expr = c.as_expr() 1131 const = R(exp(c_expr)) 1132 except ValueError: 1133 R = R.add_gens([exp(c_expr)]) 1134 p = p.set_ring(R) 1135 x = x.set_ring(R) 1136 c = c.set_ring(R) 1137 const = R(exp(c_expr)) 1138 else: 1139 try: 1140 const = R(exp(c)) 1141 except ValueError: 1142 raise DomainError("The given series can't be expanded in " 1143 "this domain.") 1144 p1 = p - c 1145 1146 # Makes use of sympy functions to evaluate the values of the cos/sin 1147 # of the constant term. 1148 return const*rs_exp(p1, x, prec) 1149 1150 if len(p) > 20: 1151 return _exp1(p, x, prec) 1152 one = R(1) 1153 n = 1 1154 c = [] 1155 for k in range(prec): 1156 c.append(one/n) 1157 k += 1 1158 n *= k 1159 1160 r = rs_series_from_list(p, c, x, prec) 1161 return r 1162 1163def _atan(p, iv, prec): 1164 """ 1165 Expansion using formula. 1166 1167 Faster on very small and univariate series. 1168 """ 1169 R = p.ring 1170 mo = R(-1) 1171 c = [-mo] 1172 p2 = rs_square(p, iv, prec) 1173 for k in range(1, prec): 1174 c.append(mo**k/(2*k + 1)) 1175 s = rs_series_from_list(p2, c, iv, prec) 1176 s = rs_mul(s, p, iv, prec) 1177 return s 1178 1179def rs_atan(p, x, prec): 1180 """ 1181 The arctangent of a series 1182 1183 Return the series expansion of the atan of ``p``, about 0. 1184 1185 Examples 1186 ======== 1187 1188 >>> from sympy.polys.domains import QQ 1189 >>> from sympy.polys.rings import ring 1190 >>> from sympy.polys.ring_series import rs_atan 1191 >>> R, x, y = ring('x, y', QQ) 1192 >>> rs_atan(x + x*y, x, 4) 1193 -1/3*x**3*y**3 - x**3*y**2 - x**3*y - 1/3*x**3 + x*y + x 1194 1195 See Also 1196 ======== 1197 1198 atan 1199 """ 1200 if rs_is_puiseux(p, x): 1201 return rs_puiseux(rs_atan, p, x, prec) 1202 R = p.ring 1203 const = 0 1204 if _has_constant_term(p, x): 1205 zm = R.zero_monom 1206 c = p[zm] 1207 if R.domain is EX: 1208 c_expr = c.as_expr() 1209 const = atan(c_expr) 1210 elif isinstance(c, PolyElement): 1211 try: 1212 c_expr = c.as_expr() 1213 const = R(atan(c_expr)) 1214 except ValueError: 1215 raise DomainError("The given series can't be expanded in " 1216 "this domain.") 1217 else: 1218 try: 1219 const = R(atan(c)) 1220 except ValueError: 1221 raise DomainError("The given series can't be expanded in " 1222 "this domain.") 1223 1224 # Instead of using a closed form formula, we differentiate atan(p) to get 1225 # `1/(1+p**2) * dp`, whose series expansion is much easier to calculate. 1226 # Finally we integrate to get back atan 1227 dp = p.diff(x) 1228 p1 = rs_square(p, x, prec) + R(1) 1229 p1 = rs_series_inversion(p1, x, prec - 1) 1230 p1 = rs_mul(dp, p1, x, prec - 1) 1231 return rs_integrate(p1, x) + const 1232 1233def rs_asin(p, x, prec): 1234 """ 1235 Arcsine of a series 1236 1237 Return the series expansion of the asin of ``p``, about 0. 1238 1239 Examples 1240 ======== 1241 1242 >>> from sympy.polys.domains import QQ 1243 >>> from sympy.polys.rings import ring 1244 >>> from sympy.polys.ring_series import rs_asin 1245 >>> R, x, y = ring('x, y', QQ) 1246 >>> rs_asin(x, x, 8) 1247 5/112*x**7 + 3/40*x**5 + 1/6*x**3 + x 1248 1249 See Also 1250 ======== 1251 1252 asin 1253 """ 1254 if rs_is_puiseux(p, x): 1255 return rs_puiseux(rs_asin, p, x, prec) 1256 if _has_constant_term(p, x): 1257 raise NotImplementedError("Polynomial must not have constant term in " 1258 "series variables") 1259 R = p.ring 1260 if x in R.gens: 1261 # get a good value 1262 if len(p) > 20: 1263 dp = rs_diff(p, x) 1264 p1 = 1 - rs_square(p, x, prec - 1) 1265 p1 = rs_nth_root(p1, -2, x, prec - 1) 1266 p1 = rs_mul(dp, p1, x, prec - 1) 1267 return rs_integrate(p1, x) 1268 one = R(1) 1269 c = [0, one, 0] 1270 for k in range(3, prec, 2): 1271 c.append((k - 2)**2*c[-2]/(k*(k - 1))) 1272 c.append(0) 1273 return rs_series_from_list(p, c, x, prec) 1274 1275 else: 1276 raise NotImplementedError 1277 1278def _tan1(p, x, prec): 1279 r""" 1280 Helper function of :func:`rs_tan`. 1281 1282 Return the series expansion of tan of a univariate series using Newton's 1283 method. It takes advantage of the fact that series expansion of atan is 1284 easier than that of tan. 1285 1286 Consider `f(x) = y - \arctan(x)` 1287 Let r be a root of f(x) found using Newton's method. 1288 Then `f(r) = 0` 1289 Or `y = \arctan(x)` where `x = \tan(y)` as required. 1290 """ 1291 R = p.ring 1292 p1 = R(0) 1293 for precx in _giant_steps(prec): 1294 tmp = p - rs_atan(p1, x, precx) 1295 tmp = rs_mul(tmp, 1 + rs_square(p1, x, precx), x, precx) 1296 p1 += tmp 1297 return p1 1298 1299def rs_tan(p, x, prec): 1300 """ 1301 Tangent of a series. 1302 1303 Return the series expansion of the tan of ``p``, about 0. 1304 1305 Examples 1306 ======== 1307 1308 >>> from sympy.polys.domains import QQ 1309 >>> from sympy.polys.rings import ring 1310 >>> from sympy.polys.ring_series import rs_tan 1311 >>> R, x, y = ring('x, y', QQ) 1312 >>> rs_tan(x + x*y, x, 4) 1313 1/3*x**3*y**3 + x**3*y**2 + x**3*y + 1/3*x**3 + x*y + x 1314 1315 See Also 1316 ======== 1317 1318 _tan1, tan 1319 """ 1320 if rs_is_puiseux(p, x): 1321 r = rs_puiseux(rs_tan, p, x, prec) 1322 return r 1323 R = p.ring 1324 const = 0 1325 c = _get_constant_term(p, x) 1326 if c: 1327 if R.domain is EX: 1328 c_expr = c.as_expr() 1329 const = tan(c_expr) 1330 elif isinstance(c, PolyElement): 1331 try: 1332 c_expr = c.as_expr() 1333 const = R(tan(c_expr)) 1334 except ValueError: 1335 R = R.add_gens([tan(c_expr, )]) 1336 p = p.set_ring(R) 1337 x = x.set_ring(R) 1338 c = c.set_ring(R) 1339 const = R(tan(c_expr)) 1340 else: 1341 try: 1342 const = R(tan(c)) 1343 except ValueError: 1344 raise DomainError("The given series can't be expanded in " 1345 "this domain.") 1346 p1 = p - c 1347 1348 # Makes use of sympy functions to evaluate the values of the cos/sin 1349 # of the constant term. 1350 t2 = rs_tan(p1, x, prec) 1351 t = rs_series_inversion(1 - const*t2, x, prec) 1352 return rs_mul(const + t2, t, x, prec) 1353 1354 if R.ngens == 1: 1355 return _tan1(p, x, prec) 1356 else: 1357 return rs_fun(p, rs_tan, x, prec) 1358 1359def rs_cot(p, x, prec): 1360 """ 1361 Cotangent of a series 1362 1363 Return the series expansion of the cot of ``p``, about 0. 1364 1365 Examples 1366 ======== 1367 1368 >>> from sympy.polys.domains import QQ 1369 >>> from sympy.polys.rings import ring 1370 >>> from sympy.polys.ring_series import rs_cot 1371 >>> R, x, y = ring('x, y', QQ) 1372 >>> rs_cot(x, x, 6) 1373 -2/945*x**5 - 1/45*x**3 - 1/3*x + x**(-1) 1374 1375 See Also 1376 ======== 1377 1378 cot 1379 """ 1380 # It can not handle series like `p = x + x*y` where the coefficient of the 1381 # linear term in the series variable is symbolic. 1382 if rs_is_puiseux(p, x): 1383 r = rs_puiseux(rs_cot, p, x, prec) 1384 return r 1385 i, m = _check_series_var(p, x, 'cot') 1386 prec1 = prec + 2*m 1387 c, s = rs_cos_sin(p, x, prec1) 1388 s = mul_xin(s, i, -m) 1389 s = rs_series_inversion(s, x, prec1) 1390 res = rs_mul(c, s, x, prec1) 1391 res = mul_xin(res, i, -m) 1392 res = rs_trunc(res, x, prec) 1393 return res 1394 1395def rs_sin(p, x, prec): 1396 """ 1397 Sine of a series 1398 1399 Return the series expansion of the sin of ``p``, about 0. 1400 1401 Examples 1402 ======== 1403 1404 >>> from sympy.polys.domains import QQ 1405 >>> from sympy.polys.rings import ring 1406 >>> from sympy.polys.ring_series import rs_sin 1407 >>> R, x, y = ring('x, y', QQ) 1408 >>> rs_sin(x + x*y, x, 4) 1409 -1/6*x**3*y**3 - 1/2*x**3*y**2 - 1/2*x**3*y - 1/6*x**3 + x*y + x 1410 >>> rs_sin(x**QQ(3, 2) + x*y**QQ(7, 5), x, 4) 1411 -1/2*x**(7/2)*y**(14/5) - 1/6*x**3*y**(21/5) + x**(3/2) + x*y**(7/5) 1412 1413 See Also 1414 ======== 1415 1416 sin 1417 """ 1418 if rs_is_puiseux(p, x): 1419 return rs_puiseux(rs_sin, p, x, prec) 1420 R = x.ring 1421 if not p: 1422 return R(0) 1423 c = _get_constant_term(p, x) 1424 if c: 1425 if R.domain is EX: 1426 c_expr = c.as_expr() 1427 t1, t2 = sin(c_expr), cos(c_expr) 1428 elif isinstance(c, PolyElement): 1429 try: 1430 c_expr = c.as_expr() 1431 t1, t2 = R(sin(c_expr)), R(cos(c_expr)) 1432 except ValueError: 1433 R = R.add_gens([sin(c_expr), cos(c_expr)]) 1434 p = p.set_ring(R) 1435 x = x.set_ring(R) 1436 c = c.set_ring(R) 1437 t1, t2 = R(sin(c_expr)), R(cos(c_expr)) 1438 else: 1439 try: 1440 t1, t2 = R(sin(c)), R(cos(c)) 1441 except ValueError: 1442 raise DomainError("The given series can't be expanded in " 1443 "this domain.") 1444 p1 = p - c 1445 1446 # Makes use of sympy cos, sin functions to evaluate the values of the 1447 # cos/sin of the constant term. 1448 return rs_sin(p1, x, prec)*t2 + rs_cos(p1, x, prec)*t1 1449 1450 # Series is calculated in terms of tan as its evaluation is fast. 1451 if len(p) > 20 and R.ngens == 1: 1452 t = rs_tan(p/2, x, prec) 1453 t2 = rs_square(t, x, prec) 1454 p1 = rs_series_inversion(1 + t2, x, prec) 1455 return rs_mul(p1, 2*t, x, prec) 1456 one = R(1) 1457 n = 1 1458 c = [0] 1459 for k in range(2, prec + 2, 2): 1460 c.append(one/n) 1461 c.append(0) 1462 n *= -k*(k + 1) 1463 return rs_series_from_list(p, c, x, prec) 1464 1465def rs_cos(p, x, prec): 1466 """ 1467 Cosine of a series 1468 1469 Return the series expansion of the cos of ``p``, about 0. 1470 1471 Examples 1472 ======== 1473 1474 >>> from sympy.polys.domains import QQ 1475 >>> from sympy.polys.rings import ring 1476 >>> from sympy.polys.ring_series import rs_cos 1477 >>> R, x, y = ring('x, y', QQ) 1478 >>> rs_cos(x + x*y, x, 4) 1479 -1/2*x**2*y**2 - x**2*y - 1/2*x**2 + 1 1480 >>> rs_cos(x + x*y, x, 4)/x**QQ(7, 5) 1481 -1/2*x**(3/5)*y**2 - x**(3/5)*y - 1/2*x**(3/5) + x**(-7/5) 1482 1483 See Also 1484 ======== 1485 1486 cos 1487 """ 1488 if rs_is_puiseux(p, x): 1489 return rs_puiseux(rs_cos, p, x, prec) 1490 R = p.ring 1491 c = _get_constant_term(p, x) 1492 if c: 1493 if R.domain is EX: 1494 c_expr = c.as_expr() 1495 _, _ = sin(c_expr), cos(c_expr) 1496 elif isinstance(c, PolyElement): 1497 try: 1498 c_expr = c.as_expr() 1499 _, _ = R(sin(c_expr)), R(cos(c_expr)) 1500 except ValueError: 1501 R = R.add_gens([sin(c_expr), cos(c_expr)]) 1502 p = p.set_ring(R) 1503 x = x.set_ring(R) 1504 c = c.set_ring(R) 1505 else: 1506 try: 1507 _, _ = R(sin(c)), R(cos(c)) 1508 except ValueError: 1509 raise DomainError("The given series can't be expanded in " 1510 "this domain.") 1511 p1 = p - c 1512 1513 # Makes use of sympy cos, sin functions to evaluate the values of the 1514 # cos/sin of the constant term. 1515 p_cos = rs_cos(p1, x, prec) 1516 p_sin = rs_sin(p1, x, prec) 1517 R = R.compose(p_cos.ring).compose(p_sin.ring) 1518 p_cos.set_ring(R) 1519 p_sin.set_ring(R) 1520 t1, t2 = R(sin(c_expr)), R(cos(c_expr)) 1521 return p_cos*t2 - p_sin*t1 1522 1523 # Series is calculated in terms of tan as its evaluation is fast. 1524 if len(p) > 20 and R.ngens == 1: 1525 t = rs_tan(p/2, x, prec) 1526 t2 = rs_square(t, x, prec) 1527 p1 = rs_series_inversion(1+t2, x, prec) 1528 return rs_mul(p1, 1 - t2, x, prec) 1529 one = R(1) 1530 n = 1 1531 c = [] 1532 for k in range(2, prec + 2, 2): 1533 c.append(one/n) 1534 c.append(0) 1535 n *= -k*(k - 1) 1536 return rs_series_from_list(p, c, x, prec) 1537 1538def rs_cos_sin(p, x, prec): 1539 r""" 1540 Return the tuple ``(rs_cos(p, x, prec)`, `rs_sin(p, x, prec))``. 1541 1542 Is faster than calling rs_cos and rs_sin separately 1543 """ 1544 if rs_is_puiseux(p, x): 1545 return rs_puiseux(rs_cos_sin, p, x, prec) 1546 t = rs_tan(p/2, x, prec) 1547 t2 = rs_square(t, x, prec) 1548 p1 = rs_series_inversion(1 + t2, x, prec) 1549 return (rs_mul(p1, 1 - t2, x, prec), rs_mul(p1, 2*t, x, prec)) 1550 1551def _atanh(p, x, prec): 1552 """ 1553 Expansion using formula 1554 1555 Faster for very small and univariate series 1556 """ 1557 R = p.ring 1558 one = R(1) 1559 c = [one] 1560 p2 = rs_square(p, x, prec) 1561 for k in range(1, prec): 1562 c.append(one/(2*k + 1)) 1563 s = rs_series_from_list(p2, c, x, prec) 1564 s = rs_mul(s, p, x, prec) 1565 return s 1566 1567def rs_atanh(p, x, prec): 1568 """ 1569 Hyperbolic arctangent of a series 1570 1571 Return the series expansion of the atanh of ``p``, about 0. 1572 1573 Examples 1574 ======== 1575 1576 >>> from sympy.polys.domains import QQ 1577 >>> from sympy.polys.rings import ring 1578 >>> from sympy.polys.ring_series import rs_atanh 1579 >>> R, x, y = ring('x, y', QQ) 1580 >>> rs_atanh(x + x*y, x, 4) 1581 1/3*x**3*y**3 + x**3*y**2 + x**3*y + 1/3*x**3 + x*y + x 1582 1583 See Also 1584 ======== 1585 1586 atanh 1587 """ 1588 if rs_is_puiseux(p, x): 1589 return rs_puiseux(rs_atanh, p, x, prec) 1590 R = p.ring 1591 const = 0 1592 if _has_constant_term(p, x): 1593 zm = R.zero_monom 1594 c = p[zm] 1595 if R.domain is EX: 1596 c_expr = c.as_expr() 1597 const = atanh(c_expr) 1598 elif isinstance(c, PolyElement): 1599 try: 1600 c_expr = c.as_expr() 1601 const = R(atanh(c_expr)) 1602 except ValueError: 1603 raise DomainError("The given series can't be expanded in " 1604 "this domain.") 1605 else: 1606 try: 1607 const = R(atanh(c)) 1608 except ValueError: 1609 raise DomainError("The given series can't be expanded in " 1610 "this domain.") 1611 1612 # Instead of using a closed form formula, we differentiate atanh(p) to get 1613 # `1/(1-p**2) * dp`, whose series expansion is much easier to calculate. 1614 # Finally we integrate to get back atanh 1615 dp = rs_diff(p, x) 1616 p1 = - rs_square(p, x, prec) + 1 1617 p1 = rs_series_inversion(p1, x, prec - 1) 1618 p1 = rs_mul(dp, p1, x, prec - 1) 1619 return rs_integrate(p1, x) + const 1620 1621def rs_sinh(p, x, prec): 1622 """ 1623 Hyperbolic sine of a series 1624 1625 Return the series expansion of the sinh of ``p``, about 0. 1626 1627 Examples 1628 ======== 1629 1630 >>> from sympy.polys.domains import QQ 1631 >>> from sympy.polys.rings import ring 1632 >>> from sympy.polys.ring_series import rs_sinh 1633 >>> R, x, y = ring('x, y', QQ) 1634 >>> rs_sinh(x + x*y, x, 4) 1635 1/6*x**3*y**3 + 1/2*x**3*y**2 + 1/2*x**3*y + 1/6*x**3 + x*y + x 1636 1637 See Also 1638 ======== 1639 1640 sinh 1641 """ 1642 if rs_is_puiseux(p, x): 1643 return rs_puiseux(rs_sinh, p, x, prec) 1644 t = rs_exp(p, x, prec) 1645 t1 = rs_series_inversion(t, x, prec) 1646 return (t - t1)/2 1647 1648def rs_cosh(p, x, prec): 1649 """ 1650 Hyperbolic cosine of a series 1651 1652 Return the series expansion of the cosh of ``p``, about 0. 1653 1654 Examples 1655 ======== 1656 1657 >>> from sympy.polys.domains import QQ 1658 >>> from sympy.polys.rings import ring 1659 >>> from sympy.polys.ring_series import rs_cosh 1660 >>> R, x, y = ring('x, y', QQ) 1661 >>> rs_cosh(x + x*y, x, 4) 1662 1/2*x**2*y**2 + x**2*y + 1/2*x**2 + 1 1663 1664 See Also 1665 ======== 1666 1667 cosh 1668 """ 1669 if rs_is_puiseux(p, x): 1670 return rs_puiseux(rs_cosh, p, x, prec) 1671 t = rs_exp(p, x, prec) 1672 t1 = rs_series_inversion(t, x, prec) 1673 return (t + t1)/2 1674 1675def _tanh(p, x, prec): 1676 r""" 1677 Helper function of :func:`rs_tanh` 1678 1679 Return the series expansion of tanh of a univariate series using Newton's 1680 method. It takes advantage of the fact that series expansion of atanh is 1681 easier than that of tanh. 1682 1683 See Also 1684 ======== 1685 1686 _tanh 1687 """ 1688 R = p.ring 1689 p1 = R(0) 1690 for precx in _giant_steps(prec): 1691 tmp = p - rs_atanh(p1, x, precx) 1692 tmp = rs_mul(tmp, 1 - rs_square(p1, x, prec), x, precx) 1693 p1 += tmp 1694 return p1 1695 1696def rs_tanh(p, x, prec): 1697 """ 1698 Hyperbolic tangent of a series 1699 1700 Return the series expansion of the tanh of ``p``, about 0. 1701 1702 Examples 1703 ======== 1704 1705 >>> from sympy.polys.domains import QQ 1706 >>> from sympy.polys.rings import ring 1707 >>> from sympy.polys.ring_series import rs_tanh 1708 >>> R, x, y = ring('x, y', QQ) 1709 >>> rs_tanh(x + x*y, x, 4) 1710 -1/3*x**3*y**3 - x**3*y**2 - x**3*y - 1/3*x**3 + x*y + x 1711 1712 See Also 1713 ======== 1714 1715 tanh 1716 """ 1717 if rs_is_puiseux(p, x): 1718 return rs_puiseux(rs_tanh, p, x, prec) 1719 R = p.ring 1720 const = 0 1721 if _has_constant_term(p, x): 1722 zm = R.zero_monom 1723 c = p[zm] 1724 if R.domain is EX: 1725 c_expr = c.as_expr() 1726 const = tanh(c_expr) 1727 elif isinstance(c, PolyElement): 1728 try: 1729 c_expr = c.as_expr() 1730 const = R(tanh(c_expr)) 1731 except ValueError: 1732 raise DomainError("The given series can't be expanded in " 1733 "this domain.") 1734 else: 1735 try: 1736 const = R(tanh(c)) 1737 except ValueError: 1738 raise DomainError("The given series can't be expanded in " 1739 "this domain.") 1740 p1 = p - c 1741 t1 = rs_tanh(p1, x, prec) 1742 t = rs_series_inversion(1 + const*t1, x, prec) 1743 return rs_mul(const + t1, t, x, prec) 1744 1745 if R.ngens == 1: 1746 return _tanh(p, x, prec) 1747 else: 1748 return rs_fun(p, _tanh, x, prec) 1749 1750def rs_newton(p, x, prec): 1751 """ 1752 Compute the truncated Newton sum of the polynomial ``p`` 1753 1754 Examples 1755 ======== 1756 1757 >>> from sympy.polys.domains import QQ 1758 >>> from sympy.polys.rings import ring 1759 >>> from sympy.polys.ring_series import rs_newton 1760 >>> R, x = ring('x', QQ) 1761 >>> p = x**2 - 2 1762 >>> rs_newton(p, x, 5) 1763 8*x**4 + 4*x**2 + 2 1764 """ 1765 deg = p.degree() 1766 p1 = _invert_monoms(p) 1767 p2 = rs_series_inversion(p1, x, prec) 1768 p3 = rs_mul(p1.diff(x), p2, x, prec) 1769 res = deg - p3*x 1770 return res 1771 1772def rs_hadamard_exp(p1, inverse=False): 1773 """ 1774 Return ``sum f_i/i!*x**i`` from ``sum f_i*x**i``, 1775 where ``x`` is the first variable. 1776 1777 If ``invers=True`` return ``sum f_i*i!*x**i`` 1778 1779 Examples 1780 ======== 1781 1782 >>> from sympy.polys.domains import QQ 1783 >>> from sympy.polys.rings import ring 1784 >>> from sympy.polys.ring_series import rs_hadamard_exp 1785 >>> R, x = ring('x', QQ) 1786 >>> p = 1 + x + x**2 + x**3 1787 >>> rs_hadamard_exp(p) 1788 1/6*x**3 + 1/2*x**2 + x + 1 1789 """ 1790 R = p1.ring 1791 if R.domain != QQ: 1792 raise NotImplementedError 1793 p = R.zero 1794 if not inverse: 1795 for exp1, v1 in p1.items(): 1796 p[exp1] = v1/int(ifac(exp1[0])) 1797 else: 1798 for exp1, v1 in p1.items(): 1799 p[exp1] = v1*int(ifac(exp1[0])) 1800 return p 1801 1802def rs_compose_add(p1, p2): 1803 """ 1804 compute the composed sum ``prod(p2(x - beta) for beta root of p1)`` 1805 1806 Examples 1807 ======== 1808 1809 >>> from sympy.polys.domains import QQ 1810 >>> from sympy.polys.rings import ring 1811 >>> from sympy.polys.ring_series import rs_compose_add 1812 >>> R, x = ring('x', QQ) 1813 >>> f = x**2 - 2 1814 >>> g = x**2 - 3 1815 >>> rs_compose_add(f, g) 1816 x**4 - 10*x**2 + 1 1817 1818 References 1819 ========== 1820 1821 .. [1] A. Bostan, P. Flajolet, B. Salvy and E. Schost 1822 "Fast Computation with Two Algebraic Numbers", 1823 (2002) Research Report 4579, Institut 1824 National de Recherche en Informatique et en Automatique 1825 """ 1826 R = p1.ring 1827 x = R.gens[0] 1828 prec = p1.degree()*p2.degree() + 1 1829 np1 = rs_newton(p1, x, prec) 1830 np1e = rs_hadamard_exp(np1) 1831 np2 = rs_newton(p2, x, prec) 1832 np2e = rs_hadamard_exp(np2) 1833 np3e = rs_mul(np1e, np2e, x, prec) 1834 np3 = rs_hadamard_exp(np3e, True) 1835 np3a = (np3[(0,)] - np3)/x 1836 q = rs_integrate(np3a, x) 1837 q = rs_exp(q, x, prec) 1838 q = _invert_monoms(q) 1839 q = q.primitive()[1] 1840 dp = p1.degree()*p2.degree() - q.degree() 1841 # `dp` is the multiplicity of the zeroes of the resultant; 1842 # these zeroes are missed in this computation so they are put here. 1843 # if p1 and p2 are monic irreducible polynomials, 1844 # there are zeroes in the resultant 1845 # if and only if p1 = p2 ; in fact in that case p1 and p2 have a 1846 # root in common, so gcd(p1, p2) != 1; being p1 and p2 irreducible 1847 # this means p1 = p2 1848 if dp: 1849 q = q*x**dp 1850 return q 1851 1852 1853_convert_func = { 1854 'sin': 'rs_sin', 1855 'cos': 'rs_cos', 1856 'exp': 'rs_exp', 1857 'tan': 'rs_tan', 1858 'log': 'rs_log' 1859 } 1860 1861def rs_min_pow(expr, series_rs, a): 1862 """Find the minimum power of `a` in the series expansion of expr""" 1863 series = 0 1864 n = 2 1865 while series == 0: 1866 series = _rs_series(expr, series_rs, a, n) 1867 n *= 2 1868 R = series.ring 1869 a = R(a) 1870 i = R.gens.index(a) 1871 return min(series, key=lambda t: t[i])[i] 1872 1873 1874def _rs_series(expr, series_rs, a, prec): 1875 # TODO Use _parallel_dict_from_expr instead of sring as sring is 1876 # inefficient. For details, read the todo in sring. 1877 args = expr.args 1878 R = series_rs.ring 1879 1880 # expr does not contain any function to be expanded 1881 if not any(arg.has(Function) for arg in args) and not expr.is_Function: 1882 return series_rs 1883 1884 if not expr.has(a): 1885 return series_rs 1886 1887 elif expr.is_Function: 1888 arg = args[0] 1889 if len(args) > 1: 1890 raise NotImplementedError 1891 R1, series = sring(arg, domain=QQ, expand=False, series=True) 1892 series_inner = _rs_series(arg, series, a, prec) 1893 1894 # Why do we need to compose these three rings? 1895 # 1896 # We want to use a simple domain (like ``QQ`` or ``RR``) but they don't 1897 # support symbolic coefficients. We need a ring that for example lets 1898 # us have `sin(1)` and `cos(1)` as coefficients if we are expanding 1899 # `sin(x + 1)`. The ``EX`` domain allows all symbolic coefficients, but 1900 # that makes it very complex and hence slow. 1901 # 1902 # To solve this problem, we add only those symbolic elements as 1903 # generators to our ring, that we need. Here, series_inner might 1904 # involve terms like `sin(4)`, `exp(a)`, etc, which are not there in 1905 # R1 or R. Hence, we compose these three rings to create one that has 1906 # the generators of all three. 1907 R = R.compose(R1).compose(series_inner.ring) 1908 series_inner = series_inner.set_ring(R) 1909 series = eval(_convert_func[str(expr.func)])(series_inner, 1910 R(a), prec) 1911 return series 1912 1913 elif expr.is_Mul: 1914 n = len(args) 1915 for arg in args: # XXX Looks redundant 1916 if not arg.is_Number: 1917 R1, _ = sring(arg, expand=False, series=True) 1918 R = R.compose(R1) 1919 min_pows = list(map(rs_min_pow, args, [R(arg) for arg in args], 1920 [a]*len(args))) 1921 sum_pows = sum(min_pows) 1922 series = R(1) 1923 1924 for i in range(n): 1925 _series = _rs_series(args[i], R(args[i]), a, prec - sum_pows + 1926 min_pows[i]) 1927 R = R.compose(_series.ring) 1928 _series = _series.set_ring(R) 1929 series = series.set_ring(R) 1930 series *= _series 1931 series = rs_trunc(series, R(a), prec) 1932 return series 1933 1934 elif expr.is_Add: 1935 n = len(args) 1936 series = R(0) 1937 for i in range(n): 1938 _series = _rs_series(args[i], R(args[i]), a, prec) 1939 R = R.compose(_series.ring) 1940 _series = _series.set_ring(R) 1941 series = series.set_ring(R) 1942 series += _series 1943 return series 1944 1945 elif expr.is_Pow: 1946 R1, _ = sring(expr.base, domain=QQ, expand=False, series=True) 1947 R = R.compose(R1) 1948 series_inner = _rs_series(expr.base, R(expr.base), a, prec) 1949 return rs_pow(series_inner, expr.exp, series_inner.ring(a), prec) 1950 1951 # The `is_constant` method is buggy hence we check it at the end. 1952 # See issue #9786 for details. 1953 elif isinstance(expr, Expr) and expr.is_constant(): 1954 return sring(expr, domain=QQ, expand=False, series=True)[1] 1955 1956 else: 1957 raise NotImplementedError 1958 1959def rs_series(expr, a, prec): 1960 """Return the series expansion of an expression about 0. 1961 1962 Parameters 1963 ========== 1964 1965 expr : :class:`Expr` 1966 a : :class:`Symbol` with respect to which expr is to be expanded 1967 prec : order of the series expansion 1968 1969 Currently supports multivariate Taylor series expansion. This is much 1970 faster that Sympy's series method as it uses sparse polynomial operations. 1971 1972 It automatically creates the simplest ring required to represent the series 1973 expansion through repeated calls to sring. 1974 1975 Examples 1976 ======== 1977 1978 >>> from sympy.polys.ring_series import rs_series 1979 >>> from sympy.functions import sin, cos, exp, tan 1980 >>> from sympy.core import symbols 1981 >>> from sympy.polys.domains import QQ 1982 >>> a, b, c = symbols('a, b, c') 1983 >>> rs_series(sin(a) + exp(a), a, 5) 1984 1/24*a**4 + 1/2*a**2 + 2*a + 1 1985 >>> series = rs_series(tan(a + b)*cos(a + c), a, 2) 1986 >>> series.as_expr() 1987 -a*sin(c)*tan(b) + a*cos(c)*tan(b)**2 + a*cos(c) + cos(c)*tan(b) 1988 >>> series = rs_series(exp(a**QQ(1,3) + a**QQ(2, 5)), a, 1) 1989 >>> series.as_expr() 1990 a**(11/15) + a**(4/5)/2 + a**(2/5) + a**(2/3)/2 + a**(1/3) + 1 1991 1992 """ 1993 R, series = sring(expr, domain=QQ, expand=False, series=True) 1994 if a not in R.symbols: 1995 R = R.add_gens([a, ]) 1996 series = series.set_ring(R) 1997 series = _rs_series(expr, series, a, prec) 1998 R = series.ring 1999 gen = R(a) 2000 prec_got = series.degree(gen) + 1 2001 2002 if prec_got >= prec: 2003 return rs_trunc(series, gen, prec) 2004 else: 2005 # increase the requested number of terms to get the desired 2006 # number keep increasing (up to 9) until the received order 2007 # is different than the original order and then predict how 2008 # many additional terms are needed 2009 for more in range(1, 9): 2010 p1 = _rs_series(expr, series, a, prec=prec + more) 2011 gen = gen.set_ring(p1.ring) 2012 new_prec = p1.degree(gen) + 1 2013 if new_prec != prec_got: 2014 prec_do = ceiling(prec + (prec - prec_got)*more/(new_prec - 2015 prec_got)) 2016 p1 = _rs_series(expr, series, a, prec=prec_do) 2017 while p1.degree(gen) + 1 < prec: 2018 p1 = _rs_series(expr, series, a, prec=prec_do) 2019 gen = gen.set_ring(p1.ring) 2020 prec_do *= 2 2021 break 2022 else: 2023 break 2024 else: 2025 raise ValueError('Could not calculate %s terms for %s' 2026 % (str(prec), expr)) 2027 return rs_trunc(p1, gen, prec) 2028