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