1"""
2=================================================
3Power Series (:mod:`numpy.polynomial.polynomial`)
4=================================================
5
6This module provides a number of objects (mostly functions) useful for
7dealing with polynomials, including a `Polynomial` class that
8encapsulates the usual arithmetic operations. (General information
9on how this module represents and works with polynomial objects is in
10the docstring for its "parent" sub-package, `numpy.polynomial`).
11
12Classes
13-------
14.. autosummary::
15 :toctree: generated/
16
17 Polynomial
18
19Constants
20---------
21.. autosummary::
22 :toctree: generated/
23
24 polydomain
25 polyzero
26 polyone
27 polyx
28
29Arithmetic
30----------
31.. autosummary::
32 :toctree: generated/
33
34 polyadd
35 polysub
36 polymulx
37 polymul
38 polydiv
39 polypow
40 polyval
41 polyval2d
42 polyval3d
43 polygrid2d
44 polygrid3d
45
46Calculus
47--------
48.. autosummary::
49 :toctree: generated/
50
51 polyder
52 polyint
53
54Misc Functions
55--------------
56.. autosummary::
57 :toctree: generated/
58
59 polyfromroots
60 polyroots
61 polyvalfromroots
62 polyvander
63 polyvander2d
64 polyvander3d
65 polycompanion
66 polyfit
67 polytrim
68 polyline
69
70See Also
71--------
72`numpy.polynomial`
73
74"""
75__all__ = [
76 'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd',
77 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval',
78 'polyvalfromroots', 'polyder', 'polyint', 'polyfromroots', 'polyvander',
79 'polyfit', 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d',
80 'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d']
81
82import numpy as np
83import numpy.linalg as la
84from numpy.core.multiarray import normalize_axis_index
85
86from . import polyutils as pu
87from ._polybase import ABCPolyBase
88
89polytrim = pu.trimcoef
90
91#
92# These are constant arrays are of integer type so as to be compatible
93# with the widest range of other types, such as Decimal.
94#
95
96# Polynomial default domain.
97polydomain = np.array([-1, 1])
98
99# Polynomial coefficients representing zero.
100polyzero = np.array([0])
101
102# Polynomial coefficients representing one.
103polyone = np.array([1])
104
105# Polynomial coefficients representing the identity x.
106polyx = np.array([0, 1])
107
108#
109# Polynomial series functions
110#
111
112
113def polyline(off, scl):
114 """
115 Returns an array representing a linear polynomial.
116
117 Parameters
118 ----------
119 off, scl : scalars
120 The "y-intercept" and "slope" of the line, respectively.
121
122 Returns
123 -------
124 y : ndarray
125 This module's representation of the linear polynomial ``off +
126 scl*x``.
127
128 See Also
129 --------
130 numpy.polynomial.chebyshev.chebline
131 numpy.polynomial.legendre.legline
132 numpy.polynomial.laguerre.lagline
133 numpy.polynomial.hermite.hermline
134 numpy.polynomial.hermite_e.hermeline
135
136 Examples
137 --------
138 >>> from numpy.polynomial import polynomial as P
139 >>> P.polyline(1,-1)
140 array([ 1, -1])
141 >>> P.polyval(1, P.polyline(1,-1)) # should be 0
142 0.0
143
144 """
145 if scl != 0:
146 return np.array([off, scl])
147 else:
148 return np.array([off])
149
150
151def polyfromroots(roots):
152 """
153 Generate a monic polynomial with given roots.
154
155 Return the coefficients of the polynomial
156
157 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
158
159 where the ``r_n`` are the roots specified in `roots`. If a zero has
160 multiplicity n, then it must appear in `roots` n times. For instance,
161 if 2 is a root of multiplicity three and 3 is a root of multiplicity 2,
162 then `roots` looks something like [2, 2, 2, 3, 3]. The roots can appear
163 in any order.
164
165 If the returned coefficients are `c`, then
166
167 .. math:: p(x) = c_0 + c_1 * x + ... + x^n
168
169 The coefficient of the last term is 1 for monic polynomials in this
170 form.
171
172 Parameters
173 ----------
174 roots : array_like
175 Sequence containing the roots.
176
177 Returns
178 -------
179 out : ndarray
180 1-D array of the polynomial's coefficients If all the roots are
181 real, then `out` is also real, otherwise it is complex. (see
182 Examples below).
183
184 See Also
185 --------
186 numpy.polynomial.chebyshev.chebfromroots
187 numpy.polynomial.legendre.legfromroots
188 numpy.polynomial.laguerre.lagfromroots
189 numpy.polynomial.hermite.hermfromroots
190 numpy.polynomial.hermite_e.hermefromroots
191
192 Notes
193 -----
194 The coefficients are determined by multiplying together linear factors
195 of the form ``(x - r_i)``, i.e.
196
197 .. math:: p(x) = (x - r_0) (x - r_1) ... (x - r_n)
198
199 where ``n == len(roots) - 1``; note that this implies that ``1`` is always
200 returned for :math:`a_n`.
201
202 Examples
203 --------
204 >>> from numpy.polynomial import polynomial as P
205 >>> P.polyfromroots((-1,0,1)) # x(x - 1)(x + 1) = x^3 - x
206 array([ 0., -1., 0., 1.])
207 >>> j = complex(0,1)
208 >>> P.polyfromroots((-j,j)) # complex returned, though values are real
209 array([1.+0.j, 0.+0.j, 1.+0.j])
210
211 """
212 return pu._fromroots(polyline, polymul, roots)
213
214
215def polyadd(c1, c2):
216 """
217 Add one polynomial to another.
218
219 Returns the sum of two polynomials `c1` + `c2`. The arguments are
220 sequences of coefficients from lowest order term to highest, i.e.,
221 [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``.
222
223 Parameters
224 ----------
225 c1, c2 : array_like
226 1-D arrays of polynomial coefficients ordered from low to high.
227
228 Returns
229 -------
230 out : ndarray
231 The coefficient array representing their sum.
232
233 See Also
234 --------
235 polysub, polymulx, polymul, polydiv, polypow
236
237 Examples
238 --------
239 >>> from numpy.polynomial import polynomial as P
240 >>> c1 = (1,2,3)
241 >>> c2 = (3,2,1)
242 >>> sum = P.polyadd(c1,c2); sum
243 array([4., 4., 4.])
244 >>> P.polyval(2, sum) # 4 + 4(2) + 4(2**2)
245 28.0
246
247 """
248 return pu._add(c1, c2)
249
250
251def polysub(c1, c2):
252 """
253 Subtract one polynomial from another.
254
255 Returns the difference of two polynomials `c1` - `c2`. The arguments
256 are sequences of coefficients from lowest order term to highest, i.e.,
257 [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``.
258
259 Parameters
260 ----------
261 c1, c2 : array_like
262 1-D arrays of polynomial coefficients ordered from low to
263 high.
264
265 Returns
266 -------
267 out : ndarray
268 Of coefficients representing their difference.
269
270 See Also
271 --------
272 polyadd, polymulx, polymul, polydiv, polypow
273
274 Examples
275 --------
276 >>> from numpy.polynomial import polynomial as P
277 >>> c1 = (1,2,3)
278 >>> c2 = (3,2,1)
279 >>> P.polysub(c1,c2)
280 array([-2., 0., 2.])
281 >>> P.polysub(c2,c1) # -P.polysub(c1,c2)
282 array([ 2., 0., -2.])
283
284 """
285 return pu._sub(c1, c2)
286
287
288def polymulx(c):
289 """Multiply a polynomial by x.
290
291 Multiply the polynomial `c` by x, where x is the independent
292 variable.
293
294
295 Parameters
296 ----------
297 c : array_like
298 1-D array of polynomial coefficients ordered from low to
299 high.
300
301 Returns
302 -------
303 out : ndarray
304 Array representing the result of the multiplication.
305
306 See Also
307 --------
308 polyadd, polysub, polymul, polydiv, polypow
309
310 Notes
311 -----
312
313 .. versionadded:: 1.5.0
314
315 """
316 # c is a trimmed copy
317 [c] = pu.as_series([c])
318 # The zero series needs special treatment
319 if len(c) == 1 and c[0] == 0:
320 return c
321
322 prd = np.empty(len(c) + 1, dtype=c.dtype)
323 prd[0] = c[0]*0
324 prd[1:] = c
325 return prd
326
327
328def polymul(c1, c2):
329 """
330 Multiply one polynomial by another.
331
332 Returns the product of two polynomials `c1` * `c2`. The arguments are
333 sequences of coefficients, from lowest order term to highest, e.g.,
334 [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2.``
335
336 Parameters
337 ----------
338 c1, c2 : array_like
339 1-D arrays of coefficients representing a polynomial, relative to the
340 "standard" basis, and ordered from lowest order term to highest.
341
342 Returns
343 -------
344 out : ndarray
345 Of the coefficients of their product.
346
347 See Also
348 --------
349 polyadd, polysub, polymulx, polydiv, polypow
350
351 Examples
352 --------
353 >>> from numpy.polynomial import polynomial as P
354 >>> c1 = (1,2,3)
355 >>> c2 = (3,2,1)
356 >>> P.polymul(c1,c2)
357 array([ 3., 8., 14., 8., 3.])
358
359 """
360 # c1, c2 are trimmed copies
361 [c1, c2] = pu.as_series([c1, c2])
362 ret = np.convolve(c1, c2)
363 return pu.trimseq(ret)
364
365
366def polydiv(c1, c2):
367 """
368 Divide one polynomial by another.
369
370 Returns the quotient-with-remainder of two polynomials `c1` / `c2`.
371 The arguments are sequences of coefficients, from lowest order term
372 to highest, e.g., [1,2,3] represents ``1 + 2*x + 3*x**2``.
373
374 Parameters
375 ----------
376 c1, c2 : array_like
377 1-D arrays of polynomial coefficients ordered from low to high.
378
379 Returns
380 -------
381 [quo, rem] : ndarrays
382 Of coefficient series representing the quotient and remainder.
383
384 See Also
385 --------
386 polyadd, polysub, polymulx, polymul, polypow
387
388 Examples
389 --------
390 >>> from numpy.polynomial import polynomial as P
391 >>> c1 = (1,2,3)
392 >>> c2 = (3,2,1)
393 >>> P.polydiv(c1,c2)
394 (array([3.]), array([-8., -4.]))
395 >>> P.polydiv(c2,c1)
396 (array([ 0.33333333]), array([ 2.66666667, 1.33333333])) # may vary
397
398 """
399 # c1, c2 are trimmed copies
400 [c1, c2] = pu.as_series([c1, c2])
401 if c2[-1] == 0:
402 raise ZeroDivisionError()
403
404 # note: this is more efficient than `pu._div(polymul, c1, c2)`
405 lc1 = len(c1)
406 lc2 = len(c2)
407 if lc1 < lc2:
408 return c1[:1]*0, c1
409 elif lc2 == 1:
410 return c1/c2[-1], c1[:1]*0
411 else:
412 dlen = lc1 - lc2
413 scl = c2[-1]
414 c2 = c2[:-1]/scl
415 i = dlen
416 j = lc1 - 1
417 while i >= 0:
418 c1[i:j] -= c2*c1[j]
419 i -= 1
420 j -= 1
421 return c1[j+1:]/scl, pu.trimseq(c1[:j+1])
422
423
424def polypow(c, pow, maxpower=None):
425 """Raise a polynomial to a power.
426
427 Returns the polynomial `c` raised to the power `pow`. The argument
428 `c` is a sequence of coefficients ordered from low to high. i.e.,
429 [1,2,3] is the series ``1 + 2*x + 3*x**2.``
430
431 Parameters
432 ----------
433 c : array_like
434 1-D array of array of series coefficients ordered from low to
435 high degree.
436 pow : integer
437 Power to which the series will be raised
438 maxpower : integer, optional
439 Maximum power allowed. This is mainly to limit growth of the series
440 to unmanageable size. Default is 16
441
442 Returns
443 -------
444 coef : ndarray
445 Power series of power.
446
447 See Also
448 --------
449 polyadd, polysub, polymulx, polymul, polydiv
450
451 Examples
452 --------
453 >>> from numpy.polynomial import polynomial as P
454 >>> P.polypow([1,2,3], 2)
455 array([ 1., 4., 10., 12., 9.])
456
457 """
458 # note: this is more efficient than `pu._pow(polymul, c1, c2)`, as it
459 # avoids calling `as_series` repeatedly
460 return pu._pow(np.convolve, c, pow, maxpower)
461
462
463def polyder(c, m=1, scl=1, axis=0):
464 """
465 Differentiate a polynomial.
466
467 Returns the polynomial coefficients `c` differentiated `m` times along
468 `axis`. At each iteration the result is multiplied by `scl` (the
469 scaling factor is for use in a linear change of variable). The
470 argument `c` is an array of coefficients from low to high degree along
471 each axis, e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``
472 while [[1,2],[1,2]] represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is
473 ``x`` and axis=1 is ``y``.
474
475 Parameters
476 ----------
477 c : array_like
478 Array of polynomial coefficients. If c is multidimensional the
479 different axis correspond to different variables with the degree
480 in each axis given by the corresponding index.
481 m : int, optional
482 Number of derivatives taken, must be non-negative. (Default: 1)
483 scl : scalar, optional
484 Each differentiation is multiplied by `scl`. The end result is
485 multiplication by ``scl**m``. This is for use in a linear change
486 of variable. (Default: 1)
487 axis : int, optional
488 Axis over which the derivative is taken. (Default: 0).
489
490 .. versionadded:: 1.7.0
491
492 Returns
493 -------
494 der : ndarray
495 Polynomial coefficients of the derivative.
496
497 See Also
498 --------
499 polyint
500
501 Examples
502 --------
503 >>> from numpy.polynomial import polynomial as P
504 >>> c = (1,2,3,4) # 1 + 2x + 3x**2 + 4x**3
505 >>> P.polyder(c) # (d/dx)(c) = 2 + 6x + 12x**2
506 array([ 2., 6., 12.])
507 >>> P.polyder(c,3) # (d**3/dx**3)(c) = 24
508 array([24.])
509 >>> P.polyder(c,scl=-1) # (d/d(-x))(c) = -2 - 6x - 12x**2
510 array([ -2., -6., -12.])
511 >>> P.polyder(c,2,-1) # (d**2/d(-x)**2)(c) = 6 + 24x
512 array([ 6., 24.])
513
514 """
515 c = np.array(c, ndmin=1, copy=True)
516 if c.dtype.char in '?bBhHiIlLqQpP':
517 # astype fails with NA
518 c = c + 0.0
519 cdt = c.dtype
520 cnt = pu._deprecate_as_int(m, "the order of derivation")
521 iaxis = pu._deprecate_as_int(axis, "the axis")
522 if cnt < 0:
523 raise ValueError("The order of derivation must be non-negative")
524 iaxis = normalize_axis_index(iaxis, c.ndim)
525
526 if cnt == 0:
527 return c
528
529 c = np.moveaxis(c, iaxis, 0)
530 n = len(c)
531 if cnt >= n:
532 c = c[:1]*0
533 else:
534 for i in range(cnt):
535 n = n - 1
536 c *= scl
537 der = np.empty((n,) + c.shape[1:], dtype=cdt)
538 for j in range(n, 0, -1):
539 der[j - 1] = j*c[j]
540 c = der
541 c = np.moveaxis(c, 0, iaxis)
542 return c
543
544
545def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
546 """
547 Integrate a polynomial.
548
549 Returns the polynomial coefficients `c` integrated `m` times from
550 `lbnd` along `axis`. At each iteration the resulting series is
551 **multiplied** by `scl` and an integration constant, `k`, is added.
552 The scaling factor is for use in a linear change of variable. ("Buyer
553 beware": note that, depending on what one is doing, one may want `scl`
554 to be the reciprocal of what one might expect; for more information,
555 see the Notes section below.) The argument `c` is an array of
556 coefficients, from low to high degree along each axis, e.g., [1,2,3]
557 represents the polynomial ``1 + 2*x + 3*x**2`` while [[1,2],[1,2]]
558 represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is ``x`` and axis=1 is
559 ``y``.
560
561 Parameters
562 ----------
563 c : array_like
564 1-D array of polynomial coefficients, ordered from low to high.
565 m : int, optional
566 Order of integration, must be positive. (Default: 1)
567 k : {[], list, scalar}, optional
568 Integration constant(s). The value of the first integral at zero
569 is the first value in the list, the value of the second integral
570 at zero is the second value, etc. If ``k == []`` (the default),
571 all constants are set to zero. If ``m == 1``, a single scalar can
572 be given instead of a list.
573 lbnd : scalar, optional
574 The lower bound of the integral. (Default: 0)
575 scl : scalar, optional
576 Following each integration the result is *multiplied* by `scl`
577 before the integration constant is added. (Default: 1)
578 axis : int, optional
579 Axis over which the integral is taken. (Default: 0).
580
581 .. versionadded:: 1.7.0
582
583 Returns
584 -------
585 S : ndarray
586 Coefficient array of the integral.
587
588 Raises
589 ------
590 ValueError
591 If ``m < 1``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
592 ``np.ndim(scl) != 0``.
593
594 See Also
595 --------
596 polyder
597
598 Notes
599 -----
600 Note that the result of each integration is *multiplied* by `scl`. Why
601 is this important to note? Say one is making a linear change of
602 variable :math:`u = ax + b` in an integral relative to `x`. Then
603 :math:`dx = du/a`, so one will need to set `scl` equal to
604 :math:`1/a` - perhaps not what one would have first thought.
605
606 Examples
607 --------
608 >>> from numpy.polynomial import polynomial as P
609 >>> c = (1,2,3)
610 >>> P.polyint(c) # should return array([0, 1, 1, 1])
611 array([0., 1., 1., 1.])
612 >>> P.polyint(c,3) # should return array([0, 0, 0, 1/6, 1/12, 1/20])
613 array([ 0. , 0. , 0. , 0.16666667, 0.08333333, # may vary
614 0.05 ])
615 >>> P.polyint(c,k=3) # should return array([3, 1, 1, 1])
616 array([3., 1., 1., 1.])
617 >>> P.polyint(c,lbnd=-2) # should return array([6, 1, 1, 1])
618 array([6., 1., 1., 1.])
619 >>> P.polyint(c,scl=-2) # should return array([0, -2, -2, -2])
620 array([ 0., -2., -2., -2.])
621
622 """
623 c = np.array(c, ndmin=1, copy=True)
624 if c.dtype.char in '?bBhHiIlLqQpP':
625 # astype doesn't preserve mask attribute.
626 c = c + 0.0
627 cdt = c.dtype
628 if not np.iterable(k):
629 k = [k]
630 cnt = pu._deprecate_as_int(m, "the order of integration")
631 iaxis = pu._deprecate_as_int(axis, "the axis")
632 if cnt < 0:
633 raise ValueError("The order of integration must be non-negative")
634 if len(k) > cnt:
635 raise ValueError("Too many integration constants")
636 if np.ndim(lbnd) != 0:
637 raise ValueError("lbnd must be a scalar.")
638 if np.ndim(scl) != 0:
639 raise ValueError("scl must be a scalar.")
640 iaxis = normalize_axis_index(iaxis, c.ndim)
641
642 if cnt == 0:
643 return c
644
645 k = list(k) + [0]*(cnt - len(k))
646 c = np.moveaxis(c, iaxis, 0)
647 for i in range(cnt):
648 n = len(c)
649 c *= scl
650 if n == 1 and np.all(c[0] == 0):
651 c[0] += k[i]
652 else:
653 tmp = np.empty((n + 1,) + c.shape[1:], dtype=cdt)
654 tmp[0] = c[0]*0
655 tmp[1] = c[0]
656 for j in range(1, n):
657 tmp[j + 1] = c[j]/(j + 1)
658 tmp[0] += k[i] - polyval(lbnd, tmp)
659 c = tmp
660 c = np.moveaxis(c, 0, iaxis)
661 return c
662
663
664def polyval(x, c, tensor=True):
665 """
666 Evaluate a polynomial at points x.
667
668 If `c` is of length `n + 1`, this function returns the value
669
670 .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n
671
672 The parameter `x` is converted to an array only if it is a tuple or a
673 list, otherwise it is treated as a scalar. In either case, either `x`
674 or its elements must support multiplication and addition both with
675 themselves and with the elements of `c`.
676
677 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
678 `c` is multidimensional, then the shape of the result depends on the
679 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
680 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
681 scalars have shape (,).
682
683 Trailing zeros in the coefficients will be used in the evaluation, so
684 they should be avoided if efficiency is a concern.
685
686 Parameters
687 ----------
688 x : array_like, compatible object
689 If `x` is a list or tuple, it is converted to an ndarray, otherwise
690 it is left unchanged and treated as a scalar. In either case, `x`
691 or its elements must support addition and multiplication with
692 with themselves and with the elements of `c`.
693 c : array_like
694 Array of coefficients ordered so that the coefficients for terms of
695 degree n are contained in c[n]. If `c` is multidimensional the
696 remaining indices enumerate multiple polynomials. In the two
697 dimensional case the coefficients may be thought of as stored in
698 the columns of `c`.
699 tensor : boolean, optional
700 If True, the shape of the coefficient array is extended with ones
701 on the right, one for each dimension of `x`. Scalars have dimension 0
702 for this action. The result is that every column of coefficients in
703 `c` is evaluated for every element of `x`. If False, `x` is broadcast
704 over the columns of `c` for the evaluation. This keyword is useful
705 when `c` is multidimensional. The default value is True.
706
707 .. versionadded:: 1.7.0
708
709 Returns
710 -------
711 values : ndarray, compatible object
712 The shape of the returned array is described above.
713
714 See Also
715 --------
716 polyval2d, polygrid2d, polyval3d, polygrid3d
717
718 Notes
719 -----
720 The evaluation uses Horner's method.
721
722 Examples
723 --------
724 >>> from numpy.polynomial.polynomial import polyval
725 >>> polyval(1, [1,2,3])
726 6.0
727 >>> a = np.arange(4).reshape(2,2)
728 >>> a
729 array([[0, 1],
730 [2, 3]])
731 >>> polyval(a, [1,2,3])
732 array([[ 1., 6.],
733 [17., 34.]])
734 >>> coef = np.arange(4).reshape(2,2) # multidimensional coefficients
735 >>> coef
736 array([[0, 1],
737 [2, 3]])
738 >>> polyval([1,2], coef, tensor=True)
739 array([[2., 4.],
740 [4., 7.]])
741 >>> polyval([1,2], coef, tensor=False)
742 array([2., 7.])
743
744 """
745 c = np.array(c, ndmin=1, copy=False)
746 if c.dtype.char in '?bBhHiIlLqQpP':
747 # astype fails with NA
748 c = c + 0.0
749 if isinstance(x, (tuple, list)):
750 x = np.asarray(x)
751 if isinstance(x, np.ndarray) and tensor:
752 c = c.reshape(c.shape + (1,)*x.ndim)
753
754 c0 = c[-1] + x*0
755 for i in range(2, len(c) + 1):
756 c0 = c[-i] + c0*x
757 return c0
758
759
760def polyvalfromroots(x, r, tensor=True):
761 """
762 Evaluate a polynomial specified by its roots at points x.
763
764 If `r` is of length `N`, this function returns the value
765
766 .. math:: p(x) = \\prod_{n=1}^{N} (x - r_n)
767
768 The parameter `x` is converted to an array only if it is a tuple or a
769 list, otherwise it is treated as a scalar. In either case, either `x`
770 or its elements must support multiplication and addition both with
771 themselves and with the elements of `r`.
772
773 If `r` is a 1-D array, then `p(x)` will have the same shape as `x`. If `r`
774 is multidimensional, then the shape of the result depends on the value of
775 `tensor`. If `tensor` is ``True`` the shape will be r.shape[1:] + x.shape;
776 that is, each polynomial is evaluated at every value of `x`. If `tensor` is
777 ``False``, the shape will be r.shape[1:]; that is, each polynomial is
778 evaluated only for the corresponding broadcast value of `x`. Note that
779 scalars have shape (,).
780
781 .. versionadded:: 1.12
782
783 Parameters
784 ----------
785 x : array_like, compatible object
786 If `x` is a list or tuple, it is converted to an ndarray, otherwise
787 it is left unchanged and treated as a scalar. In either case, `x`
788 or its elements must support addition and multiplication with
789 with themselves and with the elements of `r`.
790 r : array_like
791 Array of roots. If `r` is multidimensional the first index is the
792 root index, while the remaining indices enumerate multiple
793 polynomials. For instance, in the two dimensional case the roots
794 of each polynomial may be thought of as stored in the columns of `r`.
795 tensor : boolean, optional
796 If True, the shape of the roots array is extended with ones on the
797 right, one for each dimension of `x`. Scalars have dimension 0 for this
798 action. The result is that every column of coefficients in `r` is
799 evaluated for every element of `x`. If False, `x` is broadcast over the
800 columns of `r` for the evaluation. This keyword is useful when `r` is
801 multidimensional. The default value is True.
802
803 Returns
804 -------
805 values : ndarray, compatible object
806 The shape of the returned array is described above.
807
808 See Also
809 --------
810 polyroots, polyfromroots, polyval
811
812 Examples
813 --------
814 >>> from numpy.polynomial.polynomial import polyvalfromroots
815 >>> polyvalfromroots(1, [1,2,3])
816 0.0
817 >>> a = np.arange(4).reshape(2,2)
818 >>> a
819 array([[0, 1],
820 [2, 3]])
821 >>> polyvalfromroots(a, [-1, 0, 1])
822 array([[-0., 0.],
823 [ 6., 24.]])
824 >>> r = np.arange(-2, 2).reshape(2,2) # multidimensional coefficients
825 >>> r # each column of r defines one polynomial
826 array([[-2, -1],
827 [ 0, 1]])
828 >>> b = [-2, 1]
829 >>> polyvalfromroots(b, r, tensor=True)
830 array([[-0., 3.],
831 [ 3., 0.]])
832 >>> polyvalfromroots(b, r, tensor=False)
833 array([-0., 0.])
834 """
835 r = np.array(r, ndmin=1, copy=False)
836 if r.dtype.char in '?bBhHiIlLqQpP':
837 r = r.astype(np.double)
838 if isinstance(x, (tuple, list)):
839 x = np.asarray(x)
840 if isinstance(x, np.ndarray):
841 if tensor:
842 r = r.reshape(r.shape + (1,)*x.ndim)
843 elif x.ndim >= r.ndim:
844 raise ValueError("x.ndim must be < r.ndim when tensor == False")
845 return np.prod(x - r, axis=0)
846
847
848def polyval2d(x, y, c):
849 """
850 Evaluate a 2-D polynomial at points (x, y).
851
852 This function returns the value
853
854 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * x^i * y^j
855
856 The parameters `x` and `y` are converted to arrays only if they are
857 tuples or a lists, otherwise they are treated as a scalars and they
858 must have the same shape after conversion. In either case, either `x`
859 and `y` or their elements must support multiplication and addition both
860 with themselves and with the elements of `c`.
861
862 If `c` has fewer than two dimensions, ones are implicitly appended to
863 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
864 x.shape.
865
866 Parameters
867 ----------
868 x, y : array_like, compatible objects
869 The two dimensional series is evaluated at the points `(x, y)`,
870 where `x` and `y` must have the same shape. If `x` or `y` is a list
871 or tuple, it is first converted to an ndarray, otherwise it is left
872 unchanged and, if it isn't an ndarray, it is treated as a scalar.
873 c : array_like
874 Array of coefficients ordered so that the coefficient of the term
875 of multi-degree i,j is contained in `c[i,j]`. If `c` has
876 dimension greater than two the remaining indices enumerate multiple
877 sets of coefficients.
878
879 Returns
880 -------
881 values : ndarray, compatible object
882 The values of the two dimensional polynomial at points formed with
883 pairs of corresponding values from `x` and `y`.
884
885 See Also
886 --------
887 polyval, polygrid2d, polyval3d, polygrid3d
888
889 Notes
890 -----
891
892 .. versionadded:: 1.7.0
893
894 """
895 return pu._valnd(polyval, c, x, y)
896
897
898def polygrid2d(x, y, c):
899 """
900 Evaluate a 2-D polynomial on the Cartesian product of x and y.
901
902 This function returns the values:
903
904 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * a^i * b^j
905
906 where the points `(a, b)` consist of all pairs formed by taking
907 `a` from `x` and `b` from `y`. The resulting points form a grid with
908 `x` in the first dimension and `y` in the second.
909
910 The parameters `x` and `y` are converted to arrays only if they are
911 tuples or a lists, otherwise they are treated as a scalars. In either
912 case, either `x` and `y` or their elements must support multiplication
913 and addition both with themselves and with the elements of `c`.
914
915 If `c` has fewer than two dimensions, ones are implicitly appended to
916 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
917 x.shape + y.shape.
918
919 Parameters
920 ----------
921 x, y : array_like, compatible objects
922 The two dimensional series is evaluated at the points in the
923 Cartesian product of `x` and `y`. If `x` or `y` is a list or
924 tuple, it is first converted to an ndarray, otherwise it is left
925 unchanged and, if it isn't an ndarray, it is treated as a scalar.
926 c : array_like
927 Array of coefficients ordered so that the coefficients for terms of
928 degree i,j are contained in ``c[i,j]``. If `c` has dimension
929 greater than two the remaining indices enumerate multiple sets of
930 coefficients.
931
932 Returns
933 -------
934 values : ndarray, compatible object
935 The values of the two dimensional polynomial at points in the Cartesian
936 product of `x` and `y`.
937
938 See Also
939 --------
940 polyval, polyval2d, polyval3d, polygrid3d
941
942 Notes
943 -----
944
945 .. versionadded:: 1.7.0
946
947 """
948 return pu._gridnd(polyval, c, x, y)
949
950
951def polyval3d(x, y, z, c):
952 """
953 Evaluate a 3-D polynomial at points (x, y, z).
954
955 This function returns the values:
956
957 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * x^i * y^j * z^k
958
959 The parameters `x`, `y`, and `z` are converted to arrays only if
960 they are tuples or a lists, otherwise they are treated as a scalars and
961 they must have the same shape after conversion. In either case, either
962 `x`, `y`, and `z` or their elements must support multiplication and
963 addition both with themselves and with the elements of `c`.
964
965 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
966 shape to make it 3-D. The shape of the result will be c.shape[3:] +
967 x.shape.
968
969 Parameters
970 ----------
971 x, y, z : array_like, compatible object
972 The three dimensional series is evaluated at the points
973 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
974 any of `x`, `y`, or `z` is a list or tuple, it is first converted
975 to an ndarray, otherwise it is left unchanged and if it isn't an
976 ndarray it is treated as a scalar.
977 c : array_like
978 Array of coefficients ordered so that the coefficient of the term of
979 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
980 greater than 3 the remaining indices enumerate multiple sets of
981 coefficients.
982
983 Returns
984 -------
985 values : ndarray, compatible object
986 The values of the multidimensional polynomial on points formed with
987 triples of corresponding values from `x`, `y`, and `z`.
988
989 See Also
990 --------
991 polyval, polyval2d, polygrid2d, polygrid3d
992
993 Notes
994 -----
995
996 .. versionadded:: 1.7.0
997
998 """
999 return pu._valnd(polyval, c, x, y, z)
1000
1001
1002def polygrid3d(x, y, z, c):
1003 """
1004 Evaluate a 3-D polynomial on the Cartesian product of x, y and z.
1005
1006 This function returns the values:
1007
1008 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * a^i * b^j * c^k
1009
1010 where the points `(a, b, c)` consist of all triples formed by taking
1011 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1012 a grid with `x` in the first dimension, `y` in the second, and `z` in
1013 the third.
1014
1015 The parameters `x`, `y`, and `z` are converted to arrays only if they
1016 are tuples or a lists, otherwise they are treated as a scalars. In
1017 either case, either `x`, `y`, and `z` or their elements must support
1018 multiplication and addition both with themselves and with the elements
1019 of `c`.
1020
1021 If `c` has fewer than three dimensions, ones are implicitly appended to
1022 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1023 x.shape + y.shape + z.shape.
1024
1025 Parameters
1026 ----------
1027 x, y, z : array_like, compatible objects
1028 The three dimensional series is evaluated at the points in the
1029 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
1030 list or tuple, it is first converted to an ndarray, otherwise it is
1031 left unchanged and, if it isn't an ndarray, it is treated as a
1032 scalar.
1033 c : array_like
1034 Array of coefficients ordered so that the coefficients for terms of
1035 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1036 greater than two the remaining indices enumerate multiple sets of
1037 coefficients.
1038
1039 Returns
1040 -------
1041 values : ndarray, compatible object
1042 The values of the two dimensional polynomial at points in the Cartesian
1043 product of `x` and `y`.
1044
1045 See Also
1046 --------
1047 polyval, polyval2d, polygrid2d, polyval3d
1048
1049 Notes
1050 -----
1051
1052 .. versionadded:: 1.7.0
1053
1054 """
1055 return pu._gridnd(polyval, c, x, y, z)
1056
1057
1058def polyvander(x, deg):
1059 """Vandermonde matrix of given degree.
1060
1061 Returns the Vandermonde matrix of degree `deg` and sample points
1062 `x`. The Vandermonde matrix is defined by
1063
1064 .. math:: V[..., i] = x^i,
1065
1066 where `0 <= i <= deg`. The leading indices of `V` index the elements of
1067 `x` and the last index is the power of `x`.
1068
1069 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
1070 matrix ``V = polyvander(x, n)``, then ``np.dot(V, c)`` and
1071 ``polyval(x, c)`` are the same up to roundoff. This equivalence is
1072 useful both for least squares fitting and for the evaluation of a large
1073 number of polynomials of the same degree and sample points.
1074
1075 Parameters
1076 ----------
1077 x : array_like
1078 Array of points. The dtype is converted to float64 or complex128
1079 depending on whether any of the elements are complex. If `x` is
1080 scalar it is converted to a 1-D array.
1081 deg : int
1082 Degree of the resulting matrix.
1083
1084 Returns
1085 -------
1086 vander : ndarray.
1087 The Vandermonde matrix. The shape of the returned matrix is
1088 ``x.shape + (deg + 1,)``, where the last index is the power of `x`.
1089 The dtype will be the same as the converted `x`.
1090
1091 See Also
1092 --------
1093 polyvander2d, polyvander3d
1094
1095 """
1096 ideg = pu._deprecate_as_int(deg, "deg")
1097 if ideg < 0:
1098 raise ValueError("deg must be non-negative")
1099
1100 x = np.array(x, copy=False, ndmin=1) + 0.0
1101 dims = (ideg + 1,) + x.shape
1102 dtyp = x.dtype
1103 v = np.empty(dims, dtype=dtyp)
1104 v[0] = x*0 + 1
1105 if ideg > 0:
1106 v[1] = x
1107 for i in range(2, ideg + 1):
1108 v[i] = v[i-1]*x
1109 return np.moveaxis(v, 0, -1)
1110
1111
1112def polyvander2d(x, y, deg):
1113 """Pseudo-Vandermonde matrix of given degrees.
1114
1115 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1116 points `(x, y)`. The pseudo-Vandermonde matrix is defined by
1117
1118 .. math:: V[..., (deg[1] + 1)*i + j] = x^i * y^j,
1119
1120 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
1121 `V` index the points `(x, y)` and the last index encodes the powers of
1122 `x` and `y`.
1123
1124 If ``V = polyvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1125 correspond to the elements of a 2-D coefficient array `c` of shape
1126 (xdeg + 1, ydeg + 1) in the order
1127
1128 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
1129
1130 and ``np.dot(V, c.flat)`` and ``polyval2d(x, y, c)`` will be the same
1131 up to roundoff. This equivalence is useful both for least squares
1132 fitting and for the evaluation of a large number of 2-D polynomials
1133 of the same degrees and sample points.
1134
1135 Parameters
1136 ----------
1137 x, y : array_like
1138 Arrays of point coordinates, all of the same shape. The dtypes
1139 will be converted to either float64 or complex128 depending on
1140 whether any of the elements are complex. Scalars are converted to
1141 1-D arrays.
1142 deg : list of ints
1143 List of maximum degrees of the form [x_deg, y_deg].
1144
1145 Returns
1146 -------
1147 vander2d : ndarray
1148 The shape of the returned matrix is ``x.shape + (order,)``, where
1149 :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same
1150 as the converted `x` and `y`.
1151
1152 See Also
1153 --------
1154 polyvander, polyvander3d, polyval2d, polyval3d
1155
1156 """
1157 return pu._vander_nd_flat((polyvander, polyvander), (x, y), deg)
1158
1159
1160def polyvander3d(x, y, z, deg):
1161 """Pseudo-Vandermonde matrix of given degrees.
1162
1163 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1164 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
1165 then The pseudo-Vandermonde matrix is defined by
1166
1167 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = x^i * y^j * z^k,
1168
1169 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
1170 indices of `V` index the points `(x, y, z)` and the last index encodes
1171 the powers of `x`, `y`, and `z`.
1172
1173 If ``V = polyvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
1174 of `V` correspond to the elements of a 3-D coefficient array `c` of
1175 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
1176
1177 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
1178
1179 and ``np.dot(V, c.flat)`` and ``polyval3d(x, y, z, c)`` will be the
1180 same up to roundoff. This equivalence is useful both for least squares
1181 fitting and for the evaluation of a large number of 3-D polynomials
1182 of the same degrees and sample points.
1183
1184 Parameters
1185 ----------
1186 x, y, z : array_like
1187 Arrays of point coordinates, all of the same shape. The dtypes will
1188 be converted to either float64 or complex128 depending on whether
1189 any of the elements are complex. Scalars are converted to 1-D
1190 arrays.
1191 deg : list of ints
1192 List of maximum degrees of the form [x_deg, y_deg, z_deg].
1193
1194 Returns
1195 -------
1196 vander3d : ndarray
1197 The shape of the returned matrix is ``x.shape + (order,)``, where
1198 :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will
1199 be the same as the converted `x`, `y`, and `z`.
1200
1201 See Also
1202 --------
1203 polyvander, polyvander3d, polyval2d, polyval3d
1204
1205 Notes
1206 -----
1207
1208 .. versionadded:: 1.7.0
1209
1210 """
1211 return pu._vander_nd_flat((polyvander, polyvander, polyvander), (x, y, z), deg)
1212
1213
1214def polyfit(x, y, deg, rcond=None, full=False, w=None):
1215 """
1216 Least-squares fit of a polynomial to data.
1217
1218 Return the coefficients of a polynomial of degree `deg` that is the
1219 least squares fit to the data values `y` given at points `x`. If `y` is
1220 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
1221 fits are done, one for each column of `y`, and the resulting
1222 coefficients are stored in the corresponding columns of a 2-D return.
1223 The fitted polynomial(s) are in the form
1224
1225 .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n,
1226
1227 where `n` is `deg`.
1228
1229 Parameters
1230 ----------
1231 x : array_like, shape (`M`,)
1232 x-coordinates of the `M` sample (data) points ``(x[i], y[i])``.
1233 y : array_like, shape (`M`,) or (`M`, `K`)
1234 y-coordinates of the sample points. Several sets of sample points
1235 sharing the same x-coordinates can be (independently) fit with one
1236 call to `polyfit` by passing in for `y` a 2-D array that contains
1237 one data set per column.
1238 deg : int or 1-D array_like
1239 Degree(s) of the fitting polynomials. If `deg` is a single integer
1240 all terms up to and including the `deg`'th term are included in the
1241 fit. For NumPy versions >= 1.11.0 a list of integers specifying the
1242 degrees of the terms to include may be used instead.
1243 rcond : float, optional
1244 Relative condition number of the fit. Singular values smaller
1245 than `rcond`, relative to the largest singular value, will be
1246 ignored. The default value is ``len(x)*eps``, where `eps` is the
1247 relative precision of the platform's float type, about 2e-16 in
1248 most cases.
1249 full : bool, optional
1250 Switch determining the nature of the return value. When ``False``
1251 (the default) just the coefficients are returned; when ``True``,
1252 diagnostic information from the singular value decomposition (used
1253 to solve the fit's matrix equation) is also returned.
1254 w : array_like, shape (`M`,), optional
1255 Weights. If not None, the weight ``w[i]`` applies to the unsquared
1256 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
1257 chosen so that the errors of the products ``w[i]*y[i]`` all have the
1258 same variance. When using inverse-variance weighting, use
1259 ``w[i] = 1/sigma(y[i])``. The default value is None.
1260
1261 .. versionadded:: 1.5.0
1262
1263 Returns
1264 -------
1265 coef : ndarray, shape (`deg` + 1,) or (`deg` + 1, `K`)
1266 Polynomial coefficients ordered from low to high. If `y` was 2-D,
1267 the coefficients in column `k` of `coef` represent the polynomial
1268 fit to the data in `y`'s `k`-th column.
1269
1270 [residuals, rank, singular_values, rcond] : list
1271 These values are only returned if ``full == True``
1272
1273 - residuals -- sum of squared residuals of the least squares fit
1274 - rank -- the numerical rank of the scaled Vandermonde matrix
1275 - singular_values -- singular values of the scaled Vandermonde matrix
1276 - rcond -- value of `rcond`.
1277
1278 For more details, see `numpy.linalg.lstsq`.
1279
1280 Raises
1281 ------
1282 RankWarning
1283 Raised if the matrix in the least-squares fit is rank deficient.
1284 The warning is only raised if ``full == False``. The warnings can
1285 be turned off by:
1286
1287 >>> import warnings
1288 >>> warnings.simplefilter('ignore', np.RankWarning)
1289
1290 See Also
1291 --------
1292 numpy.polynomial.chebyshev.chebfit
1293 numpy.polynomial.legendre.legfit
1294 numpy.polynomial.laguerre.lagfit
1295 numpy.polynomial.hermite.hermfit
1296 numpy.polynomial.hermite_e.hermefit
1297 polyval : Evaluates a polynomial.
1298 polyvander : Vandermonde matrix for powers.
1299 numpy.linalg.lstsq : Computes a least-squares fit from the matrix.
1300 scipy.interpolate.UnivariateSpline : Computes spline fits.
1301
1302 Notes
1303 -----
1304 The solution is the coefficients of the polynomial `p` that minimizes
1305 the sum of the weighted squared errors
1306
1307 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
1308
1309 where the :math:`w_j` are the weights. This problem is solved by
1310 setting up the (typically) over-determined matrix equation:
1311
1312 .. math:: V(x) * c = w * y,
1313
1314 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
1315 coefficients to be solved for, `w` are the weights, and `y` are the
1316 observed values. This equation is then solved using the singular value
1317 decomposition of `V`.
1318
1319 If some of the singular values of `V` are so small that they are
1320 neglected (and `full` == ``False``), a `RankWarning` will be raised.
1321 This means that the coefficient values may be poorly determined.
1322 Fitting to a lower order polynomial will usually get rid of the warning
1323 (but may not be what you want, of course; if you have independent
1324 reason(s) for choosing the degree which isn't working, you may have to:
1325 a) reconsider those reasons, and/or b) reconsider the quality of your
1326 data). The `rcond` parameter can also be set to a value smaller than
1327 its default, but the resulting fit may be spurious and have large
1328 contributions from roundoff error.
1329
1330 Polynomial fits using double precision tend to "fail" at about
1331 (polynomial) degree 20. Fits using Chebyshev or Legendre series are
1332 generally better conditioned, but much can still depend on the
1333 distribution of the sample points and the smoothness of the data. If
1334 the quality of the fit is inadequate, splines may be a good
1335 alternative.
1336
1337 Examples
1338 --------
1339 >>> np.random.seed(123)
1340 >>> from numpy.polynomial import polynomial as P
1341 >>> x = np.linspace(-1,1,51) # x "data": [-1, -0.96, ..., 0.96, 1]
1342 >>> y = x**3 - x + np.random.randn(len(x)) # x^3 - x + Gaussian noise
1343 >>> c, stats = P.polyfit(x,y,3,full=True)
1344 >>> np.random.seed(123)
1345 >>> c # c[0], c[2] should be approx. 0, c[1] approx. -1, c[3] approx. 1
1346 array([ 0.01909725, -1.30598256, -0.00577963, 1.02644286]) # may vary
1347 >>> stats # note the large SSR, explaining the rather poor results
1348 [array([ 38.06116253]), 4, array([ 1.38446749, 1.32119158, 0.50443316, # may vary
1349 0.28853036]), 1.1324274851176597e-014]
1350
1351 Same thing without the added noise
1352
1353 >>> y = x**3 - x
1354 >>> c, stats = P.polyfit(x,y,3,full=True)
1355 >>> c # c[0], c[2] should be "very close to 0", c[1] ~= -1, c[3] ~= 1
1356 array([-6.36925336e-18, -1.00000000e+00, -4.08053781e-16, 1.00000000e+00])
1357 >>> stats # note the minuscule SSR
1358 [array([ 7.46346754e-31]), 4, array([ 1.38446749, 1.32119158, # may vary
1359 0.50443316, 0.28853036]), 1.1324274851176597e-014]
1360
1361 """
1362 return pu._fit(polyvander, x, y, deg, rcond, full, w)
1363
1364
1365def polycompanion(c):
1366 """
1367 Return the companion matrix of c.
1368
1369 The companion matrix for power series cannot be made symmetric by
1370 scaling the basis, so this function differs from those for the
1371 orthogonal polynomials.
1372
1373 Parameters
1374 ----------
1375 c : array_like
1376 1-D array of polynomial coefficients ordered from low to high
1377 degree.
1378
1379 Returns
1380 -------
1381 mat : ndarray
1382 Companion matrix of dimensions (deg, deg).
1383
1384 Notes
1385 -----
1386
1387 .. versionadded:: 1.7.0
1388
1389 """
1390 # c is a trimmed copy
1391 [c] = pu.as_series([c])
1392 if len(c) < 2:
1393 raise ValueError('Series must have maximum degree of at least 1.')
1394 if len(c) == 2:
1395 return np.array([[-c[0]/c[1]]])
1396
1397 n = len(c) - 1
1398 mat = np.zeros((n, n), dtype=c.dtype)
1399 bot = mat.reshape(-1)[n::n+1]
1400 bot[...] = 1
1401 mat[:, -1] -= c[:-1]/c[-1]
1402 return mat
1403
1404
1405def polyroots(c):
1406 """
1407 Compute the roots of a polynomial.
1408
1409 Return the roots (a.k.a. "zeros") of the polynomial
1410
1411 .. math:: p(x) = \\sum_i c[i] * x^i.
1412
1413 Parameters
1414 ----------
1415 c : 1-D array_like
1416 1-D array of polynomial coefficients.
1417
1418 Returns
1419 -------
1420 out : ndarray
1421 Array of the roots of the polynomial. If all the roots are real,
1422 then `out` is also real, otherwise it is complex.
1423
1424 See Also
1425 --------
1426 numpy.polynomial.chebyshev.chebroots
1427 numpy.polynomial.legendre.legroots
1428 numpy.polynomial.laguerre.lagroots
1429 numpy.polynomial.hermite.hermroots
1430 numpy.polynomial.hermite_e.hermeroots
1431
1432 Notes
1433 -----
1434 The root estimates are obtained as the eigenvalues of the companion
1435 matrix, Roots far from the origin of the complex plane may have large
1436 errors due to the numerical instability of the power series for such
1437 values. Roots with multiplicity greater than 1 will also show larger
1438 errors as the value of the series near such points is relatively
1439 insensitive to errors in the roots. Isolated roots near the origin can
1440 be improved by a few iterations of Newton's method.
1441
1442 Examples
1443 --------
1444 >>> import numpy.polynomial.polynomial as poly
1445 >>> poly.polyroots(poly.polyfromroots((-1,0,1)))
1446 array([-1., 0., 1.])
1447 >>> poly.polyroots(poly.polyfromroots((-1,0,1))).dtype
1448 dtype('float64')
1449 >>> j = complex(0,1)
1450 >>> poly.polyroots(poly.polyfromroots((-j,0,j)))
1451 array([ 0.00000000e+00+0.j, 0.00000000e+00+1.j, 2.77555756e-17-1.j]) # may vary
1452
1453 """
1454 # c is a trimmed copy
1455 [c] = pu.as_series([c])
1456 if len(c) < 2:
1457 return np.array([], dtype=c.dtype)
1458 if len(c) == 2:
1459 return np.array([-c[0]/c[1]])
1460
1461 # rotated companion matrix reduces error
1462 m = polycompanion(c)[::-1,::-1]
1463 r = la.eigvals(m)
1464 r.sort()
1465 return r
1466
1467
1468#
1469# polynomial class
1470#
1471
1472class Polynomial(ABCPolyBase):
1473 """A power series class.
1474
1475 The Polynomial class provides the standard Python numerical methods
1476 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
1477 attributes and methods listed in the `ABCPolyBase` documentation.
1478
1479 Parameters
1480 ----------
1481 coef : array_like
1482 Polynomial coefficients in order of increasing degree, i.e.,
1483 ``(1, 2, 3)`` give ``1 + 2*x + 3*x**2``.
1484 domain : (2,) array_like, optional
1485 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
1486 to the interval ``[window[0], window[1]]`` by shifting and scaling.
1487 The default value is [-1, 1].
1488 window : (2,) array_like, optional
1489 Window, see `domain` for its use. The default value is [-1, 1].
1490
1491 .. versionadded:: 1.6.0
1492
1493 """
1494 # Virtual Functions
1495 _add = staticmethod(polyadd)
1496 _sub = staticmethod(polysub)
1497 _mul = staticmethod(polymul)
1498 _div = staticmethod(polydiv)
1499 _pow = staticmethod(polypow)
1500 _val = staticmethod(polyval)
1501 _int = staticmethod(polyint)
1502 _der = staticmethod(polyder)
1503 _fit = staticmethod(polyfit)
1504 _line = staticmethod(polyline)
1505 _roots = staticmethod(polyroots)
1506 _fromroots = staticmethod(polyfromroots)
1507
1508 # Virtual properties
1509 domain = np.array(polydomain)
1510 window = np.array(polydomain)
1511 basis_name = None
1512
1513 @classmethod
1514 def _str_term_unicode(cls, i, arg_str):
1515 if i == '1':
1516 return f"·{arg_str}"
1517 else:
1518 return f"·{arg_str}{i.translate(cls._superscript_mapping)}"
1519
1520 @staticmethod
1521 def _str_term_ascii(i, arg_str):
1522 if i == '1':
1523 return f" {arg_str}"
1524 else:
1525 return f" {arg_str}**{i}"
1526
1527 @staticmethod
1528 def _repr_latex_term(i, arg_str, needs_parens):
1529 if needs_parens:
1530 arg_str = rf"\left({arg_str}\right)"
1531 if i == 0:
1532 return '1'
1533 elif i == 1:
1534 return arg_str
1535 else:
1536 return f"{arg_str}^{{{i}}}"