Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/numpy/lib/_polynomial_impl.py: 21%
436 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-09 06:12 +0000
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 >>> np.poly((0, 0, 0)) # Multiple root example
107 array([1., 0., 0., 0.])
109 The line above represents z**3 + 0*z**2 + 0*z + 0.
111 >>> np.poly((-1./2, 0, 1./2))
112 array([ 1. , 0. , -0.25, 0. ])
114 The line above represents z**3 - z/4
116 >>> np.poly((np.random.random(1)[0], 0, np.random.random(1)[0]))
117 array([ 1. , -0.77086955, 0.08618131, 0. ]) # random
119 Given a square array object:
121 >>> P = np.array([[0, 1./3], [-1./2, 0]])
122 >>> np.poly(P)
123 array([1. , 0. , 0.16666667])
125 Note how in all cases the leading coefficient is always 1.
127 """
128 seq_of_zeros = atleast_1d(seq_of_zeros)
129 sh = seq_of_zeros.shape
131 if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0:
132 seq_of_zeros = eigvals(seq_of_zeros)
133 elif len(sh) == 1:
134 dt = seq_of_zeros.dtype
135 # Let object arrays slip through, e.g. for arbitrary precision
136 if dt != object:
137 seq_of_zeros = seq_of_zeros.astype(mintypecode(dt.char))
138 else:
139 raise ValueError("input must be 1d or non-empty square 2d array.")
141 if len(seq_of_zeros) == 0:
142 return 1.0
143 dt = seq_of_zeros.dtype
144 a = ones((1,), dtype=dt)
145 for zero in seq_of_zeros:
146 a = NX.convolve(a, array([1, -zero], dtype=dt), mode='full')
148 if issubclass(a.dtype.type, NX.complexfloating):
149 # if complex roots are all complex conjugates, the roots are real.
150 roots = NX.asarray(seq_of_zeros, complex)
151 if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())):
152 a = a.real.copy()
154 return a
157def _roots_dispatcher(p):
158 return p
161@array_function_dispatch(_roots_dispatcher)
162def roots(p):
163 """
164 Return the roots of a polynomial with coefficients given in p.
166 .. note::
167 This forms part of the old polynomial API. Since version 1.4, the
168 new polynomial API defined in `numpy.polynomial` is preferred.
169 A summary of the differences can be found in the
170 :doc:`transition guide </reference/routines.polynomials>`.
172 The values in the rank-1 array `p` are coefficients of a polynomial.
173 If the length of `p` is n+1 then the polynomial is described by::
175 p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
177 Parameters
178 ----------
179 p : array_like
180 Rank-1 array of polynomial coefficients.
182 Returns
183 -------
184 out : ndarray
185 An array containing the roots of the polynomial.
187 Raises
188 ------
189 ValueError
190 When `p` cannot be converted to a rank-1 array.
192 See also
193 --------
194 poly : Find the coefficients of a polynomial with a given sequence
195 of roots.
196 polyval : Compute polynomial values.
197 polyfit : Least squares polynomial fit.
198 poly1d : A one-dimensional polynomial class.
200 Notes
201 -----
202 The algorithm relies on computing the eigenvalues of the
203 companion matrix [1]_.
205 References
206 ----------
207 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK:
208 Cambridge University Press, 1999, pp. 146-7.
210 Examples
211 --------
212 >>> coeff = [3.2, 2, 1]
213 >>> np.roots(coeff)
214 array([-0.3125+0.46351241j, -0.3125-0.46351241j])
216 """
217 # If input is scalar, this makes it an array
218 p = atleast_1d(p)
219 if p.ndim != 1:
220 raise ValueError("Input must be a rank-1 array.")
222 # find non-zero array entries
223 non_zero = NX.nonzero(NX.ravel(p))[0]
225 # Return an empty array if polynomial is all zeros
226 if len(non_zero) == 0:
227 return NX.array([])
229 # find the number of trailing zeros -- this is the number of roots at 0.
230 trailing_zeros = len(p) - non_zero[-1] - 1
232 # strip leading and trailing zeros
233 p = p[int(non_zero[0]):int(non_zero[-1])+1]
235 # casting: if incoming array isn't floating point, make it floating point.
236 if not issubclass(p.dtype.type, (NX.floating, NX.complexfloating)):
237 p = p.astype(float)
239 N = len(p)
240 if N > 1:
241 # build companion matrix and find its eigenvalues (the roots)
242 A = diag(NX.ones((N-2,), p.dtype), -1)
243 A[0,:] = -p[1:] / p[0]
244 roots = eigvals(A)
245 else:
246 roots = NX.array([])
248 # tack any zeros onto the back of the array
249 roots = hstack((roots, NX.zeros(trailing_zeros, roots.dtype)))
250 return roots
253def _polyint_dispatcher(p, m=None, k=None):
254 return (p,)
257@array_function_dispatch(_polyint_dispatcher)
258def polyint(p, m=1, k=None):
259 """
260 Return an antiderivative (indefinite integral) of a polynomial.
262 .. note::
263 This forms part of the old polynomial API. Since version 1.4, the
264 new polynomial API defined in `numpy.polynomial` is preferred.
265 A summary of the differences can be found in the
266 :doc:`transition guide </reference/routines.polynomials>`.
268 The returned order `m` antiderivative `P` of polynomial `p` satisfies
269 :math:`\\frac{d^m}{dx^m}P(x) = p(x)` and is defined up to `m - 1`
270 integration constants `k`. The constants determine the low-order
271 polynomial part
273 .. math:: \\frac{k_{m-1}}{0!} x^0 + \\ldots + \\frac{k_0}{(m-1)!}x^{m-1}
275 of `P` so that :math:`P^{(j)}(0) = k_{m-j-1}`.
277 Parameters
278 ----------
279 p : array_like or poly1d
280 Polynomial to integrate.
281 A sequence is interpreted as polynomial coefficients, see `poly1d`.
282 m : int, optional
283 Order of the antiderivative. (Default: 1)
284 k : list of `m` scalars or scalar, optional
285 Integration constants. They are given in the order of integration:
286 those corresponding to highest-order terms come first.
288 If ``None`` (default), all constants are assumed to be zero.
289 If `m = 1`, a single scalar can be given instead of a list.
291 See Also
292 --------
293 polyder : derivative of a polynomial
294 poly1d.integ : equivalent method
296 Examples
297 --------
298 The defining property of the antiderivative:
300 >>> p = np.poly1d([1,1,1])
301 >>> P = np.polyint(p)
302 >>> P
303 poly1d([ 0.33333333, 0.5 , 1. , 0. ]) # may vary
304 >>> np.polyder(P) == p
305 True
307 The integration constants default to zero, but can be specified:
309 >>> P = np.polyint(p, 3)
310 >>> P(0)
311 0.0
312 >>> np.polyder(P)(0)
313 0.0
314 >>> np.polyder(P, 2)(0)
315 0.0
316 >>> P = np.polyint(p, 3, k=[6,5,3])
317 >>> P
318 poly1d([ 0.01666667, 0.04166667, 0.16666667, 3. , 5. , 3. ]) # may vary
320 Note that 3 = 6 / 2!, and that the constants are given in the order of
321 integrations. Constant of the highest-order polynomial term comes first:
323 >>> np.polyder(P, 2)(0)
324 6.0
325 >>> np.polyder(P, 1)(0)
326 5.0
327 >>> P(0)
328 3.0
330 """
331 m = int(m)
332 if m < 0:
333 raise ValueError("Order of integral must be positive (see polyder)")
334 if k is None:
335 k = NX.zeros(m, float)
336 k = atleast_1d(k)
337 if len(k) == 1 and m > 1:
338 k = k[0]*NX.ones(m, float)
339 if len(k) < m:
340 raise ValueError(
341 "k must be a scalar or a rank-1 array of length 1 or >m.")
343 truepoly = isinstance(p, poly1d)
344 p = NX.asarray(p)
345 if m == 0:
346 if truepoly:
347 return poly1d(p)
348 return p
349 else:
350 # Note: this must work also with object and integer arrays
351 y = NX.concatenate((p.__truediv__(NX.arange(len(p), 0, -1)), [k[0]]))
352 val = polyint(y, m - 1, k=k[1:])
353 if truepoly:
354 return poly1d(val)
355 return val
358def _polyder_dispatcher(p, m=None):
359 return (p,)
362@array_function_dispatch(_polyder_dispatcher)
363def polyder(p, m=1):
364 """
365 Return the derivative of the specified order of a polynomial.
367 .. note::
368 This forms part of the old polynomial API. Since version 1.4, the
369 new polynomial API defined in `numpy.polynomial` is preferred.
370 A summary of the differences can be found in the
371 :doc:`transition guide </reference/routines.polynomials>`.
373 Parameters
374 ----------
375 p : poly1d or sequence
376 Polynomial to differentiate.
377 A sequence is interpreted as polynomial coefficients, see `poly1d`.
378 m : int, optional
379 Order of differentiation (default: 1)
381 Returns
382 -------
383 der : poly1d
384 A new polynomial representing the derivative.
386 See Also
387 --------
388 polyint : Anti-derivative of a polynomial.
389 poly1d : Class for one-dimensional polynomials.
391 Examples
392 --------
393 The derivative of the polynomial :math:`x^3 + x^2 + x^1 + 1` is:
395 >>> p = np.poly1d([1,1,1,1])
396 >>> p2 = np.polyder(p)
397 >>> p2
398 poly1d([3, 2, 1])
400 which evaluates to:
402 >>> p2(2.)
403 17.0
405 We can verify this, approximating the derivative with
406 ``(f(x + h) - f(x))/h``:
408 >>> (p(2. + 0.001) - p(2.)) / 0.001
409 17.007000999997857
411 The fourth-order derivative of a 3rd-order polynomial is zero:
413 >>> np.polyder(p, 2)
414 poly1d([6, 2])
415 >>> np.polyder(p, 3)
416 poly1d([6])
417 >>> np.polyder(p, 4)
418 poly1d([0])
420 """
421 m = int(m)
422 if m < 0:
423 raise ValueError("Order of derivative must be positive (see polyint)")
425 truepoly = isinstance(p, poly1d)
426 p = NX.asarray(p)
427 n = len(p) - 1
428 y = p[:-1] * NX.arange(n, 0, -1)
429 if m == 0:
430 val = p
431 else:
432 val = polyder(y, m - 1)
433 if truepoly:
434 val = poly1d(val)
435 return val
438def _polyfit_dispatcher(x, y, deg, rcond=None, full=None, w=None, cov=None):
439 return (x, y, w)
442@array_function_dispatch(_polyfit_dispatcher)
443def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
444 """
445 Least squares polynomial fit.
447 .. note::
448 This forms part of the old polynomial API. Since version 1.4, the
449 new polynomial API defined in `numpy.polynomial` is preferred.
450 A summary of the differences can be found in the
451 :doc:`transition guide </reference/routines.polynomials>`.
453 Fit a polynomial ``p(x) = p[0] * x**deg + ... + p[deg]`` of degree `deg`
454 to points `(x, y)`. Returns a vector of coefficients `p` that minimises
455 the squared error in the order `deg`, `deg-1`, ... `0`.
457 The `Polynomial.fit <numpy.polynomial.polynomial.Polynomial.fit>` class
458 method is recommended for new code as it is more stable numerically. See
459 the documentation of the method for more information.
461 Parameters
462 ----------
463 x : array_like, shape (M,)
464 x-coordinates of the M sample points ``(x[i], y[i])``.
465 y : array_like, shape (M,) or (M, K)
466 y-coordinates of the sample points. Several data sets of sample
467 points sharing the same x-coordinates can be fitted at once by
468 passing in a 2D-array that contains one dataset per column.
469 deg : int
470 Degree of the fitting polynomial
471 rcond : float, optional
472 Relative condition number of the fit. Singular values smaller than
473 this relative to the largest singular value will be ignored. The
474 default value is len(x)*eps, where eps is the relative precision of
475 the float type, about 2e-16 in most cases.
476 full : bool, optional
477 Switch determining nature of return value. When it is False (the
478 default) just the coefficients are returned, when True diagnostic
479 information from the singular value decomposition is also returned.
480 w : array_like, shape (M,), optional
481 Weights. If not None, the weight ``w[i]`` applies to the unsquared
482 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
483 chosen so that the errors of the products ``w[i]*y[i]`` all have the
484 same variance. When using inverse-variance weighting, use
485 ``w[i] = 1/sigma(y[i])``. The default value is None.
486 cov : bool or str, optional
487 If given and not `False`, return not just the estimate but also its
488 covariance matrix. By default, the covariance are scaled by
489 chi2/dof, where dof = M - (deg + 1), i.e., the weights are presumed
490 to be unreliable except in a relative sense and everything is scaled
491 such that the reduced chi2 is unity. This scaling is omitted if
492 ``cov='unscaled'``, as is relevant for the case that the weights are
493 w = 1/sigma, with sigma known to be a reliable estimate of the
494 uncertainty.
496 Returns
497 -------
498 p : ndarray, shape (deg + 1,) or (deg + 1, K)
499 Polynomial coefficients, highest power first. If `y` was 2-D, the
500 coefficients for `k`-th data set are in ``p[:,k]``.
502 residuals, rank, singular_values, rcond
503 These values are only returned if ``full == True``
505 - residuals -- sum of squared residuals of the least squares fit
506 - rank -- the effective rank of the scaled Vandermonde
507 coefficient matrix
508 - singular_values -- singular values of the scaled Vandermonde
509 coefficient matrix
510 - rcond -- value of `rcond`.
512 For more details, see `numpy.linalg.lstsq`.
514 V : ndarray, shape (deg + 1, deg + 1) or (deg + 1, deg + 1, K)
515 Present only if ``full == False`` and ``cov == True``. The covariance
516 matrix of the polynomial coefficient estimates. The diagonal of
517 this matrix are the variance estimates for each coefficient. If y
518 is a 2-D array, then the covariance matrix for the `k`-th data set
519 are in ``V[:,:,k]``
522 Warns
523 -----
524 RankWarning
525 The rank of the coefficient matrix in the least-squares fit is
526 deficient. The warning is only raised if ``full == False``.
528 The warnings can be turned off by
530 >>> import warnings
531 >>> warnings.simplefilter('ignore', np.exceptions.RankWarning)
533 See Also
534 --------
535 polyval : Compute polynomial values.
536 linalg.lstsq : Computes a least-squares fit.
537 scipy.interpolate.UnivariateSpline : Computes spline fits.
539 Notes
540 -----
541 The solution minimizes the squared error
543 .. math::
544 E = \\sum_{j=0}^k |p(x_j) - y_j|^2
546 in the equations::
548 x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0]
549 x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1]
550 ...
551 x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k]
553 The coefficient matrix of the coefficients `p` is a Vandermonde matrix.
555 `polyfit` issues a `~exceptions.RankWarning` when the least-squares fit is
556 badly conditioned. This implies that the best fit is not well-defined due
557 to numerical error. The results may be improved by lowering the polynomial
558 degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter
559 can also be set to a value smaller than its default, but the resulting
560 fit may be spurious: including contributions from the small singular
561 values can add numerical noise to the result.
563 Note that fitting polynomial coefficients is inherently badly conditioned
564 when the degree of the polynomial is large or the interval of sample points
565 is badly centered. The quality of the fit should always be checked in these
566 cases. When polynomial fits are not satisfactory, splines may be a good
567 alternative.
569 References
570 ----------
571 .. [1] Wikipedia, "Curve fitting",
572 https://en.wikipedia.org/wiki/Curve_fitting
573 .. [2] Wikipedia, "Polynomial interpolation",
574 https://en.wikipedia.org/wiki/Polynomial_interpolation
576 Examples
577 --------
578 >>> import warnings
579 >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
580 >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
581 >>> z = np.polyfit(x, y, 3)
582 >>> z
583 array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) # may vary
585 It is convenient to use `poly1d` objects for dealing with polynomials:
587 >>> p = np.poly1d(z)
588 >>> p(0.5)
589 0.6143849206349179 # may vary
590 >>> p(3.5)
591 -0.34732142857143039 # may vary
592 >>> p(10)
593 22.579365079365115 # may vary
595 High-order polynomials may oscillate wildly:
597 >>> with warnings.catch_warnings():
598 ... warnings.simplefilter('ignore', np.exceptions.RankWarning)
599 ... p30 = np.poly1d(np.polyfit(x, y, 30))
600 ...
601 >>> p30(4)
602 -0.80000000000000204 # may vary
603 >>> p30(5)
604 -0.99999999999999445 # may vary
605 >>> p30(4.5)
606 -0.10547061179440398 # may vary
608 Illustration:
610 >>> import matplotlib.pyplot as plt
611 >>> xp = np.linspace(-2, 6, 100)
612 >>> _ = plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--')
613 >>> plt.ylim(-2,2)
614 (-2, 2)
615 >>> plt.show()
617 """
618 order = int(deg) + 1
619 x = NX.asarray(x) + 0.0
620 y = NX.asarray(y) + 0.0
622 # check arguments.
623 if deg < 0:
624 raise ValueError("expected deg >= 0")
625 if x.ndim != 1:
626 raise TypeError("expected 1D vector for x")
627 if x.size == 0:
628 raise TypeError("expected non-empty vector for x")
629 if y.ndim < 1 or y.ndim > 2:
630 raise TypeError("expected 1D or 2D array for y")
631 if x.shape[0] != y.shape[0]:
632 raise TypeError("expected x and y to have same length")
634 # set rcond
635 if rcond is None:
636 rcond = len(x)*finfo(x.dtype).eps
638 # set up least squares equation for powers of x
639 lhs = vander(x, order)
640 rhs = y
642 # apply weighting
643 if w is not None:
644 w = NX.asarray(w) + 0.0
645 if w.ndim != 1:
646 raise TypeError("expected a 1-d array for weights")
647 if w.shape[0] != y.shape[0]:
648 raise TypeError("expected w and y to have the same length")
649 lhs *= w[:, NX.newaxis]
650 if rhs.ndim == 2:
651 rhs *= w[:, NX.newaxis]
652 else:
653 rhs *= w
655 # scale lhs to improve condition number and solve
656 scale = NX.sqrt((lhs*lhs).sum(axis=0))
657 lhs /= scale
658 c, resids, rank, s = lstsq(lhs, rhs, rcond)
659 c = (c.T/scale).T # broadcast scale coefficients
661 # warn on rank reduction, which indicates an ill conditioned matrix
662 if rank != order and not full:
663 msg = "Polyfit may be poorly conditioned"
664 warnings.warn(msg, RankWarning, stacklevel=2)
666 if full:
667 return c, resids, rank, s, rcond
668 elif cov:
669 Vbase = inv(dot(lhs.T, lhs))
670 Vbase /= NX.outer(scale, scale)
671 if cov == "unscaled":
672 fac = 1
673 else:
674 if len(x) <= order:
675 raise ValueError("the number of data points must exceed order "
676 "to scale the covariance matrix")
677 # note, this used to be: fac = resids / (len(x) - order - 2.0)
678 # it was deciced that the "- 2" (originally justified by "Bayesian
679 # uncertainty analysis") is not what the user expects
680 # (see gh-11196 and gh-11197)
681 fac = resids / (len(x) - order)
682 if y.ndim == 1:
683 return c, Vbase * fac
684 else:
685 return c, Vbase[:,:, NX.newaxis] * fac
686 else:
687 return c
690def _polyval_dispatcher(p, x):
691 return (p, x)
694@array_function_dispatch(_polyval_dispatcher)
695def polyval(p, x):
696 """
697 Evaluate a polynomial at specific values.
699 .. note::
700 This forms part of the old polynomial API. Since version 1.4, the
701 new polynomial API defined in `numpy.polynomial` is preferred.
702 A summary of the differences can be found in the
703 :doc:`transition guide </reference/routines.polynomials>`.
705 If `p` is of length N, this function returns the value::
707 p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]
709 If `x` is a sequence, then ``p(x)`` is returned for each element of ``x``.
710 If `x` is another polynomial then the composite polynomial ``p(x(t))``
711 is returned.
713 Parameters
714 ----------
715 p : array_like or poly1d object
716 1D array of polynomial coefficients (including coefficients equal
717 to zero) from highest degree to the constant term, or an
718 instance of poly1d.
719 x : array_like or poly1d object
720 A number, an array of numbers, or an instance of poly1d, at
721 which to evaluate `p`.
723 Returns
724 -------
725 values : ndarray or poly1d
726 If `x` is a poly1d instance, the result is the composition of the two
727 polynomials, i.e., `x` is "substituted" in `p` and the simplified
728 result is returned. In addition, the type of `x` - array_like or
729 poly1d - governs the type of the output: `x` array_like => `values`
730 array_like, `x` a poly1d object => `values` is also.
732 See Also
733 --------
734 poly1d: A polynomial class.
736 Notes
737 -----
738 Horner's scheme [1]_ is used to evaluate the polynomial. Even so,
739 for polynomials of high degree the values may be inaccurate due to
740 rounding errors. Use carefully.
742 If `x` is a subtype of `ndarray` the return value will be of the same type.
744 References
745 ----------
746 .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng.
747 trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand
748 Reinhold Co., 1985, pg. 720.
750 Examples
751 --------
752 >>> np.polyval([3,0,1], 5) # 3 * 5**2 + 0 * 5**1 + 1
753 76
754 >>> np.polyval([3,0,1], np.poly1d(5))
755 poly1d([76])
756 >>> np.polyval(np.poly1d([3,0,1]), 5)
757 76
758 >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5))
759 poly1d([76])
761 """
762 p = NX.asarray(p)
763 if isinstance(x, poly1d):
764 y = 0
765 else:
766 x = NX.asanyarray(x)
767 y = NX.zeros_like(x)
768 for pv in p:
769 y = y * x + pv
770 return y
773def _binary_op_dispatcher(a1, a2):
774 return (a1, a2)
777@array_function_dispatch(_binary_op_dispatcher)
778def polyadd(a1, a2):
779 """
780 Find the sum of two polynomials.
782 .. note::
783 This forms part of the old polynomial API. Since version 1.4, the
784 new polynomial API defined in `numpy.polynomial` is preferred.
785 A summary of the differences can be found in the
786 :doc:`transition guide </reference/routines.polynomials>`.
788 Returns the polynomial resulting from the sum of two input polynomials.
789 Each input must be either a poly1d object or a 1D sequence of polynomial
790 coefficients, from highest to lowest degree.
792 Parameters
793 ----------
794 a1, a2 : array_like or poly1d object
795 Input polynomials.
797 Returns
798 -------
799 out : ndarray or poly1d object
800 The sum of the inputs. If either input is a poly1d object, then the
801 output is also a poly1d object. Otherwise, it is a 1D array of
802 polynomial coefficients from highest to lowest degree.
804 See Also
805 --------
806 poly1d : A one-dimensional polynomial class.
807 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
809 Examples
810 --------
811 >>> np.polyadd([1, 2], [9, 5, 4])
812 array([9, 6, 6])
814 Using poly1d objects:
816 >>> p1 = np.poly1d([1, 2])
817 >>> p2 = np.poly1d([9, 5, 4])
818 >>> print(p1)
819 1 x + 2
820 >>> print(p2)
821 2
822 9 x + 5 x + 4
823 >>> print(np.polyadd(p1, p2))
824 2
825 9 x + 6 x + 6
827 """
828 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
829 a1 = atleast_1d(a1)
830 a2 = atleast_1d(a2)
831 diff = len(a2) - len(a1)
832 if diff == 0:
833 val = a1 + a2
834 elif diff > 0:
835 zr = NX.zeros(diff, a1.dtype)
836 val = NX.concatenate((zr, a1)) + a2
837 else:
838 zr = NX.zeros(abs(diff), a2.dtype)
839 val = a1 + NX.concatenate((zr, a2))
840 if truepoly:
841 val = poly1d(val)
842 return val
845@array_function_dispatch(_binary_op_dispatcher)
846def polysub(a1, a2):
847 """
848 Difference (subtraction) of two polynomials.
850 .. note::
851 This forms part of the old polynomial API. Since version 1.4, the
852 new polynomial API defined in `numpy.polynomial` is preferred.
853 A summary of the differences can be found in the
854 :doc:`transition guide </reference/routines.polynomials>`.
856 Given two polynomials `a1` and `a2`, returns ``a1 - a2``.
857 `a1` and `a2` can be either array_like sequences of the polynomials'
858 coefficients (including coefficients equal to zero), or `poly1d` objects.
860 Parameters
861 ----------
862 a1, a2 : array_like or poly1d
863 Minuend and subtrahend polynomials, respectively.
865 Returns
866 -------
867 out : ndarray or poly1d
868 Array or `poly1d` object of the difference polynomial's coefficients.
870 See Also
871 --------
872 polyval, polydiv, polymul, polyadd
874 Examples
875 --------
876 .. math:: (2 x^2 + 10 x - 2) - (3 x^2 + 10 x -4) = (-x^2 + 2)
878 >>> np.polysub([2, 10, -2], [3, 10, -4])
879 array([-1, 0, 2])
881 """
882 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
883 a1 = atleast_1d(a1)
884 a2 = atleast_1d(a2)
885 diff = len(a2) - len(a1)
886 if diff == 0:
887 val = a1 - a2
888 elif diff > 0:
889 zr = NX.zeros(diff, a1.dtype)
890 val = NX.concatenate((zr, a1)) - a2
891 else:
892 zr = NX.zeros(abs(diff), a2.dtype)
893 val = a1 - NX.concatenate((zr, a2))
894 if truepoly:
895 val = poly1d(val)
896 return val
899@array_function_dispatch(_binary_op_dispatcher)
900def polymul(a1, a2):
901 """
902 Find the product of two polynomials.
904 .. note::
905 This forms part of the old polynomial API. Since version 1.4, the
906 new polynomial API defined in `numpy.polynomial` is preferred.
907 A summary of the differences can be found in the
908 :doc:`transition guide </reference/routines.polynomials>`.
910 Finds the polynomial resulting from the multiplication of the two input
911 polynomials. Each input must be either a poly1d object or a 1D sequence
912 of polynomial coefficients, from highest to lowest degree.
914 Parameters
915 ----------
916 a1, a2 : array_like or poly1d object
917 Input polynomials.
919 Returns
920 -------
921 out : ndarray or poly1d object
922 The polynomial resulting from the multiplication of the inputs. If
923 either inputs is a poly1d object, then the output is also a poly1d
924 object. Otherwise, it is a 1D array of polynomial coefficients from
925 highest to lowest degree.
927 See Also
928 --------
929 poly1d : A one-dimensional polynomial class.
930 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
931 convolve : Array convolution. Same output as polymul, but has parameter
932 for overlap mode.
934 Examples
935 --------
936 >>> np.polymul([1, 2, 3], [9, 5, 1])
937 array([ 9, 23, 38, 17, 3])
939 Using poly1d objects:
941 >>> p1 = np.poly1d([1, 2, 3])
942 >>> p2 = np.poly1d([9, 5, 1])
943 >>> print(p1)
944 2
945 1 x + 2 x + 3
946 >>> print(p2)
947 2
948 9 x + 5 x + 1
949 >>> print(np.polymul(p1, p2))
950 4 3 2
951 9 x + 23 x + 38 x + 17 x + 3
953 """
954 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
955 a1, a2 = poly1d(a1), poly1d(a2)
956 val = NX.convolve(a1, a2)
957 if truepoly:
958 val = poly1d(val)
959 return val
962def _polydiv_dispatcher(u, v):
963 return (u, v)
966@array_function_dispatch(_polydiv_dispatcher)
967def polydiv(u, v):
968 """
969 Returns the quotient and remainder of polynomial division.
971 .. note::
972 This forms part of the old polynomial API. Since version 1.4, the
973 new polynomial API defined in `numpy.polynomial` is preferred.
974 A summary of the differences can be found in the
975 :doc:`transition guide </reference/routines.polynomials>`.
977 The input arrays are the coefficients (including any coefficients
978 equal to zero) of the "numerator" (dividend) and "denominator"
979 (divisor) polynomials, respectively.
981 Parameters
982 ----------
983 u : array_like or poly1d
984 Dividend polynomial's coefficients.
986 v : array_like or poly1d
987 Divisor polynomial's coefficients.
989 Returns
990 -------
991 q : ndarray
992 Coefficients, including those equal to zero, of the quotient.
993 r : ndarray
994 Coefficients, including those equal to zero, of the remainder.
996 See Also
997 --------
998 poly, polyadd, polyder, polydiv, polyfit, polyint, polymul, polysub
999 polyval
1001 Notes
1002 -----
1003 Both `u` and `v` must be 0-d or 1-d (ndim = 0 or 1), but `u.ndim` need
1004 not equal `v.ndim`. In other words, all four possible combinations -
1005 ``u.ndim = v.ndim = 0``, ``u.ndim = v.ndim = 1``,
1006 ``u.ndim = 1, v.ndim = 0``, and ``u.ndim = 0, v.ndim = 1`` - work.
1008 Examples
1009 --------
1010 .. math:: \\frac{3x^2 + 5x + 2}{2x + 1} = 1.5x + 1.75, remainder 0.25
1012 >>> x = np.array([3.0, 5.0, 2.0])
1013 >>> y = np.array([2.0, 1.0])
1014 >>> np.polydiv(x, y)
1015 (array([1.5 , 1.75]), array([0.25]))
1017 """
1018 truepoly = (isinstance(u, poly1d) or isinstance(v, poly1d))
1019 u = atleast_1d(u) + 0.0
1020 v = atleast_1d(v) + 0.0
1021 # w has the common type
1022 w = u[0] + v[0]
1023 m = len(u) - 1
1024 n = len(v) - 1
1025 scale = 1. / v[0]
1026 q = NX.zeros((max(m - n + 1, 1),), w.dtype)
1027 r = u.astype(w.dtype)
1028 for k in range(0, m-n+1):
1029 d = scale * r[k]
1030 q[k] = d
1031 r[k:k+n+1] -= d*v
1032 while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
1033 r = r[1:]
1034 if truepoly:
1035 return poly1d(q), poly1d(r)
1036 return q, r
1038_poly_mat = re.compile(r"\*\*([0-9]*)")
1039def _raise_power(astr, wrap=70):
1040 n = 0
1041 line1 = ''
1042 line2 = ''
1043 output = ' '
1044 while True:
1045 mat = _poly_mat.search(astr, n)
1046 if mat is None:
1047 break
1048 span = mat.span()
1049 power = mat.groups()[0]
1050 partstr = astr[n:span[0]]
1051 n = span[1]
1052 toadd2 = partstr + ' '*(len(power)-1)
1053 toadd1 = ' '*(len(partstr)-1) + power
1054 if ((len(line2) + len(toadd2) > wrap) or
1055 (len(line1) + len(toadd1) > wrap)):
1056 output += line1 + "\n" + line2 + "\n "
1057 line1 = toadd1
1058 line2 = toadd2
1059 else:
1060 line2 += partstr + ' '*(len(power)-1)
1061 line1 += ' '*(len(partstr)-1) + power
1062 output += line1 + "\n" + line2
1063 return output + astr[n:]
1066@set_module('numpy')
1067class poly1d:
1068 """
1069 A one-dimensional polynomial class.
1071 .. note::
1072 This forms part of the old polynomial API. Since version 1.4, the
1073 new polynomial API defined in `numpy.polynomial` is preferred.
1074 A summary of the differences can be found in the
1075 :doc:`transition guide </reference/routines.polynomials>`.
1077 A convenience class, used to encapsulate "natural" operations on
1078 polynomials so that said operations may take on their customary
1079 form in code (see Examples).
1081 Parameters
1082 ----------
1083 c_or_r : array_like
1084 The polynomial's coefficients, in decreasing powers, or if
1085 the value of the second parameter is True, the polynomial's
1086 roots (values where the polynomial evaluates to 0). For example,
1087 ``poly1d([1, 2, 3])`` returns an object that represents
1088 :math:`x^2 + 2x + 3`, whereas ``poly1d([1, 2, 3], True)`` returns
1089 one that represents :math:`(x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x -6`.
1090 r : bool, optional
1091 If True, `c_or_r` specifies the polynomial's roots; the default
1092 is False.
1093 variable : str, optional
1094 Changes the variable used when printing `p` from `x` to `variable`
1095 (see Examples).
1097 Examples
1098 --------
1099 Construct the polynomial :math:`x^2 + 2x + 3`:
1101 >>> p = np.poly1d([1, 2, 3])
1102 >>> print(np.poly1d(p))
1103 2
1104 1 x + 2 x + 3
1106 Evaluate the polynomial at :math:`x = 0.5`:
1108 >>> p(0.5)
1109 4.25
1111 Find the roots:
1113 >>> p.r
1114 array([-1.+1.41421356j, -1.-1.41421356j])
1115 >>> p(p.r)
1116 array([ -4.44089210e-16+0.j, -4.44089210e-16+0.j]) # may vary
1118 These numbers in the previous line represent (0, 0) to machine precision
1120 Show the coefficients:
1122 >>> p.c
1123 array([1, 2, 3])
1125 Display the order (the leading zero-coefficients are removed):
1127 >>> p.order
1128 2
1130 Show the coefficient of the k-th power in the polynomial
1131 (which is equivalent to ``p.c[-(i+1)]``):
1133 >>> p[1]
1134 2
1136 Polynomials can be added, subtracted, multiplied, and divided
1137 (returns quotient and remainder):
1139 >>> p * p
1140 poly1d([ 1, 4, 10, 12, 9])
1142 >>> (p**3 + 4) / p
1143 (poly1d([ 1., 4., 10., 12., 9.]), poly1d([4.]))
1145 ``asarray(p)`` gives the coefficient array, so polynomials can be
1146 used in all functions that accept arrays:
1148 >>> p**2 # square of polynomial
1149 poly1d([ 1, 4, 10, 12, 9])
1151 >>> np.square(p) # square of individual coefficients
1152 array([1, 4, 9])
1154 The variable used in the string representation of `p` can be modified,
1155 using the `variable` parameter:
1157 >>> p = np.poly1d([1,2,3], variable='z')
1158 >>> print(p)
1159 2
1160 1 z + 2 z + 3
1162 Construct a polynomial from its roots:
1164 >>> np.poly1d([1, 2], True)
1165 poly1d([ 1., -3., 2.])
1167 This is the same polynomial as obtained by:
1169 >>> np.poly1d([1, -1]) * np.poly1d([1, -2])
1170 poly1d([ 1, -3, 2])
1172 """
1173 __hash__ = None
1175 @property
1176 def coeffs(self):
1177 """ The polynomial coefficients """
1178 return self._coeffs
1180 @coeffs.setter
1181 def coeffs(self, value):
1182 # allowing this makes p.coeffs *= 2 legal
1183 if value is not self._coeffs:
1184 raise AttributeError("Cannot set attribute")
1186 @property
1187 def variable(self):
1188 """ The name of the polynomial variable """
1189 return self._variable
1191 # calculated attributes
1192 @property
1193 def order(self):
1194 """ The order or degree of the polynomial """
1195 return len(self._coeffs) - 1
1197 @property
1198 def roots(self):
1199 """ The roots of the polynomial, where self(x) == 0 """
1200 return roots(self._coeffs)
1202 # our internal _coeffs property need to be backed by __dict__['coeffs'] for
1203 # scipy to work correctly.
1204 @property
1205 def _coeffs(self):
1206 return self.__dict__['coeffs']
1207 @_coeffs.setter
1208 def _coeffs(self, coeffs):
1209 self.__dict__['coeffs'] = coeffs
1211 # alias attributes
1212 r = roots
1213 c = coef = coefficients = coeffs
1214 o = order
1216 def __init__(self, c_or_r, r=False, variable=None):
1217 if isinstance(c_or_r, poly1d):
1218 self._variable = c_or_r._variable
1219 self._coeffs = c_or_r._coeffs
1221 if set(c_or_r.__dict__) - set(self.__dict__):
1222 msg = ("In the future extra properties will not be copied "
1223 "across when constructing one poly1d from another")
1224 warnings.warn(msg, FutureWarning, stacklevel=2)
1225 self.__dict__.update(c_or_r.__dict__)
1227 if variable is not None:
1228 self._variable = variable
1229 return
1230 if r:
1231 c_or_r = poly(c_or_r)
1232 c_or_r = atleast_1d(c_or_r)
1233 if c_or_r.ndim > 1:
1234 raise ValueError("Polynomial must be 1d only.")
1235 c_or_r = trim_zeros(c_or_r, trim='f')
1236 if len(c_or_r) == 0:
1237 c_or_r = NX.array([0], dtype=c_or_r.dtype)
1238 self._coeffs = c_or_r
1239 if variable is None:
1240 variable = 'x'
1241 self._variable = variable
1243 def __array__(self, t=None, copy=None):
1244 if t:
1245 return NX.asarray(self.coeffs, t, copy=copy)
1246 else:
1247 return NX.asarray(self.coeffs, copy=copy)
1249 def __repr__(self):
1250 vals = repr(self.coeffs)
1251 vals = vals[6:-1]
1252 return "poly1d(%s)" % vals
1254 def __len__(self):
1255 return self.order
1257 def __str__(self):
1258 thestr = "0"
1259 var = self.variable
1261 # Remove leading zeros
1262 coeffs = self.coeffs[NX.logical_or.accumulate(self.coeffs != 0)]
1263 N = len(coeffs)-1
1265 def fmt_float(q):
1266 s = '%.4g' % q
1267 if s.endswith('.0000'):
1268 s = s[:-5]
1269 return s
1271 for k, coeff in enumerate(coeffs):
1272 if not iscomplex(coeff):
1273 coefstr = fmt_float(real(coeff))
1274 elif real(coeff) == 0:
1275 coefstr = '%sj' % fmt_float(imag(coeff))
1276 else:
1277 coefstr = '(%s + %sj)' % (fmt_float(real(coeff)),
1278 fmt_float(imag(coeff)))
1280 power = (N-k)
1281 if power == 0:
1282 if coefstr != '0':
1283 newstr = '%s' % (coefstr,)
1284 else:
1285 if k == 0:
1286 newstr = '0'
1287 else:
1288 newstr = ''
1289 elif power == 1:
1290 if coefstr == '0':
1291 newstr = ''
1292 elif coefstr == 'b':
1293 newstr = var
1294 else:
1295 newstr = '%s %s' % (coefstr, var)
1296 else:
1297 if coefstr == '0':
1298 newstr = ''
1299 elif coefstr == 'b':
1300 newstr = '%s**%d' % (var, power,)
1301 else:
1302 newstr = '%s %s**%d' % (coefstr, var, power)
1304 if k > 0:
1305 if newstr != '':
1306 if newstr.startswith('-'):
1307 thestr = "%s - %s" % (thestr, newstr[1:])
1308 else:
1309 thestr = "%s + %s" % (thestr, newstr)
1310 else:
1311 thestr = newstr
1312 return _raise_power(thestr)
1314 def __call__(self, val):
1315 return polyval(self.coeffs, val)
1317 def __neg__(self):
1318 return poly1d(-self.coeffs)
1320 def __pos__(self):
1321 return self
1323 def __mul__(self, other):
1324 if isscalar(other):
1325 return poly1d(self.coeffs * other)
1326 else:
1327 other = poly1d(other)
1328 return poly1d(polymul(self.coeffs, other.coeffs))
1330 def __rmul__(self, other):
1331 if isscalar(other):
1332 return poly1d(other * self.coeffs)
1333 else:
1334 other = poly1d(other)
1335 return poly1d(polymul(self.coeffs, other.coeffs))
1337 def __add__(self, other):
1338 other = poly1d(other)
1339 return poly1d(polyadd(self.coeffs, other.coeffs))
1341 def __radd__(self, other):
1342 other = poly1d(other)
1343 return poly1d(polyadd(self.coeffs, other.coeffs))
1345 def __pow__(self, val):
1346 if not isscalar(val) or int(val) != val or val < 0:
1347 raise ValueError("Power to non-negative integers only.")
1348 res = [1]
1349 for _ in range(val):
1350 res = polymul(self.coeffs, res)
1351 return poly1d(res)
1353 def __sub__(self, other):
1354 other = poly1d(other)
1355 return poly1d(polysub(self.coeffs, other.coeffs))
1357 def __rsub__(self, other):
1358 other = poly1d(other)
1359 return poly1d(polysub(other.coeffs, self.coeffs))
1361 def __div__(self, other):
1362 if isscalar(other):
1363 return poly1d(self.coeffs/other)
1364 else:
1365 other = poly1d(other)
1366 return polydiv(self, other)
1368 __truediv__ = __div__
1370 def __rdiv__(self, other):
1371 if isscalar(other):
1372 return poly1d(other/self.coeffs)
1373 else:
1374 other = poly1d(other)
1375 return polydiv(other, self)
1377 __rtruediv__ = __rdiv__
1379 def __eq__(self, other):
1380 if not isinstance(other, poly1d):
1381 return NotImplemented
1382 if self.coeffs.shape != other.coeffs.shape:
1383 return False
1384 return (self.coeffs == other.coeffs).all()
1386 def __ne__(self, other):
1387 if not isinstance(other, poly1d):
1388 return NotImplemented
1389 return not self.__eq__(other)
1392 def __getitem__(self, val):
1393 ind = self.order - val
1394 if val > self.order:
1395 return self.coeffs.dtype.type(0)
1396 if val < 0:
1397 return self.coeffs.dtype.type(0)
1398 return self.coeffs[ind]
1400 def __setitem__(self, key, val):
1401 ind = self.order - key
1402 if key < 0:
1403 raise ValueError("Does not support negative powers.")
1404 if key > self.order:
1405 zr = NX.zeros(key-self.order, self.coeffs.dtype)
1406 self._coeffs = NX.concatenate((zr, self.coeffs))
1407 ind = 0
1408 self._coeffs[ind] = val
1409 return
1411 def __iter__(self):
1412 return iter(self.coeffs)
1414 def integ(self, m=1, k=0):
1415 """
1416 Return an antiderivative (indefinite integral) of this polynomial.
1418 Refer to `polyint` for full documentation.
1420 See Also
1421 --------
1422 polyint : equivalent function
1424 """
1425 return poly1d(polyint(self.coeffs, m=m, k=k))
1427 def deriv(self, m=1):
1428 """
1429 Return a derivative of this polynomial.
1431 Refer to `polyder` for full documentation.
1433 See Also
1434 --------
1435 polyder : equivalent function
1437 """
1438 return poly1d(polyder(self.coeffs, m=m))
1440# Stuff to do on module import
1442warnings.simplefilter('always', RankWarning)