Coverage for /usr/lib/python3/dist-packages/sympy/polys/ring_series.py: 6%
967 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Power series evaluation and manipulation using sparse Polynomials
3Implementing a new function
4---------------------------
6There are a few things to be kept in mind when adding a new function here::
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:
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.
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 does not 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).
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.
28Extending rs_series
29-------------------
31To make a function work with rs_series you need to do two things::
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.
40Look at rs_sin and rs_series for further reference.
42"""
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.functions import sin, cos, tan, atan, exp, atanh, tanh, log, ceiling
53from sympy.utilities.misc import as_int
54from mpmath.libmp.libintmath import giant_steps
55import math
58def _invert_monoms(p1):
59 """
60 Compute ``x**n * p1(1/x)`` for a univariate polynomial ``p1`` in ``x``.
62 Examples
63 ========
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
73 See Also
74 ========
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 mvi, cvi in zip(mv, cv):
86 p[(deg - mvi[0],)] = cvi
87 return p
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
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)``
101 Examples
102 ========
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
123def rs_is_puiseux(p, x):
124 """
125 Test if ``p`` is Puiseux series in ``x``.
127 Raise an exception if it has a negative power in ``x``.
129 Examples
130 ========
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
148def rs_puiseux(f, p, x, prec):
149 """
150 Return the puiseux series for `f(p, x, prec)`.
152 To be used when function ``f`` is implemented only for regular series.
154 Examples
155 ========
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
187def rs_puiseux2(f, p, q, x, prec):
188 """
189 Return the puiseux series for `f(p, q, x, prec)`.
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
212def rs_mul(p1, p2, x, prec):
213 """
214 Return the product of the given two series, modulo ``O(x**prec)``.
216 ``x`` is the series variable or its position in the generators.
218 Examples
219 ========
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('p2 must be a polynomial')
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
260 p.strip_zero()
261 return p
263def rs_square(p1, x, prec):
264 """
265 Square the series modulo ``O(x**prec)``
267 Examples
268 ========
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
303def rs_pow(p1, n, x, prec):
304 """
305 Return ``p1**n`` modulo ``O(x**prec)``
307 Examples
308 ========
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
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
357def rs_subs(p, rules, x, prec):
358 """
359 Substitution with truncation according to the mapping in ``rules``.
361 Return a series with precision ``prec`` in the generator ``x``
363 Note that substitutions are not done one after the other
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
375 which differs from
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
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.
387 Examples
388 ========
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
425def _has_constant_term(p, x):
426 """
427 Check if ``p`` has a constant term in ``x``
429 Examples
430 ========
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
451def _get_constant_term(p, x):
452 """Return constant term in p with respect to x
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
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
478def _series_inversion1(p, x, prec):
479 """
480 Univariate series inversion ``1/p`` modulo ``O(x**prec)``.
482 The Newton method is used.
484 Examples
485 ========
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]
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)
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
524def rs_series_inversion(p, x, prec):
525 """
526 Multivariate series inversion ``1/p`` modulo ``O(x**prec)``.
528 Examples
529 ========
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")
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
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
575def rs_series_reversion(p, x, n, y):
576 r"""
577 Reversion of a series.
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.
582 $f(x) = \sum_{k=2}^{n-1} a_kx_k$
584 Parameters
585 ==========
587 a_k : Can depend polynomially on other variables, not indicated.
588 x : Variable with name x.
589 y : Variable with name y.
591 Returns
592 =======
594 Solve $p = y$, that is, given $ax + f(x) - y = 0$,
595 find the solution $x = r(y)$ up to $O(y^n)$.
597 Algorithm
598 =========
600 If $r_i$ is the solution at order $i$, then:
601 $ar_i + f(r_i) - y = O\left(y^{i + 1}\right)$
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)$
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$
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$
614 Examples
615 ========
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
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)`.
651 It reduces the number of multiplications by summing concurrently.
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`
657 Examples
658 ========
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
674 """
676 # TODO: Add this when it is documented in Sphinx
677 """
678 See Also
679 ========
681 sympy.polys.rings.PolyRing.compose
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
736def rs_diff(p, x):
737 """
738 Return partial derivative of ``p`` with respect to ``x``.
740 Parameters
741 ==========
743 x : :class:`~.PolyElement` with respect to which ``p`` is differentiated.
745 Examples
746 ========
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
768def rs_integrate(p, x):
769 """
770 Integrate ``p`` with respect to ``x``.
772 Parameters
773 ==========
775 x : :class:`~.PolyElement` with respect to which ``p`` is integrated.
777 Examples
778 ========
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)
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
800def rs_fun(p, f, *args):
801 r"""
802 Function of a multivariate series computed by substitution.
804 The case with f method name is used to compute `rs\_tan` and `rs\_nth\_root`
805 of a multivariate series:
807 `rs\_fun(p, tan, iv, prec)`
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
813 Parameters
814 ==========
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.
821 Examples
822 ========
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
856def mul_xin(p, i, n):
857 r"""
858 Return `p*x_i**n`.
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
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
889def _nth_root1(p, n, x, prec):
890 """
891 Univariate series expansion of the nth root of ``p``.
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)
924def rs_nth_root(p, n, x, prec):
925 """
926 Multivariate series expansion of the nth root of ``p``.
928 Parameters
929 ==========
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.
939 Notes
940 =====
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 cannot yet handle
945 roots of symbols.
947 Examples
948 ========
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
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 cannot 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 cannot 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
1000def rs_log(p, x, prec):
1001 """
1002 The Logarithm of ``p`` modulo ``O(x**prec)``.
1004 Notes
1005 =====
1007 Truncation of ``integral dx p**-1*d p/dx`` is used.
1009 Examples
1010 ========
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 cannot be expanded in "
1049 "this domain.")
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
1057def rs_LambertW(p, x, prec):
1058 """
1059 Calculate the series expansion of the principal branch of the Lambert W
1060 function.
1062 Examples
1063 ========
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
1072 See Also
1073 ========
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
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
1106def rs_exp(p, x, prec):
1107 """
1108 Exponentiation of a series modulo ``O(x**prec)``
1110 Examples
1111 ========
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 cannot be expanded in "
1143 "this domain.")
1144 p1 = p - c
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)
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
1160 r = rs_series_from_list(p, c, x, prec)
1161 return r
1163def _atan(p, iv, prec):
1164 """
1165 Expansion using formula.
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
1179def rs_atan(p, x, prec):
1180 """
1181 The arctangent of a series
1183 Return the series expansion of the atan of ``p``, about 0.
1185 Examples
1186 ========
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
1195 See Also
1196 ========
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 cannot be expanded in "
1216 "this domain.")
1217 else:
1218 try:
1219 const = R(atan(c))
1220 except ValueError:
1221 raise DomainError("The given series cannot be expanded in "
1222 "this domain.")
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
1233def rs_asin(p, x, prec):
1234 """
1235 Arcsine of a series
1237 Return the series expansion of the asin of ``p``, about 0.
1239 Examples
1240 ========
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
1249 See Also
1250 ========
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)
1275 else:
1276 raise NotImplementedError
1278def _tan1(p, x, prec):
1279 r"""
1280 Helper function of :func:`rs_tan`.
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.
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
1299def rs_tan(p, x, prec):
1300 """
1301 Tangent of a series.
1303 Return the series expansion of the tan of ``p``, about 0.
1305 Examples
1306 ========
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
1315 See Also
1316 ========
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 cannot be expanded in "
1345 "this domain.")
1346 p1 = p - c
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)
1354 if R.ngens == 1:
1355 return _tan1(p, x, prec)
1356 else:
1357 return rs_fun(p, rs_tan, x, prec)
1359def rs_cot(p, x, prec):
1360 """
1361 Cotangent of a series
1363 Return the series expansion of the cot of ``p``, about 0.
1365 Examples
1366 ========
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)
1375 See Also
1376 ========
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
1395def rs_sin(p, x, prec):
1396 """
1397 Sine of a series
1399 Return the series expansion of the sin of ``p``, about 0.
1401 Examples
1402 ========
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)
1413 See Also
1414 ========
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 cannot be expanded in "
1443 "this domain.")
1444 p1 = p - c
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
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)
1465def rs_cos(p, x, prec):
1466 """
1467 Cosine of a series
1469 Return the series expansion of the cos of ``p``, about 0.
1471 Examples
1472 ========
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)
1483 See Also
1484 ========
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 cannot be expanded in "
1510 "this domain.")
1511 p1 = p - c
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
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)
1538def rs_cos_sin(p, x, prec):
1539 r"""
1540 Return the tuple ``(rs_cos(p, x, prec)`, `rs_sin(p, x, prec))``.
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))
1551def _atanh(p, x, prec):
1552 """
1553 Expansion using formula
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
1567def rs_atanh(p, x, prec):
1568 """
1569 Hyperbolic arctangent of a series
1571 Return the series expansion of the atanh of ``p``, about 0.
1573 Examples
1574 ========
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
1583 See Also
1584 ========
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 cannot be expanded in "
1604 "this domain.")
1605 else:
1606 try:
1607 const = R(atanh(c))
1608 except ValueError:
1609 raise DomainError("The given series cannot be expanded in "
1610 "this domain.")
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
1621def rs_sinh(p, x, prec):
1622 """
1623 Hyperbolic sine of a series
1625 Return the series expansion of the sinh of ``p``, about 0.
1627 Examples
1628 ========
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
1637 See Also
1638 ========
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
1648def rs_cosh(p, x, prec):
1649 """
1650 Hyperbolic cosine of a series
1652 Return the series expansion of the cosh of ``p``, about 0.
1654 Examples
1655 ========
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
1664 See Also
1665 ========
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
1675def _tanh(p, x, prec):
1676 r"""
1677 Helper function of :func:`rs_tanh`
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.
1683 See Also
1684 ========
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
1696def rs_tanh(p, x, prec):
1697 """
1698 Hyperbolic tangent of a series
1700 Return the series expansion of the tanh of ``p``, about 0.
1702 Examples
1703 ========
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
1712 See Also
1713 ========
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 cannot be expanded in "
1733 "this domain.")
1734 else:
1735 try:
1736 const = R(tanh(c))
1737 except ValueError:
1738 raise DomainError("The given series cannot 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)
1745 if R.ngens == 1:
1746 return _tanh(p, x, prec)
1747 else:
1748 return rs_fun(p, _tanh, x, prec)
1750def rs_newton(p, x, prec):
1751 """
1752 Compute the truncated Newton sum of the polynomial ``p``
1754 Examples
1755 ========
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
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.
1777 If ``invers=True`` return ``sum f_i*i!*x**i``
1779 Examples
1780 ========
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
1802def rs_compose_add(p1, p2):
1803 """
1804 compute the composed sum ``prod(p2(x - beta) for beta root of p1)``
1806 Examples
1807 ========
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
1818 References
1819 ==========
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
1853_convert_func = {
1854 'sin': 'rs_sin',
1855 'cos': 'rs_cos',
1856 'exp': 'rs_exp',
1857 'tan': 'rs_tan',
1858 'log': 'rs_log'
1859 }
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]
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
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
1884 if not expr.has(a):
1885 return series_rs
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)
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
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)
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
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
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)
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]
1956 else:
1957 raise NotImplementedError
1959def rs_series(expr, a, prec):
1960 """Return the series expansion of an expression about 0.
1962 Parameters
1963 ==========
1965 expr : :class:`Expr`
1966 a : :class:`Symbol` with respect to which expr is to be expanded
1967 prec : order of the series expansion
1969 Currently supports multivariate Taylor series expansion. This is much
1970 faster that SymPy's series method as it uses sparse polynomial operations.
1972 It automatically creates the simplest ring required to represent the series
1973 expansion through repeated calls to sring.
1975 Examples
1976 ========
1978 >>> from sympy.polys.ring_series import rs_series
1979 >>> from sympy import sin, cos, exp, tan, symbols, QQ
1980 >>> a, b, c = symbols('a, b, c')
1981 >>> rs_series(sin(a) + exp(a), a, 5)
1982 1/24*a**4 + 1/2*a**2 + 2*a + 1
1983 >>> series = rs_series(tan(a + b)*cos(a + c), a, 2)
1984 >>> series.as_expr()
1985 -a*sin(c)*tan(b) + a*cos(c)*tan(b)**2 + a*cos(c) + cos(c)*tan(b)
1986 >>> series = rs_series(exp(a**QQ(1,3) + a**QQ(2, 5)), a, 1)
1987 >>> series.as_expr()
1988 a**(11/15) + a**(4/5)/2 + a**(2/5) + a**(2/3)/2 + a**(1/3) + 1
1990 """
1991 R, series = sring(expr, domain=QQ, expand=False, series=True)
1992 if a not in R.symbols:
1993 R = R.add_gens([a, ])
1994 series = series.set_ring(R)
1995 series = _rs_series(expr, series, a, prec)
1996 R = series.ring
1997 gen = R(a)
1998 prec_got = series.degree(gen) + 1
2000 if prec_got >= prec:
2001 return rs_trunc(series, gen, prec)
2002 else:
2003 # increase the requested number of terms to get the desired
2004 # number keep increasing (up to 9) until the received order
2005 # is different than the original order and then predict how
2006 # many additional terms are needed
2007 for more in range(1, 9):
2008 p1 = _rs_series(expr, series, a, prec=prec + more)
2009 gen = gen.set_ring(p1.ring)
2010 new_prec = p1.degree(gen) + 1
2011 if new_prec != prec_got:
2012 prec_do = ceiling(prec + (prec - prec_got)*more/(new_prec -
2013 prec_got))
2014 p1 = _rs_series(expr, series, a, prec=prec_do)
2015 while p1.degree(gen) + 1 < prec:
2016 p1 = _rs_series(expr, series, a, prec=prec_do)
2017 gen = gen.set_ring(p1.ring)
2018 prec_do *= 2
2019 break
2020 else:
2021 break
2022 else:
2023 raise ValueError('Could not calculate %s terms for %s'
2024 % (str(prec), expr))
2025 return rs_trunc(p1, gen, prec)