Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Functions to operate on polynomials.
4"""
5__all__ = ['poly', 'roots', 'polyint', 'polyder', 'polyadd',
6 'polysub', 'polymul', 'polydiv', 'polyval', 'poly1d',
7 'polyfit']
9import functools
10import re
11import warnings
13import numpy._core.numeric as NX
14from numpy._core import (
15 abs,
16 array,
17 atleast_1d,
18 dot,
19 finfo,
20 hstack,
21 isscalar,
22 ones,
23 overrides,
24)
25from numpy._utils import set_module
26from numpy.exceptions import RankWarning
27from numpy.lib._function_base_impl import trim_zeros
28from numpy.lib._twodim_base_impl import diag, vander
29from numpy.lib._type_check_impl import imag, iscomplex, mintypecode, real
30from numpy.linalg import eigvals, inv, lstsq
32array_function_dispatch = functools.partial(
33 overrides.array_function_dispatch, module='numpy')
36def _poly_dispatcher(seq_of_zeros):
37 return seq_of_zeros
40@array_function_dispatch(_poly_dispatcher)
41def poly(seq_of_zeros):
42 """
43 Find the coefficients of a polynomial with the given sequence of roots.
45 .. note::
46 This forms part of the old polynomial API. Since version 1.4, the
47 new polynomial API defined in `numpy.polynomial` is preferred.
48 A summary of the differences can be found in the
49 :doc:`transition guide </reference/routines.polynomials>`.
51 Returns the coefficients of the polynomial whose leading coefficient
52 is one for the given sequence of zeros (multiple roots must be included
53 in the sequence as many times as their multiplicity; see Examples).
54 A square matrix (or array, which will be treated as a matrix) can also
55 be given, in which case the coefficients of the characteristic polynomial
56 of the matrix are returned.
58 Parameters
59 ----------
60 seq_of_zeros : array_like, shape (N,) or (N, N)
61 A sequence of polynomial roots, or a square array or matrix object.
63 Returns
64 -------
65 c : ndarray
66 1D array of polynomial coefficients from highest to lowest degree:
68 ``c[0] * x**(N) + c[1] * x**(N-1) + ... + c[N-1] * x + c[N]``
69 where c[0] always equals 1.
71 Raises
72 ------
73 ValueError
74 If input is the wrong shape (the input must be a 1-D or square
75 2-D array).
77 See Also
78 --------
79 polyval : Compute polynomial values.
80 roots : Return the roots of a polynomial.
81 polyfit : Least squares polynomial fit.
82 poly1d : A one-dimensional polynomial class.
84 Notes
85 -----
86 Specifying the roots of a polynomial still leaves one degree of
87 freedom, typically represented by an undetermined leading
88 coefficient. [1]_ In the case of this function, that coefficient -
89 the first one in the returned array - is always taken as one. (If
90 for some reason you have one other point, the only automatic way
91 presently to leverage that information is to use ``polyfit``.)
93 The characteristic polynomial, :math:`p_a(t)`, of an `n`-by-`n`
94 matrix **A** is given by
96 :math:`p_a(t) = \\mathrm{det}(t\\, \\mathbf{I} - \\mathbf{A})`,
98 where **I** is the `n`-by-`n` identity matrix. [2]_
100 References
101 ----------
102 .. [1] M. Sullivan and M. Sullivan, III, "Algebra and Trigonometry,
103 Enhanced With Graphing Utilities," Prentice-Hall, pg. 318, 1996.
105 .. [2] G. Strang, "Linear Algebra and Its Applications, 2nd Edition,"
106 Academic Press, pg. 182, 1980.
108 Examples
109 --------
111 Given a sequence of a polynomial's zeros:
113 >>> import numpy as np
115 >>> np.poly((0, 0, 0)) # Multiple root example
116 array([1., 0., 0., 0.])
118 The line above represents z**3 + 0*z**2 + 0*z + 0.
120 >>> np.poly((-1./2, 0, 1./2))
121 array([ 1. , 0. , -0.25, 0. ])
123 The line above represents z**3 - z/4
125 >>> np.poly((np.random.random(1)[0], 0, np.random.random(1)[0]))
126 array([ 1. , -0.77086955, 0.08618131, 0. ]) # random
128 Given a square array object:
130 >>> P = np.array([[0, 1./3], [-1./2, 0]])
131 >>> np.poly(P)
132 array([1. , 0. , 0.16666667])
134 Note how in all cases the leading coefficient is always 1.
136 """
137 seq_of_zeros = atleast_1d(seq_of_zeros)
138 sh = seq_of_zeros.shape
140 if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0:
141 seq_of_zeros = eigvals(seq_of_zeros)
142 elif len(sh) == 1:
143 dt = seq_of_zeros.dtype
144 # Let object arrays slip through, e.g. for arbitrary precision
145 if dt != object:
146 seq_of_zeros = seq_of_zeros.astype(mintypecode(dt.char))
147 else:
148 raise ValueError("input must be 1d or non-empty square 2d array.")
150 if len(seq_of_zeros) == 0:
151 return 1.0
152 dt = seq_of_zeros.dtype
153 a = ones((1,), dtype=dt)
154 for zero in seq_of_zeros:
155 a = NX.convolve(a, array([1, -zero], dtype=dt), mode='full')
157 if issubclass(a.dtype.type, NX.complexfloating):
158 # if complex roots are all complex conjugates, the roots are real.
159 roots = NX.asarray(seq_of_zeros, complex)
160 if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())):
161 a = a.real.copy()
163 return a
166def _roots_dispatcher(p):
167 return p
170@array_function_dispatch(_roots_dispatcher)
171def roots(p):
172 """
173 Return the roots of a polynomial with coefficients given in p.
175 .. note::
176 This forms part of the old polynomial API. Since version 1.4, the
177 new polynomial API defined in `numpy.polynomial` is preferred.
178 A summary of the differences can be found in the
179 :doc:`transition guide </reference/routines.polynomials>`.
181 The values in the rank-1 array `p` are coefficients of a polynomial.
182 If the length of `p` is n+1 then the polynomial is described by::
184 p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
186 Parameters
187 ----------
188 p : array_like
189 Rank-1 array of polynomial coefficients.
191 Returns
192 -------
193 out : ndarray
194 An array containing the roots of the polynomial.
196 Raises
197 ------
198 ValueError
199 When `p` cannot be converted to a rank-1 array.
201 See also
202 --------
203 poly : Find the coefficients of a polynomial with a given sequence
204 of roots.
205 polyval : Compute polynomial values.
206 polyfit : Least squares polynomial fit.
207 poly1d : A one-dimensional polynomial class.
209 Notes
210 -----
211 The algorithm relies on computing the eigenvalues of the
212 companion matrix [1]_.
214 References
215 ----------
216 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK:
217 Cambridge University Press, 1999, pp. 146-7.
219 Examples
220 --------
221 >>> import numpy as np
222 >>> coeff = [3.2, 2, 1]
223 >>> np.roots(coeff)
224 array([-0.3125+0.46351241j, -0.3125-0.46351241j])
226 """
227 # If input is scalar, this makes it an array
228 p = atleast_1d(p)
229 if p.ndim != 1:
230 raise ValueError("Input must be a rank-1 array.")
232 # find non-zero array entries
233 non_zero = NX.nonzero(NX.ravel(p))[0]
235 # Return an empty array if polynomial is all zeros
236 if len(non_zero) == 0:
237 return NX.array([])
239 # find the number of trailing zeros -- this is the number of roots at 0.
240 trailing_zeros = len(p) - non_zero[-1] - 1
242 # strip leading and trailing zeros
243 p = p[int(non_zero[0]):int(non_zero[-1]) + 1]
245 # casting: if incoming array isn't floating point, make it floating point.
246 if not issubclass(p.dtype.type, (NX.floating, NX.complexfloating)):
247 p = p.astype(float)
249 N = len(p)
250 if N > 1:
251 # build companion matrix and find its eigenvalues (the roots)
252 A = diag(NX.ones((N - 2,), p.dtype), -1)
253 A[0, :] = -p[1:] / p[0]
254 roots = eigvals(A)
255 else:
256 roots = NX.array([])
258 # tack any zeros onto the back of the array
259 roots = hstack((roots, NX.zeros(trailing_zeros, roots.dtype)))
260 return roots
263def _polyint_dispatcher(p, m=None, k=None):
264 return (p,)
267@array_function_dispatch(_polyint_dispatcher)
268def polyint(p, m=1, k=None):
269 """
270 Return an antiderivative (indefinite integral) of a polynomial.
272 .. note::
273 This forms part of the old polynomial API. Since version 1.4, the
274 new polynomial API defined in `numpy.polynomial` is preferred.
275 A summary of the differences can be found in the
276 :doc:`transition guide </reference/routines.polynomials>`.
278 The returned order `m` antiderivative `P` of polynomial `p` satisfies
279 :math:`\\frac{d^m}{dx^m}P(x) = p(x)` and is defined up to `m - 1`
280 integration constants `k`. The constants determine the low-order
281 polynomial part
283 .. math:: \\frac{k_{m-1}}{0!} x^0 + \\ldots + \\frac{k_0}{(m-1)!}x^{m-1}
285 of `P` so that :math:`P^{(j)}(0) = k_{m-j-1}`.
287 Parameters
288 ----------
289 p : array_like or poly1d
290 Polynomial to integrate.
291 A sequence is interpreted as polynomial coefficients, see `poly1d`.
292 m : int, optional
293 Order of the antiderivative. (Default: 1)
294 k : list of `m` scalars or scalar, optional
295 Integration constants. They are given in the order of integration:
296 those corresponding to highest-order terms come first.
298 If ``None`` (default), all constants are assumed to be zero.
299 If `m = 1`, a single scalar can be given instead of a list.
301 See Also
302 --------
303 polyder : derivative of a polynomial
304 poly1d.integ : equivalent method
306 Examples
307 --------
309 The defining property of the antiderivative:
311 >>> import numpy as np
313 >>> p = np.poly1d([1,1,1])
314 >>> P = np.polyint(p)
315 >>> P
316 poly1d([ 0.33333333, 0.5 , 1. , 0. ]) # may vary
317 >>> np.polyder(P) == p
318 True
320 The integration constants default to zero, but can be specified:
322 >>> P = np.polyint(p, 3)
323 >>> P(0)
324 0.0
325 >>> np.polyder(P)(0)
326 0.0
327 >>> np.polyder(P, 2)(0)
328 0.0
329 >>> P = np.polyint(p, 3, k=[6,5,3])
330 >>> P
331 poly1d([ 0.01666667, 0.04166667, 0.16666667, 3. , 5. , 3. ]) # may vary
333 Note that 3 = 6 / 2!, and that the constants are given in the order of
334 integrations. Constant of the highest-order polynomial term comes first:
336 >>> np.polyder(P, 2)(0)
337 6.0
338 >>> np.polyder(P, 1)(0)
339 5.0
340 >>> P(0)
341 3.0
343 """
344 m = int(m)
345 if m < 0:
346 raise ValueError("Order of integral must be positive (see polyder)")
347 if k is None:
348 k = NX.zeros(m, float)
349 k = atleast_1d(k)
350 if len(k) == 1 and m > 1:
351 k = k[0] * NX.ones(m, float)
352 if len(k) < m:
353 raise ValueError(
354 "k must be a scalar or a rank-1 array of length 1 or >m.")
356 truepoly = isinstance(p, poly1d)
357 p = NX.asarray(p)
358 if m == 0:
359 if truepoly:
360 return poly1d(p)
361 return p
362 else:
363 # Note: this must work also with object and integer arrays
364 y = NX.concatenate((p.__truediv__(NX.arange(len(p), 0, -1)), [k[0]]))
365 val = polyint(y, m - 1, k=k[1:])
366 if truepoly:
367 return poly1d(val)
368 return val
371def _polyder_dispatcher(p, m=None):
372 return (p,)
375@array_function_dispatch(_polyder_dispatcher)
376def polyder(p, m=1):
377 """
378 Return the derivative of the specified order of a polynomial.
380 .. note::
381 This forms part of the old polynomial API. Since version 1.4, the
382 new polynomial API defined in `numpy.polynomial` is preferred.
383 A summary of the differences can be found in the
384 :doc:`transition guide </reference/routines.polynomials>`.
386 Parameters
387 ----------
388 p : poly1d or sequence
389 Polynomial to differentiate.
390 A sequence is interpreted as polynomial coefficients, see `poly1d`.
391 m : int, optional
392 Order of differentiation (default: 1)
394 Returns
395 -------
396 der : poly1d
397 A new polynomial representing the derivative.
399 See Also
400 --------
401 polyint : Anti-derivative of a polynomial.
402 poly1d : Class for one-dimensional polynomials.
404 Examples
405 --------
407 The derivative of the polynomial :math:`x^3 + x^2 + x^1 + 1` is:
409 >>> import numpy as np
411 >>> p = np.poly1d([1,1,1,1])
412 >>> p2 = np.polyder(p)
413 >>> p2
414 poly1d([3, 2, 1])
416 which evaluates to:
418 >>> p2(2.)
419 17.0
421 We can verify this, approximating the derivative with
422 ``(f(x + h) - f(x))/h``:
424 >>> (p(2. + 0.001) - p(2.)) / 0.001
425 17.007000999997857
427 The fourth-order derivative of a 3rd-order polynomial is zero:
429 >>> np.polyder(p, 2)
430 poly1d([6, 2])
431 >>> np.polyder(p, 3)
432 poly1d([6])
433 >>> np.polyder(p, 4)
434 poly1d([0])
436 """
437 m = int(m)
438 if m < 0:
439 raise ValueError("Order of derivative must be positive (see polyint)")
441 truepoly = isinstance(p, poly1d)
442 p = NX.asarray(p)
443 n = len(p) - 1
444 y = p[:-1] * NX.arange(n, 0, -1)
445 if m == 0:
446 val = p
447 else:
448 val = polyder(y, m - 1)
449 if truepoly:
450 val = poly1d(val)
451 return val
454def _polyfit_dispatcher(x, y, deg, rcond=None, full=None, w=None, cov=None):
455 return (x, y, w)
458@array_function_dispatch(_polyfit_dispatcher)
459def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
460 """
461 Least squares polynomial fit.
463 .. note::
464 This forms part of the old polynomial API. Since version 1.4, the
465 new polynomial API defined in `numpy.polynomial` is preferred.
466 A summary of the differences can be found in the
467 :doc:`transition guide </reference/routines.polynomials>`.
469 Fit a polynomial ``p[0] * x**deg + ... + p[deg]`` of degree `deg`
470 to points `(x, y)`. Returns a vector of coefficients `p` that minimises
471 the squared error in the order `deg`, `deg-1`, ... `0`.
473 The `Polynomial.fit <numpy.polynomial.polynomial.Polynomial.fit>` class
474 method is recommended for new code as it is more stable numerically. See
475 the documentation of the method for more information.
477 Parameters
478 ----------
479 x : array_like, shape (M,)
480 x-coordinates of the M sample points ``(x[i], y[i])``.
481 y : array_like, shape (M,) or (M, K)
482 y-coordinates of the sample points. Several data sets of sample
483 points sharing the same x-coordinates can be fitted at once by
484 passing in a 2D-array that contains one dataset per column.
485 deg : int
486 Degree of the fitting polynomial
487 rcond : float, optional
488 Relative condition number of the fit. Singular values smaller than
489 this relative to the largest singular value will be ignored. The
490 default value is len(x)*eps, where eps is the relative precision of
491 the float type, about 2e-16 in most cases.
492 full : bool, optional
493 Switch determining nature of return value. When it is False (the
494 default) just the coefficients are returned, when True diagnostic
495 information from the singular value decomposition is also returned.
496 w : array_like, shape (M,), optional
497 Weights. If not None, the weight ``w[i]`` applies to the unsquared
498 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
499 chosen so that the errors of the products ``w[i]*y[i]`` all have the
500 same variance. When using inverse-variance weighting, use
501 ``w[i] = 1/sigma(y[i])``. The default value is None.
502 cov : bool or str, optional
503 If given and not `False`, return not just the estimate but also its
504 covariance matrix. By default, the covariance are scaled by
505 chi2/dof, where dof = M - (deg + 1), i.e., the weights are presumed
506 to be unreliable except in a relative sense and everything is scaled
507 such that the reduced chi2 is unity. This scaling is omitted if
508 ``cov='unscaled'``, as is relevant for the case that the weights are
509 w = 1/sigma, with sigma known to be a reliable estimate of the
510 uncertainty.
512 Returns
513 -------
514 p : ndarray, shape (deg + 1,) or (deg + 1, K)
515 Polynomial coefficients, highest power first. If `y` was 2-D, the
516 coefficients for `k`-th data set are in ``p[:,k]``.
518 residuals, rank, singular_values, rcond
519 These values are only returned if ``full == True``
521 - residuals -- sum of squared residuals of the least squares fit
522 - rank -- the effective rank of the scaled Vandermonde
523 coefficient matrix
524 - singular_values -- singular values of the scaled Vandermonde
525 coefficient matrix
526 - rcond -- value of `rcond`.
528 For more details, see `numpy.linalg.lstsq`.
530 V : ndarray, shape (deg + 1, deg + 1) or (deg + 1, deg + 1, K)
531 Present only if ``full == False`` and ``cov == True``. The covariance
532 matrix of the polynomial coefficient estimates. The diagonal of
533 this matrix are the variance estimates for each coefficient. If y
534 is a 2-D array, then the covariance matrix for the `k`-th data set
535 are in ``V[:,:,k]``
538 Warns
539 -----
540 RankWarning
541 The rank of the coefficient matrix in the least-squares fit is
542 deficient. The warning is only raised if ``full == False``.
544 The warnings can be turned off by
546 >>> import warnings
547 >>> warnings.simplefilter('ignore', np.exceptions.RankWarning)
549 See Also
550 --------
551 polyval : Compute polynomial values.
552 linalg.lstsq : Computes a least-squares fit.
553 scipy.interpolate.UnivariateSpline : Computes spline fits.
555 Notes
556 -----
557 The solution minimizes the squared error
559 .. math::
560 E = \\sum_{j=0}^k |p(x_j) - y_j|^2
562 in the equations::
564 x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0]
565 x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1]
566 ...
567 x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k]
569 The coefficient matrix of the coefficients `p` is a Vandermonde matrix.
571 `polyfit` issues a `~exceptions.RankWarning` when the least-squares fit is
572 badly conditioned. This implies that the best fit is not well-defined due
573 to numerical error. The results may be improved by lowering the polynomial
574 degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter
575 can also be set to a value smaller than its default, but the resulting
576 fit may be spurious: including contributions from the small singular
577 values can add numerical noise to the result.
579 Note that fitting polynomial coefficients is inherently badly conditioned
580 when the degree of the polynomial is large or the interval of sample points
581 is badly centered. The quality of the fit should always be checked in these
582 cases. When polynomial fits are not satisfactory, splines may be a good
583 alternative.
585 References
586 ----------
587 .. [1] Wikipedia, "Curve fitting",
588 https://en.wikipedia.org/wiki/Curve_fitting
589 .. [2] Wikipedia, "Polynomial interpolation",
590 https://en.wikipedia.org/wiki/Polynomial_interpolation
592 Examples
593 --------
594 >>> import numpy as np
595 >>> import warnings
596 >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
597 >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
598 >>> z = np.polyfit(x, y, 3)
599 >>> z
600 array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) # may vary
602 It is convenient to use `poly1d` objects for dealing with polynomials:
604 >>> p = np.poly1d(z)
605 >>> p(0.5)
606 0.6143849206349179 # may vary
607 >>> p(3.5)
608 -0.34732142857143039 # may vary
609 >>> p(10)
610 22.579365079365115 # may vary
612 High-order polynomials may oscillate wildly:
614 >>> with warnings.catch_warnings():
615 ... warnings.simplefilter('ignore', np.exceptions.RankWarning)
616 ... p30 = np.poly1d(np.polyfit(x, y, 30))
617 ...
618 >>> p30(4)
619 -0.80000000000000204 # may vary
620 >>> p30(5)
621 -0.99999999999999445 # may vary
622 >>> p30(4.5)
623 -0.10547061179440398 # may vary
625 Illustration:
627 >>> import matplotlib.pyplot as plt
628 >>> xp = np.linspace(-2, 6, 100)
629 >>> _ = plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--')
630 >>> plt.ylim(-2,2)
631 (-2, 2)
632 >>> plt.show()
634 """
635 order = int(deg) + 1
636 x = NX.asarray(x) + 0.0
637 y = NX.asarray(y) + 0.0
639 # check arguments.
640 if deg < 0:
641 raise ValueError("expected deg >= 0")
642 if x.ndim != 1:
643 raise TypeError("expected 1D vector for x")
644 if x.size == 0:
645 raise TypeError("expected non-empty vector for x")
646 if y.ndim < 1 or y.ndim > 2:
647 raise TypeError("expected 1D or 2D array for y")
648 if x.shape[0] != y.shape[0]:
649 raise TypeError("expected x and y to have same length")
651 # set rcond
652 if rcond is None:
653 rcond = len(x) * finfo(x.dtype).eps
655 # set up least squares equation for powers of x
656 lhs = vander(x, order)
657 rhs = y
659 # apply weighting
660 if w is not None:
661 w = NX.asarray(w) + 0.0
662 if w.ndim != 1:
663 raise TypeError("expected a 1-d array for weights")
664 if w.shape[0] != y.shape[0]:
665 raise TypeError("expected w and y to have the same length")
666 lhs *= w[:, NX.newaxis]
667 if rhs.ndim == 2:
668 rhs *= w[:, NX.newaxis]
669 else:
670 rhs *= w
672 # scale lhs to improve condition number and solve
673 scale = NX.sqrt((lhs * lhs).sum(axis=0))
674 lhs /= scale
675 c, resids, rank, s = lstsq(lhs, rhs, rcond)
676 c = (c.T / scale).T # broadcast scale coefficients
678 # warn on rank reduction, which indicates an ill conditioned matrix
679 if rank != order and not full:
680 msg = "Polyfit may be poorly conditioned"
681 warnings.warn(msg, RankWarning, stacklevel=2)
683 if full:
684 return c, resids, rank, s, rcond
685 elif cov:
686 Vbase = inv(dot(lhs.T, lhs))
687 Vbase /= NX.outer(scale, scale)
688 if cov == "unscaled":
689 fac = 1
690 else:
691 if len(x) <= order:
692 raise ValueError("the number of data points must exceed order "
693 "to scale the covariance matrix")
694 # note, this used to be: fac = resids / (len(x) - order - 2.0)
695 # it was decided that the "- 2" (originally justified by "Bayesian
696 # uncertainty analysis") is not what the user expects
697 # (see gh-11196 and gh-11197)
698 fac = resids / (len(x) - order)
699 if y.ndim == 1:
700 return c, Vbase * fac
701 else:
702 return c, Vbase[:, :, NX.newaxis] * fac
703 else:
704 return c
707def _polyval_dispatcher(p, x):
708 return (p, x)
711@array_function_dispatch(_polyval_dispatcher)
712def polyval(p, x):
713 """
714 Evaluate a polynomial at specific values.
716 .. note::
717 This forms part of the old polynomial API. Since version 1.4, the
718 new polynomial API defined in `numpy.polynomial` is preferred.
719 A summary of the differences can be found in the
720 :doc:`transition guide </reference/routines.polynomials>`.
722 If `p` is of length N, this function returns the value::
724 p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]
726 If `x` is a sequence, then ``p(x)`` is returned for each element of ``x``.
727 If `x` is another polynomial then the composite polynomial ``p(x(t))``
728 is returned.
730 Parameters
731 ----------
732 p : array_like or poly1d object
733 1D array of polynomial coefficients (including coefficients equal
734 to zero) from highest degree to the constant term, or an
735 instance of poly1d.
736 x : array_like or poly1d object
737 A number, an array of numbers, or an instance of poly1d, at
738 which to evaluate `p`.
740 Returns
741 -------
742 values : ndarray or poly1d
743 If `x` is a poly1d instance, the result is the composition of the two
744 polynomials, i.e., `x` is "substituted" in `p` and the simplified
745 result is returned. In addition, the type of `x` - array_like or
746 poly1d - governs the type of the output: `x` array_like => `values`
747 array_like, `x` a poly1d object => `values` is also.
749 See Also
750 --------
751 poly1d: A polynomial class.
753 Notes
754 -----
755 Horner's scheme [1]_ is used to evaluate the polynomial. Even so,
756 for polynomials of high degree the values may be inaccurate due to
757 rounding errors. Use carefully.
759 If `x` is a subtype of `ndarray` the return value will be of the same type.
761 References
762 ----------
763 .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng.
764 trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand
765 Reinhold Co., 1985, pg. 720.
767 Examples
768 --------
769 >>> import numpy as np
770 >>> np.polyval([3,0,1], 5) # 3 * 5**2 + 0 * 5**1 + 1
771 76
772 >>> np.polyval([3,0,1], np.poly1d(5))
773 poly1d([76])
774 >>> np.polyval(np.poly1d([3,0,1]), 5)
775 76
776 >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5))
777 poly1d([76])
779 """
780 p = NX.asarray(p)
781 if isinstance(x, poly1d):
782 y = 0
783 else:
784 x = NX.asanyarray(x)
785 y = NX.zeros_like(x)
786 for pv in p:
787 y = y * x + pv
788 return y
791def _binary_op_dispatcher(a1, a2):
792 return (a1, a2)
795@array_function_dispatch(_binary_op_dispatcher)
796def polyadd(a1, a2):
797 """
798 Find the sum of two polynomials.
800 .. note::
801 This forms part of the old polynomial API. Since version 1.4, the
802 new polynomial API defined in `numpy.polynomial` is preferred.
803 A summary of the differences can be found in the
804 :doc:`transition guide </reference/routines.polynomials>`.
806 Returns the polynomial resulting from the sum of two input polynomials.
807 Each input must be either a poly1d object or a 1D sequence of polynomial
808 coefficients, from highest to lowest degree.
810 Parameters
811 ----------
812 a1, a2 : array_like or poly1d object
813 Input polynomials.
815 Returns
816 -------
817 out : ndarray or poly1d object
818 The sum of the inputs. If either input is a poly1d object, then the
819 output is also a poly1d object. Otherwise, it is a 1D array of
820 polynomial coefficients from highest to lowest degree.
822 See Also
823 --------
824 poly1d : A one-dimensional polynomial class.
825 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
827 Examples
828 --------
829 >>> import numpy as np
830 >>> np.polyadd([1, 2], [9, 5, 4])
831 array([9, 6, 6])
833 Using poly1d objects:
835 >>> p1 = np.poly1d([1, 2])
836 >>> p2 = np.poly1d([9, 5, 4])
837 >>> print(p1)
838 1 x + 2
839 >>> print(p2)
840 2
841 9 x + 5 x + 4
842 >>> print(np.polyadd(p1, p2))
843 2
844 9 x + 6 x + 6
846 """
847 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
848 a1 = atleast_1d(a1)
849 a2 = atleast_1d(a2)
850 diff = len(a2) - len(a1)
851 if diff == 0:
852 val = a1 + a2
853 elif diff > 0:
854 zr = NX.zeros(diff, a1.dtype)
855 val = NX.concatenate((zr, a1)) + a2
856 else:
857 zr = NX.zeros(abs(diff), a2.dtype)
858 val = a1 + NX.concatenate((zr, a2))
859 if truepoly:
860 val = poly1d(val)
861 return val
864@array_function_dispatch(_binary_op_dispatcher)
865def polysub(a1, a2):
866 """
867 Difference (subtraction) of two polynomials.
869 .. note::
870 This forms part of the old polynomial API. Since version 1.4, the
871 new polynomial API defined in `numpy.polynomial` is preferred.
872 A summary of the differences can be found in the
873 :doc:`transition guide </reference/routines.polynomials>`.
875 Given two polynomials `a1` and `a2`, returns ``a1 - a2``.
876 `a1` and `a2` can be either array_like sequences of the polynomials'
877 coefficients (including coefficients equal to zero), or `poly1d` objects.
879 Parameters
880 ----------
881 a1, a2 : array_like or poly1d
882 Minuend and subtrahend polynomials, respectively.
884 Returns
885 -------
886 out : ndarray or poly1d
887 Array or `poly1d` object of the difference polynomial's coefficients.
889 See Also
890 --------
891 polyval, polydiv, polymul, polyadd
893 Examples
894 --------
896 .. math:: (2 x^2 + 10 x - 2) - (3 x^2 + 10 x -4) = (-x^2 + 2)
898 >>> import numpy as np
900 >>> np.polysub([2, 10, -2], [3, 10, -4])
901 array([-1, 0, 2])
903 """
904 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
905 a1 = atleast_1d(a1)
906 a2 = atleast_1d(a2)
907 diff = len(a2) - len(a1)
908 if diff == 0:
909 val = a1 - a2
910 elif diff > 0:
911 zr = NX.zeros(diff, a1.dtype)
912 val = NX.concatenate((zr, a1)) - a2
913 else:
914 zr = NX.zeros(abs(diff), a2.dtype)
915 val = a1 - NX.concatenate((zr, a2))
916 if truepoly:
917 val = poly1d(val)
918 return val
921@array_function_dispatch(_binary_op_dispatcher)
922def polymul(a1, a2):
923 """
924 Find the product of two polynomials.
926 .. note::
927 This forms part of the old polynomial API. Since version 1.4, the
928 new polynomial API defined in `numpy.polynomial` is preferred.
929 A summary of the differences can be found in the
930 :doc:`transition guide </reference/routines.polynomials>`.
932 Finds the polynomial resulting from the multiplication of the two input
933 polynomials. Each input must be either a poly1d object or a 1D sequence
934 of polynomial coefficients, from highest to lowest degree.
936 Parameters
937 ----------
938 a1, a2 : array_like or poly1d object
939 Input polynomials.
941 Returns
942 -------
943 out : ndarray or poly1d object
944 The polynomial resulting from the multiplication of the inputs. If
945 either inputs is a poly1d object, then the output is also a poly1d
946 object. Otherwise, it is a 1D array of polynomial coefficients from
947 highest to lowest degree.
949 See Also
950 --------
951 poly1d : A one-dimensional polynomial class.
952 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
953 convolve : Array convolution. Same output as polymul, but has parameter
954 for overlap mode.
956 Examples
957 --------
958 >>> import numpy as np
959 >>> np.polymul([1, 2, 3], [9, 5, 1])
960 array([ 9, 23, 38, 17, 3])
962 Using poly1d objects:
964 >>> p1 = np.poly1d([1, 2, 3])
965 >>> p2 = np.poly1d([9, 5, 1])
966 >>> print(p1)
967 2
968 1 x + 2 x + 3
969 >>> print(p2)
970 2
971 9 x + 5 x + 1
972 >>> print(np.polymul(p1, p2))
973 4 3 2
974 9 x + 23 x + 38 x + 17 x + 3
976 """
977 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
978 a1, a2 = poly1d(a1), poly1d(a2)
979 val = NX.convolve(a1, a2)
980 if truepoly:
981 val = poly1d(val)
982 return val
985def _polydiv_dispatcher(u, v):
986 return (u, v)
989@array_function_dispatch(_polydiv_dispatcher)
990def polydiv(u, v):
991 """
992 Returns the quotient and remainder of polynomial division.
994 .. note::
995 This forms part of the old polynomial API. Since version 1.4, the
996 new polynomial API defined in `numpy.polynomial` is preferred.
997 A summary of the differences can be found in the
998 :doc:`transition guide </reference/routines.polynomials>`.
1000 The input arrays are the coefficients (including any coefficients
1001 equal to zero) of the "numerator" (dividend) and "denominator"
1002 (divisor) polynomials, respectively.
1004 Parameters
1005 ----------
1006 u : array_like or poly1d
1007 Dividend polynomial's coefficients.
1009 v : array_like or poly1d
1010 Divisor polynomial's coefficients.
1012 Returns
1013 -------
1014 q : ndarray
1015 Coefficients, including those equal to zero, of the quotient.
1016 r : ndarray
1017 Coefficients, including those equal to zero, of the remainder.
1019 See Also
1020 --------
1021 poly, polyadd, polyder, polydiv, polyfit, polyint, polymul, polysub
1022 polyval
1024 Notes
1025 -----
1026 Both `u` and `v` must be 0-d or 1-d (ndim = 0 or 1), but `u.ndim` need
1027 not equal `v.ndim`. In other words, all four possible combinations -
1028 ``u.ndim = v.ndim = 0``, ``u.ndim = v.ndim = 1``,
1029 ``u.ndim = 1, v.ndim = 0``, and ``u.ndim = 0, v.ndim = 1`` - work.
1031 Examples
1032 --------
1034 .. math:: \\frac{3x^2 + 5x + 2}{2x + 1} = 1.5x + 1.75, remainder 0.25
1036 >>> import numpy as np
1038 >>> x = np.array([3.0, 5.0, 2.0])
1039 >>> y = np.array([2.0, 1.0])
1040 >>> np.polydiv(x, y)
1041 (array([1.5 , 1.75]), array([0.25]))
1043 """
1044 truepoly = (isinstance(u, poly1d) or isinstance(v, poly1d))
1045 u = atleast_1d(u) + 0.0
1046 v = atleast_1d(v) + 0.0
1047 # w has the common type
1048 w = u[0] + v[0]
1049 m = len(u) - 1
1050 n = len(v) - 1
1051 scale = 1. / v[0]
1052 q = NX.zeros((max(m - n + 1, 1),), w.dtype)
1053 r = u.astype(w.dtype)
1054 for k in range(m - n + 1):
1055 d = scale * r[k]
1056 q[k] = d
1057 r[k:k + n + 1] -= d * v
1058 while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
1059 r = r[1:]
1060 if truepoly:
1061 return poly1d(q), poly1d(r)
1062 return q, r
1065_poly_mat = re.compile(r"\*\*([0-9]*)")
1066def _raise_power(astr, wrap=70):
1067 n = 0
1068 line1 = ''
1069 line2 = ''
1070 output = ' '
1071 while True:
1072 mat = _poly_mat.search(astr, n)
1073 if mat is None:
1074 break
1075 span = mat.span()
1076 power = mat.groups()[0]
1077 partstr = astr[n:span[0]]
1078 n = span[1]
1079 toadd2 = partstr + ' ' * (len(power) - 1)
1080 toadd1 = ' ' * (len(partstr) - 1) + power
1081 if ((len(line2) + len(toadd2) > wrap) or
1082 (len(line1) + len(toadd1) > wrap)):
1083 output += line1 + "\n" + line2 + "\n "
1084 line1 = toadd1
1085 line2 = toadd2
1086 else:
1087 line2 += partstr + ' ' * (len(power) - 1)
1088 line1 += ' ' * (len(partstr) - 1) + power
1089 output += line1 + "\n" + line2
1090 return output + astr[n:]
1093@set_module('numpy')
1094class poly1d:
1095 """
1096 A one-dimensional polynomial class.
1098 .. note::
1099 This forms part of the old polynomial API. Since version 1.4, the
1100 new polynomial API defined in `numpy.polynomial` is preferred.
1101 A summary of the differences can be found in the
1102 :doc:`transition guide </reference/routines.polynomials>`.
1104 A convenience class, used to encapsulate "natural" operations on
1105 polynomials so that said operations may take on their customary
1106 form in code (see Examples).
1108 Parameters
1109 ----------
1110 c_or_r : array_like
1111 The polynomial's coefficients, in decreasing powers, or if
1112 the value of the second parameter is True, the polynomial's
1113 roots (values where the polynomial evaluates to 0). For example,
1114 ``poly1d([1, 2, 3])`` returns an object that represents
1115 :math:`x^2 + 2x + 3`, whereas ``poly1d([1, 2, 3], True)`` returns
1116 one that represents :math:`(x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x -6`.
1117 r : bool, optional
1118 If True, `c_or_r` specifies the polynomial's roots; the default
1119 is False.
1120 variable : str, optional
1121 Changes the variable used when printing `p` from `x` to `variable`
1122 (see Examples).
1124 Examples
1125 --------
1126 >>> import numpy as np
1128 Construct the polynomial :math:`x^2 + 2x + 3`:
1130 >>> import numpy as np
1132 >>> p = np.poly1d([1, 2, 3])
1133 >>> print(np.poly1d(p))
1134 2
1135 1 x + 2 x + 3
1137 Evaluate the polynomial at :math:`x = 0.5`:
1139 >>> p(0.5)
1140 4.25
1142 Find the roots:
1144 >>> p.r
1145 array([-1.+1.41421356j, -1.-1.41421356j])
1146 >>> p(p.r)
1147 array([ -4.44089210e-16+0.j, -4.44089210e-16+0.j]) # may vary
1149 These numbers in the previous line represent (0, 0) to machine precision
1151 Show the coefficients:
1153 >>> p.c
1154 array([1, 2, 3])
1156 Display the order (the leading zero-coefficients are removed):
1158 >>> p.order
1159 2
1161 Show the coefficient of the k-th power in the polynomial
1162 (which is equivalent to ``p.c[-(i+1)]``):
1164 >>> p[1]
1165 2
1167 Polynomials can be added, subtracted, multiplied, and divided
1168 (returns quotient and remainder):
1170 >>> p * p
1171 poly1d([ 1, 4, 10, 12, 9])
1173 >>> (p**3 + 4) / p
1174 (poly1d([ 1., 4., 10., 12., 9.]), poly1d([4.]))
1176 ``asarray(p)`` gives the coefficient array, so polynomials can be
1177 used in all functions that accept arrays:
1179 >>> p**2 # square of polynomial
1180 poly1d([ 1, 4, 10, 12, 9])
1182 >>> np.square(p) # square of individual coefficients
1183 array([1, 4, 9])
1185 The variable used in the string representation of `p` can be modified,
1186 using the `variable` parameter:
1188 >>> p = np.poly1d([1,2,3], variable='z')
1189 >>> print(p)
1190 2
1191 1 z + 2 z + 3
1193 Construct a polynomial from its roots:
1195 >>> np.poly1d([1, 2], True)
1196 poly1d([ 1., -3., 2.])
1198 This is the same polynomial as obtained by:
1200 >>> np.poly1d([1, -1]) * np.poly1d([1, -2])
1201 poly1d([ 1, -3, 2])
1203 """
1204 __hash__ = None
1206 @property
1207 def coeffs(self):
1208 """ The polynomial coefficients """
1209 return self._coeffs
1211 @coeffs.setter
1212 def coeffs(self, value):
1213 # allowing this makes p.coeffs *= 2 legal
1214 if value is not self._coeffs:
1215 raise AttributeError("Cannot set attribute")
1217 @property
1218 def variable(self):
1219 """ The name of the polynomial variable """
1220 return self._variable
1222 # calculated attributes
1223 @property
1224 def order(self):
1225 """ The order or degree of the polynomial """
1226 return len(self._coeffs) - 1
1228 @property
1229 def roots(self):
1230 """ The roots of the polynomial, where self(x) == 0 """
1231 return roots(self._coeffs)
1233 # our internal _coeffs property need to be backed by __dict__['coeffs'] for
1234 # scipy to work correctly.
1235 @property
1236 def _coeffs(self):
1237 return self.__dict__['coeffs']
1239 @_coeffs.setter
1240 def _coeffs(self, coeffs):
1241 self.__dict__['coeffs'] = coeffs
1243 # alias attributes
1244 r = roots
1245 c = coef = coefficients = coeffs
1246 o = order
1248 def __init__(self, c_or_r, r=False, variable=None):
1249 if isinstance(c_or_r, poly1d):
1250 self._variable = c_or_r._variable
1251 self._coeffs = c_or_r._coeffs
1253 if set(c_or_r.__dict__) - set(self.__dict__):
1254 msg = ("In the future extra properties will not be copied "
1255 "across when constructing one poly1d from another")
1256 warnings.warn(msg, FutureWarning, stacklevel=2)
1257 self.__dict__.update(c_or_r.__dict__)
1259 if variable is not None:
1260 self._variable = variable
1261 return
1262 if r:
1263 c_or_r = poly(c_or_r)
1264 c_or_r = atleast_1d(c_or_r)
1265 if c_or_r.ndim > 1:
1266 raise ValueError("Polynomial must be 1d only.")
1267 c_or_r = trim_zeros(c_or_r, trim='f')
1268 if len(c_or_r) == 0:
1269 c_or_r = NX.array([0], dtype=c_or_r.dtype)
1270 self._coeffs = c_or_r
1271 if variable is None:
1272 variable = 'x'
1273 self._variable = variable
1275 def __array__(self, t=None, copy=None):
1276 if t:
1277 return NX.asarray(self.coeffs, t, copy=copy)
1278 else:
1279 return NX.asarray(self.coeffs, copy=copy)
1281 def __repr__(self):
1282 vals = repr(self.coeffs)
1283 vals = vals[6:-1]
1284 return f"poly1d({vals})"
1286 def __len__(self):
1287 return self.order
1289 def __str__(self):
1290 thestr = "0"
1291 var = self.variable
1293 # Remove leading zeros
1294 coeffs = self.coeffs[NX.logical_or.accumulate(self.coeffs != 0)]
1295 N = len(coeffs) - 1
1297 def fmt_float(q):
1298 s = f'{q:.4g}'
1299 s = s.removesuffix('.0000')
1300 return s
1302 for k, coeff in enumerate(coeffs):
1303 if not iscomplex(coeff):
1304 coefstr = fmt_float(real(coeff))
1305 elif real(coeff) == 0:
1306 coefstr = f'{fmt_float(imag(coeff))}j'
1307 else:
1308 coefstr = f'({fmt_float(real(coeff))} + {fmt_float(imag(coeff))}j)'
1310 power = (N - k)
1311 if power == 0:
1312 if coefstr != '0':
1313 newstr = f'{coefstr}'
1314 elif k == 0:
1315 newstr = '0'
1316 else:
1317 newstr = ''
1318 elif power == 1:
1319 if coefstr == '0':
1320 newstr = ''
1321 elif coefstr == 'b':
1322 newstr = var
1323 else:
1324 newstr = f'{coefstr} {var}'
1325 elif coefstr == '0':
1326 newstr = ''
1327 elif coefstr == 'b':
1328 newstr = '%s**%d' % (var, power,)
1329 else:
1330 newstr = '%s %s**%d' % (coefstr, var, power)
1332 if k > 0:
1333 if newstr != '':
1334 if newstr.startswith('-'):
1335 thestr = f"{thestr} - {newstr[1:]}"
1336 else:
1337 thestr = f"{thestr} + {newstr}"
1338 else:
1339 thestr = newstr
1340 return _raise_power(thestr)
1342 def __call__(self, val):
1343 return polyval(self.coeffs, val)
1345 def __neg__(self):
1346 return poly1d(-self.coeffs)
1348 def __pos__(self):
1349 return self
1351 def __mul__(self, other):
1352 if isscalar(other):
1353 return poly1d(self.coeffs * other)
1354 else:
1355 other = poly1d(other)
1356 return poly1d(polymul(self.coeffs, other.coeffs))
1358 def __rmul__(self, other):
1359 if isscalar(other):
1360 return poly1d(other * self.coeffs)
1361 else:
1362 other = poly1d(other)
1363 return poly1d(polymul(self.coeffs, other.coeffs))
1365 def __add__(self, other):
1366 other = poly1d(other)
1367 return poly1d(polyadd(self.coeffs, other.coeffs))
1369 def __radd__(self, other):
1370 other = poly1d(other)
1371 return poly1d(polyadd(self.coeffs, other.coeffs))
1373 def __pow__(self, val):
1374 if not isscalar(val) or int(val) != val or val < 0:
1375 raise ValueError("Power to non-negative integers only.")
1376 res = [1]
1377 for _ in range(val):
1378 res = polymul(self.coeffs, res)
1379 return poly1d(res)
1381 def __sub__(self, other):
1382 other = poly1d(other)
1383 return poly1d(polysub(self.coeffs, other.coeffs))
1385 def __rsub__(self, other):
1386 other = poly1d(other)
1387 return poly1d(polysub(other.coeffs, self.coeffs))
1389 def __truediv__(self, other):
1390 if isscalar(other):
1391 return poly1d(self.coeffs / other)
1392 else:
1393 other = poly1d(other)
1394 return polydiv(self, other)
1396 def __rtruediv__(self, other):
1397 if isscalar(other):
1398 return poly1d(other / self.coeffs)
1399 else:
1400 other = poly1d(other)
1401 return polydiv(other, self)
1403 def __eq__(self, other):
1404 if not isinstance(other, poly1d):
1405 return NotImplemented
1406 if self.coeffs.shape != other.coeffs.shape:
1407 return False
1408 return (self.coeffs == other.coeffs).all()
1410 def __ne__(self, other):
1411 if not isinstance(other, poly1d):
1412 return NotImplemented
1413 return not self.__eq__(other)
1415 def __getitem__(self, val):
1416 ind = self.order - val
1417 if val > self.order:
1418 return self.coeffs.dtype.type(0)
1419 if val < 0:
1420 return self.coeffs.dtype.type(0)
1421 return self.coeffs[ind]
1423 def __setitem__(self, key, val):
1424 ind = self.order - key
1425 if key < 0:
1426 raise ValueError("Does not support negative powers.")
1427 if key > self.order:
1428 zr = NX.zeros(key - self.order, self.coeffs.dtype)
1429 self._coeffs = NX.concatenate((zr, self.coeffs))
1430 ind = 0
1431 self._coeffs[ind] = val
1433 def __iter__(self):
1434 return iter(self.coeffs)
1436 def integ(self, m=1, k=0):
1437 """
1438 Return an antiderivative (indefinite integral) of this polynomial.
1440 Refer to `polyint` for full documentation.
1442 See Also
1443 --------
1444 polyint : equivalent function
1446 """
1447 return poly1d(polyint(self.coeffs, m=m, k=k))
1449 def deriv(self, m=1):
1450 """
1451 Return a derivative of this polynomial.
1453 Refer to `polyder` for full documentation.
1455 See Also
1456 --------
1457 polyder : equivalent function
1459 """
1460 return poly1d(polyder(self.coeffs, m=m))
1462# Stuff to do on module import
1465warnings.simplefilter('always', RankWarning)