Coverage for /var/srv/projects/api.amasfac.comuna18.com/tmp/venv/lib/python3.9/site-packages/numpy/lib/polynomial.py: 18%
397 statements
« prev ^ index » next coverage.py v6.4.4, created at 2023-07-17 14:22 -0600
« prev ^ index » next coverage.py v6.4.4, created at 2023-07-17 14:22 -0600
1"""
2Functions to operate on polynomials.
4"""
5__all__ = ['poly', 'roots', 'polyint', 'polyder', 'polyadd',
6 'polysub', 'polymul', 'polydiv', 'polyval', 'poly1d',
7 'polyfit', 'RankWarning']
9import functools
10import re
11import warnings
12import numpy.core.numeric as NX
14from numpy.core import (isscalar, abs, finfo, atleast_1d, hstack, dot, array,
15 ones)
16from numpy.core import overrides
17from numpy.core.overrides import set_module
18from numpy.lib.twodim_base import diag, vander
19from numpy.lib.function_base import trim_zeros
20from numpy.lib.type_check import iscomplex, real, imag, mintypecode
21from numpy.linalg import eigvals, lstsq, inv
24array_function_dispatch = functools.partial(
25 overrides.array_function_dispatch, module='numpy')
28@set_module('numpy')
29class RankWarning(UserWarning):
30 """
31 Issued by `polyfit` when the Vandermonde matrix is rank deficient.
33 For more information, a way to suppress the warning, and an example of
34 `RankWarning` being issued, see `polyfit`.
36 """
37 pass
40def _poly_dispatcher(seq_of_zeros):
41 return seq_of_zeros
44@array_function_dispatch(_poly_dispatcher)
45def poly(seq_of_zeros):
46 """
47 Find the coefficients of a polynomial with the given sequence of roots.
49 .. note::
50 This forms part of the old polynomial API. Since version 1.4, the
51 new polynomial API defined in `numpy.polynomial` is preferred.
52 A summary of the differences can be found in the
53 :doc:`transition guide </reference/routines.polynomials>`.
55 Returns the coefficients of the polynomial whose leading coefficient
56 is one for the given sequence of zeros (multiple roots must be included
57 in the sequence as many times as their multiplicity; see Examples).
58 A square matrix (or array, which will be treated as a matrix) can also
59 be given, in which case the coefficients of the characteristic polynomial
60 of the matrix are returned.
62 Parameters
63 ----------
64 seq_of_zeros : array_like, shape (N,) or (N, N)
65 A sequence of polynomial roots, or a square array or matrix object.
67 Returns
68 -------
69 c : ndarray
70 1D array of polynomial coefficients from highest to lowest degree:
72 ``c[0] * x**(N) + c[1] * x**(N-1) + ... + c[N-1] * x + c[N]``
73 where c[0] always equals 1.
75 Raises
76 ------
77 ValueError
78 If input is the wrong shape (the input must be a 1-D or square
79 2-D array).
81 See Also
82 --------
83 polyval : Compute polynomial values.
84 roots : Return the roots of a polynomial.
85 polyfit : Least squares polynomial fit.
86 poly1d : A one-dimensional polynomial class.
88 Notes
89 -----
90 Specifying the roots of a polynomial still leaves one degree of
91 freedom, typically represented by an undetermined leading
92 coefficient. [1]_ In the case of this function, that coefficient -
93 the first one in the returned array - is always taken as one. (If
94 for some reason you have one other point, the only automatic way
95 presently to leverage that information is to use ``polyfit``.)
97 The characteristic polynomial, :math:`p_a(t)`, of an `n`-by-`n`
98 matrix **A** is given by
100 :math:`p_a(t) = \\mathrm{det}(t\\, \\mathbf{I} - \\mathbf{A})`,
102 where **I** is the `n`-by-`n` identity matrix. [2]_
104 References
105 ----------
106 .. [1] M. Sullivan and M. Sullivan, III, "Algebra and Trignometry,
107 Enhanced With Graphing Utilities," Prentice-Hall, pg. 318, 1996.
109 .. [2] G. Strang, "Linear Algebra and Its Applications, 2nd Edition,"
110 Academic Press, pg. 182, 1980.
112 Examples
113 --------
114 Given a sequence of a polynomial's zeros:
116 >>> np.poly((0, 0, 0)) # Multiple root example
117 array([1., 0., 0., 0.])
119 The line above represents z**3 + 0*z**2 + 0*z + 0.
121 >>> np.poly((-1./2, 0, 1./2))
122 array([ 1. , 0. , -0.25, 0. ])
124 The line above represents z**3 - z/4
126 >>> np.poly((np.random.random(1)[0], 0, np.random.random(1)[0]))
127 array([ 1. , -0.77086955, 0.08618131, 0. ]) # random
129 Given a square array object:
131 >>> P = np.array([[0, 1./3], [-1./2, 0]])
132 >>> np.poly(P)
133 array([1. , 0. , 0.16666667])
135 Note how in all cases the leading coefficient is always 1.
137 """
138 seq_of_zeros = atleast_1d(seq_of_zeros)
139 sh = seq_of_zeros.shape
141 if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0:
142 seq_of_zeros = eigvals(seq_of_zeros)
143 elif len(sh) == 1:
144 dt = seq_of_zeros.dtype
145 # Let object arrays slip through, e.g. for arbitrary precision
146 if dt != object:
147 seq_of_zeros = seq_of_zeros.astype(mintypecode(dt.char))
148 else:
149 raise ValueError("input must be 1d or non-empty square 2d array.")
151 if len(seq_of_zeros) == 0:
152 return 1.0
153 dt = seq_of_zeros.dtype
154 a = ones((1,), dtype=dt)
155 for zero in seq_of_zeros:
156 a = NX.convolve(a, array([1, -zero], dtype=dt), mode='full')
158 if issubclass(a.dtype.type, NX.complexfloating):
159 # if complex roots are all complex conjugates, the roots are real.
160 roots = NX.asarray(seq_of_zeros, complex)
161 if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())):
162 a = a.real.copy()
164 return a
167def _roots_dispatcher(p):
168 return p
171@array_function_dispatch(_roots_dispatcher)
172def roots(p):
173 """
174 Return the roots of a polynomial with coefficients given in p.
176 .. note::
177 This forms part of the old polynomial API. Since version 1.4, the
178 new polynomial API defined in `numpy.polynomial` is preferred.
179 A summary of the differences can be found in the
180 :doc:`transition guide </reference/routines.polynomials>`.
182 The values in the rank-1 array `p` are coefficients of a polynomial.
183 If the length of `p` is n+1 then the polynomial is described by::
185 p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
187 Parameters
188 ----------
189 p : array_like
190 Rank-1 array of polynomial coefficients.
192 Returns
193 -------
194 out : ndarray
195 An array containing the roots of the polynomial.
197 Raises
198 ------
199 ValueError
200 When `p` cannot be converted to a rank-1 array.
202 See also
203 --------
204 poly : Find the coefficients of a polynomial with a given sequence
205 of roots.
206 polyval : Compute polynomial values.
207 polyfit : Least squares polynomial fit.
208 poly1d : A one-dimensional polynomial class.
210 Notes
211 -----
212 The algorithm relies on computing the eigenvalues of the
213 companion matrix [1]_.
215 References
216 ----------
217 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK:
218 Cambridge University Press, 1999, pp. 146-7.
220 Examples
221 --------
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 --------
308 The defining property of the antiderivative:
310 >>> p = np.poly1d([1,1,1])
311 >>> P = np.polyint(p)
312 >>> P
313 poly1d([ 0.33333333, 0.5 , 1. , 0. ]) # may vary
314 >>> np.polyder(P) == p
315 True
317 The integration constants default to zero, but can be specified:
319 >>> P = np.polyint(p, 3)
320 >>> P(0)
321 0.0
322 >>> np.polyder(P)(0)
323 0.0
324 >>> np.polyder(P, 2)(0)
325 0.0
326 >>> P = np.polyint(p, 3, k=[6,5,3])
327 >>> P
328 poly1d([ 0.01666667, 0.04166667, 0.16666667, 3. , 5. , 3. ]) # may vary
330 Note that 3 = 6 / 2!, and that the constants are given in the order of
331 integrations. Constant of the highest-order polynomial term comes first:
333 >>> np.polyder(P, 2)(0)
334 6.0
335 >>> np.polyder(P, 1)(0)
336 5.0
337 >>> P(0)
338 3.0
340 """
341 m = int(m)
342 if m < 0:
343 raise ValueError("Order of integral must be positive (see polyder)")
344 if k is None:
345 k = NX.zeros(m, float)
346 k = atleast_1d(k)
347 if len(k) == 1 and m > 1:
348 k = k[0]*NX.ones(m, float)
349 if len(k) < m:
350 raise ValueError(
351 "k must be a scalar or a rank-1 array of length 1 or >m.")
353 truepoly = isinstance(p, poly1d)
354 p = NX.asarray(p)
355 if m == 0:
356 if truepoly:
357 return poly1d(p)
358 return p
359 else:
360 # Note: this must work also with object and integer arrays
361 y = NX.concatenate((p.__truediv__(NX.arange(len(p), 0, -1)), [k[0]]))
362 val = polyint(y, m - 1, k=k[1:])
363 if truepoly:
364 return poly1d(val)
365 return val
368def _polyder_dispatcher(p, m=None):
369 return (p,)
372@array_function_dispatch(_polyder_dispatcher)
373def polyder(p, m=1):
374 """
375 Return the derivative of the specified order of a polynomial.
377 .. note::
378 This forms part of the old polynomial API. Since version 1.4, the
379 new polynomial API defined in `numpy.polynomial` is preferred.
380 A summary of the differences can be found in the
381 :doc:`transition guide </reference/routines.polynomials>`.
383 Parameters
384 ----------
385 p : poly1d or sequence
386 Polynomial to differentiate.
387 A sequence is interpreted as polynomial coefficients, see `poly1d`.
388 m : int, optional
389 Order of differentiation (default: 1)
391 Returns
392 -------
393 der : poly1d
394 A new polynomial representing the derivative.
396 See Also
397 --------
398 polyint : Anti-derivative of a polynomial.
399 poly1d : Class for one-dimensional polynomials.
401 Examples
402 --------
403 The derivative of the polynomial :math:`x^3 + x^2 + x^1 + 1` is:
405 >>> p = np.poly1d([1,1,1,1])
406 >>> p2 = np.polyder(p)
407 >>> p2
408 poly1d([3, 2, 1])
410 which evaluates to:
412 >>> p2(2.)
413 17.0
415 We can verify this, approximating the derivative with
416 ``(f(x + h) - f(x))/h``:
418 >>> (p(2. + 0.001) - p(2.)) / 0.001
419 17.007000999997857
421 The fourth-order derivative of a 3rd-order polynomial is zero:
423 >>> np.polyder(p, 2)
424 poly1d([6, 2])
425 >>> np.polyder(p, 3)
426 poly1d([6])
427 >>> np.polyder(p, 4)
428 poly1d([0])
430 """
431 m = int(m)
432 if m < 0:
433 raise ValueError("Order of derivative must be positive (see polyint)")
435 truepoly = isinstance(p, poly1d)
436 p = NX.asarray(p)
437 n = len(p) - 1
438 y = p[:-1] * NX.arange(n, 0, -1)
439 if m == 0:
440 val = p
441 else:
442 val = polyder(y, m - 1)
443 if truepoly:
444 val = poly1d(val)
445 return val
448def _polyfit_dispatcher(x, y, deg, rcond=None, full=None, w=None, cov=None):
449 return (x, y, w)
452@array_function_dispatch(_polyfit_dispatcher)
453def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
454 """
455 Least squares polynomial fit.
457 .. note::
458 This forms part of the old polynomial API. Since version 1.4, the
459 new polynomial API defined in `numpy.polynomial` is preferred.
460 A summary of the differences can be found in the
461 :doc:`transition guide </reference/routines.polynomials>`.
463 Fit a polynomial ``p(x) = p[0] * x**deg + ... + p[deg]`` of degree `deg`
464 to points `(x, y)`. Returns a vector of coefficients `p` that minimises
465 the squared error in the order `deg`, `deg-1`, ... `0`.
467 The `Polynomial.fit <numpy.polynomial.polynomial.Polynomial.fit>` class
468 method is recommended for new code as it is more stable numerically. See
469 the documentation of the method for more information.
471 Parameters
472 ----------
473 x : array_like, shape (M,)
474 x-coordinates of the M sample points ``(x[i], y[i])``.
475 y : array_like, shape (M,) or (M, K)
476 y-coordinates of the sample points. Several data sets of sample
477 points sharing the same x-coordinates can be fitted at once by
478 passing in a 2D-array that contains one dataset per column.
479 deg : int
480 Degree of the fitting polynomial
481 rcond : float, optional
482 Relative condition number of the fit. Singular values smaller than
483 this relative to the largest singular value will be ignored. The
484 default value is len(x)*eps, where eps is the relative precision of
485 the float type, about 2e-16 in most cases.
486 full : bool, optional
487 Switch determining nature of return value. When it is False (the
488 default) just the coefficients are returned, when True diagnostic
489 information from the singular value decomposition is also returned.
490 w : array_like, shape (M,), optional
491 Weights. If not None, the weight ``w[i]`` applies to the unsquared
492 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
493 chosen so that the errors of the products ``w[i]*y[i]`` all have the
494 same variance. When using inverse-variance weighting, use
495 ``w[i] = 1/sigma(y[i])``. The default value is None.
496 cov : bool or str, optional
497 If given and not `False`, return not just the estimate but also its
498 covariance matrix. By default, the covariance are scaled by
499 chi2/dof, where dof = M - (deg + 1), i.e., the weights are presumed
500 to be unreliable except in a relative sense and everything is scaled
501 such that the reduced chi2 is unity. This scaling is omitted if
502 ``cov='unscaled'``, as is relevant for the case that the weights are
503 w = 1/sigma, with sigma known to be a reliable estimate of the
504 uncertainty.
506 Returns
507 -------
508 p : ndarray, shape (deg + 1,) or (deg + 1, K)
509 Polynomial coefficients, highest power first. If `y` was 2-D, the
510 coefficients for `k`-th data set are in ``p[:,k]``.
512 residuals, rank, singular_values, rcond
513 These values are only returned if ``full == True``
515 - residuals -- sum of squared residuals of the least squares fit
516 - rank -- the effective rank of the scaled Vandermonde
517 coefficient matrix
518 - singular_values -- singular values of the scaled Vandermonde
519 coefficient matrix
520 - rcond -- value of `rcond`.
522 For more details, see `numpy.linalg.lstsq`.
524 V : ndarray, shape (M,M) or (M,M,K)
525 Present only if ``full == False`` and ``cov == True``. The covariance
526 matrix of the polynomial coefficient estimates. The diagonal of
527 this matrix are the variance estimates for each coefficient. If y
528 is a 2-D array, then the covariance matrix for the `k`-th data set
529 are in ``V[:,:,k]``
532 Warns
533 -----
534 RankWarning
535 The rank of the coefficient matrix in the least-squares fit is
536 deficient. The warning is only raised if ``full == False``.
538 The warnings can be turned off by
540 >>> import warnings
541 >>> warnings.simplefilter('ignore', np.RankWarning)
543 See Also
544 --------
545 polyval : Compute polynomial values.
546 linalg.lstsq : Computes a least-squares fit.
547 scipy.interpolate.UnivariateSpline : Computes spline fits.
549 Notes
550 -----
551 The solution minimizes the squared error
553 .. math::
554 E = \\sum_{j=0}^k |p(x_j) - y_j|^2
556 in the equations::
558 x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0]
559 x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1]
560 ...
561 x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k]
563 The coefficient matrix of the coefficients `p` is a Vandermonde matrix.
565 `polyfit` issues a `RankWarning` when the least-squares fit is badly
566 conditioned. This implies that the best fit is not well-defined due
567 to numerical error. The results may be improved by lowering the polynomial
568 degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter
569 can also be set to a value smaller than its default, but the resulting
570 fit may be spurious: including contributions from the small singular
571 values can add numerical noise to the result.
573 Note that fitting polynomial coefficients is inherently badly conditioned
574 when the degree of the polynomial is large or the interval of sample points
575 is badly centered. The quality of the fit should always be checked in these
576 cases. When polynomial fits are not satisfactory, splines may be a good
577 alternative.
579 References
580 ----------
581 .. [1] Wikipedia, "Curve fitting",
582 https://en.wikipedia.org/wiki/Curve_fitting
583 .. [2] Wikipedia, "Polynomial interpolation",
584 https://en.wikipedia.org/wiki/Polynomial_interpolation
586 Examples
587 --------
588 >>> import warnings
589 >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
590 >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
591 >>> z = np.polyfit(x, y, 3)
592 >>> z
593 array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) # may vary
595 It is convenient to use `poly1d` objects for dealing with polynomials:
597 >>> p = np.poly1d(z)
598 >>> p(0.5)
599 0.6143849206349179 # may vary
600 >>> p(3.5)
601 -0.34732142857143039 # may vary
602 >>> p(10)
603 22.579365079365115 # may vary
605 High-order polynomials may oscillate wildly:
607 >>> with warnings.catch_warnings():
608 ... warnings.simplefilter('ignore', np.RankWarning)
609 ... p30 = np.poly1d(np.polyfit(x, y, 30))
610 ...
611 >>> p30(4)
612 -0.80000000000000204 # may vary
613 >>> p30(5)
614 -0.99999999999999445 # may vary
615 >>> p30(4.5)
616 -0.10547061179440398 # may vary
618 Illustration:
620 >>> import matplotlib.pyplot as plt
621 >>> xp = np.linspace(-2, 6, 100)
622 >>> _ = plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--')
623 >>> plt.ylim(-2,2)
624 (-2, 2)
625 >>> plt.show()
627 """
628 order = int(deg) + 1
629 x = NX.asarray(x) + 0.0
630 y = NX.asarray(y) + 0.0
632 # check arguments.
633 if deg < 0:
634 raise ValueError("expected deg >= 0")
635 if x.ndim != 1:
636 raise TypeError("expected 1D vector for x")
637 if x.size == 0:
638 raise TypeError("expected non-empty vector for x")
639 if y.ndim < 1 or y.ndim > 2:
640 raise TypeError("expected 1D or 2D array for y")
641 if x.shape[0] != y.shape[0]:
642 raise TypeError("expected x and y to have same length")
644 # set rcond
645 if rcond is None:
646 rcond = len(x)*finfo(x.dtype).eps
648 # set up least squares equation for powers of x
649 lhs = vander(x, order)
650 rhs = y
652 # apply weighting
653 if w is not None:
654 w = NX.asarray(w) + 0.0
655 if w.ndim != 1:
656 raise TypeError("expected a 1-d array for weights")
657 if w.shape[0] != y.shape[0]:
658 raise TypeError("expected w and y to have the same length")
659 lhs *= w[:, NX.newaxis]
660 if rhs.ndim == 2:
661 rhs *= w[:, NX.newaxis]
662 else:
663 rhs *= w
665 # scale lhs to improve condition number and solve
666 scale = NX.sqrt((lhs*lhs).sum(axis=0))
667 lhs /= scale
668 c, resids, rank, s = lstsq(lhs, rhs, rcond)
669 c = (c.T/scale).T # broadcast scale coefficients
671 # warn on rank reduction, which indicates an ill conditioned matrix
672 if rank != order and not full:
673 msg = "Polyfit may be poorly conditioned"
674 warnings.warn(msg, RankWarning, stacklevel=4)
676 if full:
677 return c, resids, rank, s, rcond
678 elif cov:
679 Vbase = inv(dot(lhs.T, lhs))
680 Vbase /= NX.outer(scale, scale)
681 if cov == "unscaled":
682 fac = 1
683 else:
684 if len(x) <= order:
685 raise ValueError("the number of data points must exceed order "
686 "to scale the covariance matrix")
687 # note, this used to be: fac = resids / (len(x) - order - 2.0)
688 # it was deciced that the "- 2" (originally justified by "Bayesian
689 # uncertainty analysis") is not what the user expects
690 # (see gh-11196 and gh-11197)
691 fac = resids / (len(x) - order)
692 if y.ndim == 1:
693 return c, Vbase * fac
694 else:
695 return c, Vbase[:,:, NX.newaxis] * fac
696 else:
697 return c
700def _polyval_dispatcher(p, x):
701 return (p, x)
704@array_function_dispatch(_polyval_dispatcher)
705def polyval(p, x):
706 """
707 Evaluate a polynomial at specific values.
709 .. note::
710 This forms part of the old polynomial API. Since version 1.4, the
711 new polynomial API defined in `numpy.polynomial` is preferred.
712 A summary of the differences can be found in the
713 :doc:`transition guide </reference/routines.polynomials>`.
715 If `p` is of length N, this function returns the value:
717 ``p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]``
719 If `x` is a sequence, then ``p(x)`` is returned for each element of ``x``.
720 If `x` is another polynomial then the composite polynomial ``p(x(t))``
721 is returned.
723 Parameters
724 ----------
725 p : array_like or poly1d object
726 1D array of polynomial coefficients (including coefficients equal
727 to zero) from highest degree to the constant term, or an
728 instance of poly1d.
729 x : array_like or poly1d object
730 A number, an array of numbers, or an instance of poly1d, at
731 which to evaluate `p`.
733 Returns
734 -------
735 values : ndarray or poly1d
736 If `x` is a poly1d instance, the result is the composition of the two
737 polynomials, i.e., `x` is "substituted" in `p` and the simplified
738 result is returned. In addition, the type of `x` - array_like or
739 poly1d - governs the type of the output: `x` array_like => `values`
740 array_like, `x` a poly1d object => `values` is also.
742 See Also
743 --------
744 poly1d: A polynomial class.
746 Notes
747 -----
748 Horner's scheme [1]_ is used to evaluate the polynomial. Even so,
749 for polynomials of high degree the values may be inaccurate due to
750 rounding errors. Use carefully.
752 If `x` is a subtype of `ndarray` the return value will be of the same type.
754 References
755 ----------
756 .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng.
757 trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand
758 Reinhold Co., 1985, pg. 720.
760 Examples
761 --------
762 >>> np.polyval([3,0,1], 5) # 3 * 5**2 + 0 * 5**1 + 1
763 76
764 >>> np.polyval([3,0,1], np.poly1d(5))
765 poly1d([76])
766 >>> np.polyval(np.poly1d([3,0,1]), 5)
767 76
768 >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5))
769 poly1d([76])
771 """
772 p = NX.asarray(p)
773 if isinstance(x, poly1d):
774 y = 0
775 else:
776 x = NX.asanyarray(x)
777 y = NX.zeros_like(x)
778 for pv in p:
779 y = y * x + pv
780 return y
783def _binary_op_dispatcher(a1, a2):
784 return (a1, a2)
787@array_function_dispatch(_binary_op_dispatcher)
788def polyadd(a1, a2):
789 """
790 Find the sum of two polynomials.
792 .. note::
793 This forms part of the old polynomial API. Since version 1.4, the
794 new polynomial API defined in `numpy.polynomial` is preferred.
795 A summary of the differences can be found in the
796 :doc:`transition guide </reference/routines.polynomials>`.
798 Returns the polynomial resulting from the sum of two input polynomials.
799 Each input must be either a poly1d object or a 1D sequence of polynomial
800 coefficients, from highest to lowest degree.
802 Parameters
803 ----------
804 a1, a2 : array_like or poly1d object
805 Input polynomials.
807 Returns
808 -------
809 out : ndarray or poly1d object
810 The sum of the inputs. If either input is a poly1d object, then the
811 output is also a poly1d object. Otherwise, it is a 1D array of
812 polynomial coefficients from highest to lowest degree.
814 See Also
815 --------
816 poly1d : A one-dimensional polynomial class.
817 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
819 Examples
820 --------
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 >>> np.polysub([2, 10, -2], [3, 10, -4])
889 array([-1, 0, 2])
891 """
892 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
893 a1 = atleast_1d(a1)
894 a2 = atleast_1d(a2)
895 diff = len(a2) - len(a1)
896 if diff == 0:
897 val = a1 - a2
898 elif diff > 0:
899 zr = NX.zeros(diff, a1.dtype)
900 val = NX.concatenate((zr, a1)) - a2
901 else:
902 zr = NX.zeros(abs(diff), a2.dtype)
903 val = a1 - NX.concatenate((zr, a2))
904 if truepoly:
905 val = poly1d(val)
906 return val
909@array_function_dispatch(_binary_op_dispatcher)
910def polymul(a1, a2):
911 """
912 Find the product of two polynomials.
914 .. note::
915 This forms part of the old polynomial API. Since version 1.4, the
916 new polynomial API defined in `numpy.polynomial` is preferred.
917 A summary of the differences can be found in the
918 :doc:`transition guide </reference/routines.polynomials>`.
920 Finds the polynomial resulting from the multiplication of the two input
921 polynomials. Each input must be either a poly1d object or a 1D sequence
922 of polynomial coefficients, from highest to lowest degree.
924 Parameters
925 ----------
926 a1, a2 : array_like or poly1d object
927 Input polynomials.
929 Returns
930 -------
931 out : ndarray or poly1d object
932 The polynomial resulting from the multiplication of the inputs. If
933 either inputs is a poly1d object, then the output is also a poly1d
934 object. Otherwise, it is a 1D array of polynomial coefficients from
935 highest to lowest degree.
937 See Also
938 --------
939 poly1d : A one-dimensional polynomial class.
940 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
941 convolve : Array convolution. Same output as polymul, but has parameter
942 for overlap mode.
944 Examples
945 --------
946 >>> np.polymul([1, 2, 3], [9, 5, 1])
947 array([ 9, 23, 38, 17, 3])
949 Using poly1d objects:
951 >>> p1 = np.poly1d([1, 2, 3])
952 >>> p2 = np.poly1d([9, 5, 1])
953 >>> print(p1)
954 2
955 1 x + 2 x + 3
956 >>> print(p2)
957 2
958 9 x + 5 x + 1
959 >>> print(np.polymul(p1, p2))
960 4 3 2
961 9 x + 23 x + 38 x + 17 x + 3
963 """
964 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
965 a1, a2 = poly1d(a1), poly1d(a2)
966 val = NX.convolve(a1, a2)
967 if truepoly:
968 val = poly1d(val)
969 return val
972def _polydiv_dispatcher(u, v):
973 return (u, v)
976@array_function_dispatch(_polydiv_dispatcher)
977def polydiv(u, v):
978 """
979 Returns the quotient and remainder of polynomial division.
981 .. note::
982 This forms part of the old polynomial API. Since version 1.4, the
983 new polynomial API defined in `numpy.polynomial` is preferred.
984 A summary of the differences can be found in the
985 :doc:`transition guide </reference/routines.polynomials>`.
987 The input arrays are the coefficients (including any coefficients
988 equal to zero) of the "numerator" (dividend) and "denominator"
989 (divisor) polynomials, respectively.
991 Parameters
992 ----------
993 u : array_like or poly1d
994 Dividend polynomial's coefficients.
996 v : array_like or poly1d
997 Divisor polynomial's coefficients.
999 Returns
1000 -------
1001 q : ndarray
1002 Coefficients, including those equal to zero, of the quotient.
1003 r : ndarray
1004 Coefficients, including those equal to zero, of the remainder.
1006 See Also
1007 --------
1008 poly, polyadd, polyder, polydiv, polyfit, polyint, polymul, polysub
1009 polyval
1011 Notes
1012 -----
1013 Both `u` and `v` must be 0-d or 1-d (ndim = 0 or 1), but `u.ndim` need
1014 not equal `v.ndim`. In other words, all four possible combinations -
1015 ``u.ndim = v.ndim = 0``, ``u.ndim = v.ndim = 1``,
1016 ``u.ndim = 1, v.ndim = 0``, and ``u.ndim = 0, v.ndim = 1`` - work.
1018 Examples
1019 --------
1020 .. math:: \\frac{3x^2 + 5x + 2}{2x + 1} = 1.5x + 1.75, remainder 0.25
1022 >>> x = np.array([3.0, 5.0, 2.0])
1023 >>> y = np.array([2.0, 1.0])
1024 >>> np.polydiv(x, y)
1025 (array([1.5 , 1.75]), array([0.25]))
1027 """
1028 truepoly = (isinstance(u, poly1d) or isinstance(v, poly1d))
1029 u = atleast_1d(u) + 0.0
1030 v = atleast_1d(v) + 0.0
1031 # w has the common type
1032 w = u[0] + v[0]
1033 m = len(u) - 1
1034 n = len(v) - 1
1035 scale = 1. / v[0]
1036 q = NX.zeros((max(m - n + 1, 1),), w.dtype)
1037 r = u.astype(w.dtype)
1038 for k in range(0, m-n+1):
1039 d = scale * r[k]
1040 q[k] = d
1041 r[k:k+n+1] -= d*v
1042 while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
1043 r = r[1:]
1044 if truepoly:
1045 return poly1d(q), poly1d(r)
1046 return q, r
1048_poly_mat = re.compile(r"\*\*([0-9]*)")
1049def _raise_power(astr, wrap=70):
1050 n = 0
1051 line1 = ''
1052 line2 = ''
1053 output = ' '
1054 while True:
1055 mat = _poly_mat.search(astr, n)
1056 if mat is None:
1057 break
1058 span = mat.span()
1059 power = mat.groups()[0]
1060 partstr = astr[n:span[0]]
1061 n = span[1]
1062 toadd2 = partstr + ' '*(len(power)-1)
1063 toadd1 = ' '*(len(partstr)-1) + power
1064 if ((len(line2) + len(toadd2) > wrap) or
1065 (len(line1) + len(toadd1) > wrap)):
1066 output += line1 + "\n" + line2 + "\n "
1067 line1 = toadd1
1068 line2 = toadd2
1069 else:
1070 line2 += partstr + ' '*(len(power)-1)
1071 line1 += ' '*(len(partstr)-1) + power
1072 output += line1 + "\n" + line2
1073 return output + astr[n:]
1076@set_module('numpy')
1077class poly1d:
1078 """
1079 A one-dimensional polynomial class.
1081 .. note::
1082 This forms part of the old polynomial API. Since version 1.4, the
1083 new polynomial API defined in `numpy.polynomial` is preferred.
1084 A summary of the differences can be found in the
1085 :doc:`transition guide </reference/routines.polynomials>`.
1087 A convenience class, used to encapsulate "natural" operations on
1088 polynomials so that said operations may take on their customary
1089 form in code (see Examples).
1091 Parameters
1092 ----------
1093 c_or_r : array_like
1094 The polynomial's coefficients, in decreasing powers, or if
1095 the value of the second parameter is True, the polynomial's
1096 roots (values where the polynomial evaluates to 0). For example,
1097 ``poly1d([1, 2, 3])`` returns an object that represents
1098 :math:`x^2 + 2x + 3`, whereas ``poly1d([1, 2, 3], True)`` returns
1099 one that represents :math:`(x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x -6`.
1100 r : bool, optional
1101 If True, `c_or_r` specifies the polynomial's roots; the default
1102 is False.
1103 variable : str, optional
1104 Changes the variable used when printing `p` from `x` to `variable`
1105 (see Examples).
1107 Examples
1108 --------
1109 Construct the polynomial :math:`x^2 + 2x + 3`:
1111 >>> p = np.poly1d([1, 2, 3])
1112 >>> print(np.poly1d(p))
1113 2
1114 1 x + 2 x + 3
1116 Evaluate the polynomial at :math:`x = 0.5`:
1118 >>> p(0.5)
1119 4.25
1121 Find the roots:
1123 >>> p.r
1124 array([-1.+1.41421356j, -1.-1.41421356j])
1125 >>> p(p.r)
1126 array([ -4.44089210e-16+0.j, -4.44089210e-16+0.j]) # may vary
1128 These numbers in the previous line represent (0, 0) to machine precision
1130 Show the coefficients:
1132 >>> p.c
1133 array([1, 2, 3])
1135 Display the order (the leading zero-coefficients are removed):
1137 >>> p.order
1138 2
1140 Show the coefficient of the k-th power in the polynomial
1141 (which is equivalent to ``p.c[-(i+1)]``):
1143 >>> p[1]
1144 2
1146 Polynomials can be added, subtracted, multiplied, and divided
1147 (returns quotient and remainder):
1149 >>> p * p
1150 poly1d([ 1, 4, 10, 12, 9])
1152 >>> (p**3 + 4) / p
1153 (poly1d([ 1., 4., 10., 12., 9.]), poly1d([4.]))
1155 ``asarray(p)`` gives the coefficient array, so polynomials can be
1156 used in all functions that accept arrays:
1158 >>> p**2 # square of polynomial
1159 poly1d([ 1, 4, 10, 12, 9])
1161 >>> np.square(p) # square of individual coefficients
1162 array([1, 4, 9])
1164 The variable used in the string representation of `p` can be modified,
1165 using the `variable` parameter:
1167 >>> p = np.poly1d([1,2,3], variable='z')
1168 >>> print(p)
1169 2
1170 1 z + 2 z + 3
1172 Construct a polynomial from its roots:
1174 >>> np.poly1d([1, 2], True)
1175 poly1d([ 1., -3., 2.])
1177 This is the same polynomial as obtained by:
1179 >>> np.poly1d([1, -1]) * np.poly1d([1, -2])
1180 poly1d([ 1, -3, 2])
1182 """
1183 __hash__ = None
1185 @property
1186 def coeffs(self):
1187 """ The polynomial coefficients """
1188 return self._coeffs
1190 @coeffs.setter
1191 def coeffs(self, value):
1192 # allowing this makes p.coeffs *= 2 legal
1193 if value is not self._coeffs:
1194 raise AttributeError("Cannot set attribute")
1196 @property
1197 def variable(self):
1198 """ The name of the polynomial variable """
1199 return self._variable
1201 # calculated attributes
1202 @property
1203 def order(self):
1204 """ The order or degree of the polynomial """
1205 return len(self._coeffs) - 1
1207 @property
1208 def roots(self):
1209 """ The roots of the polynomial, where self(x) == 0 """
1210 return roots(self._coeffs)
1212 # our internal _coeffs property need to be backed by __dict__['coeffs'] for
1213 # scipy to work correctly.
1214 @property
1215 def _coeffs(self):
1216 return self.__dict__['coeffs']
1217 @_coeffs.setter
1218 def _coeffs(self, coeffs):
1219 self.__dict__['coeffs'] = coeffs
1221 # alias attributes
1222 r = roots
1223 c = coef = coefficients = coeffs
1224 o = order
1226 def __init__(self, c_or_r, r=False, variable=None):
1227 if isinstance(c_or_r, poly1d):
1228 self._variable = c_or_r._variable
1229 self._coeffs = c_or_r._coeffs
1231 if set(c_or_r.__dict__) - set(self.__dict__):
1232 msg = ("In the future extra properties will not be copied "
1233 "across when constructing one poly1d from another")
1234 warnings.warn(msg, FutureWarning, stacklevel=2)
1235 self.__dict__.update(c_or_r.__dict__)
1237 if variable is not None:
1238 self._variable = variable
1239 return
1240 if r:
1241 c_or_r = poly(c_or_r)
1242 c_or_r = atleast_1d(c_or_r)
1243 if c_or_r.ndim > 1:
1244 raise ValueError("Polynomial must be 1d only.")
1245 c_or_r = trim_zeros(c_or_r, trim='f')
1246 if len(c_or_r) == 0:
1247 c_or_r = NX.array([0], dtype=c_or_r.dtype)
1248 self._coeffs = c_or_r
1249 if variable is None:
1250 variable = 'x'
1251 self._variable = variable
1253 def __array__(self, t=None):
1254 if t:
1255 return NX.asarray(self.coeffs, t)
1256 else:
1257 return NX.asarray(self.coeffs)
1259 def __repr__(self):
1260 vals = repr(self.coeffs)
1261 vals = vals[6:-1]
1262 return "poly1d(%s)" % vals
1264 def __len__(self):
1265 return self.order
1267 def __str__(self):
1268 thestr = "0"
1269 var = self.variable
1271 # Remove leading zeros
1272 coeffs = self.coeffs[NX.logical_or.accumulate(self.coeffs != 0)]
1273 N = len(coeffs)-1
1275 def fmt_float(q):
1276 s = '%.4g' % q
1277 if s.endswith('.0000'):
1278 s = s[:-5]
1279 return s
1281 for k, coeff in enumerate(coeffs):
1282 if not iscomplex(coeff):
1283 coefstr = fmt_float(real(coeff))
1284 elif real(coeff) == 0:
1285 coefstr = '%sj' % fmt_float(imag(coeff))
1286 else:
1287 coefstr = '(%s + %sj)' % (fmt_float(real(coeff)),
1288 fmt_float(imag(coeff)))
1290 power = (N-k)
1291 if power == 0:
1292 if coefstr != '0':
1293 newstr = '%s' % (coefstr,)
1294 else:
1295 if k == 0:
1296 newstr = '0'
1297 else:
1298 newstr = ''
1299 elif power == 1:
1300 if coefstr == '0':
1301 newstr = ''
1302 elif coefstr == 'b':
1303 newstr = var
1304 else:
1305 newstr = '%s %s' % (coefstr, var)
1306 else:
1307 if coefstr == '0':
1308 newstr = ''
1309 elif coefstr == 'b':
1310 newstr = '%s**%d' % (var, power,)
1311 else:
1312 newstr = '%s %s**%d' % (coefstr, var, power)
1314 if k > 0:
1315 if newstr != '':
1316 if newstr.startswith('-'):
1317 thestr = "%s - %s" % (thestr, newstr[1:])
1318 else:
1319 thestr = "%s + %s" % (thestr, newstr)
1320 else:
1321 thestr = newstr
1322 return _raise_power(thestr)
1324 def __call__(self, val):
1325 return polyval(self.coeffs, val)
1327 def __neg__(self):
1328 return poly1d(-self.coeffs)
1330 def __pos__(self):
1331 return self
1333 def __mul__(self, other):
1334 if isscalar(other):
1335 return poly1d(self.coeffs * other)
1336 else:
1337 other = poly1d(other)
1338 return poly1d(polymul(self.coeffs, other.coeffs))
1340 def __rmul__(self, other):
1341 if isscalar(other):
1342 return poly1d(other * self.coeffs)
1343 else:
1344 other = poly1d(other)
1345 return poly1d(polymul(self.coeffs, other.coeffs))
1347 def __add__(self, other):
1348 other = poly1d(other)
1349 return poly1d(polyadd(self.coeffs, other.coeffs))
1351 def __radd__(self, other):
1352 other = poly1d(other)
1353 return poly1d(polyadd(self.coeffs, other.coeffs))
1355 def __pow__(self, val):
1356 if not isscalar(val) or int(val) != val or val < 0:
1357 raise ValueError("Power to non-negative integers only.")
1358 res = [1]
1359 for _ in range(val):
1360 res = polymul(self.coeffs, res)
1361 return poly1d(res)
1363 def __sub__(self, other):
1364 other = poly1d(other)
1365 return poly1d(polysub(self.coeffs, other.coeffs))
1367 def __rsub__(self, other):
1368 other = poly1d(other)
1369 return poly1d(polysub(other.coeffs, self.coeffs))
1371 def __div__(self, other):
1372 if isscalar(other):
1373 return poly1d(self.coeffs/other)
1374 else:
1375 other = poly1d(other)
1376 return polydiv(self, other)
1378 __truediv__ = __div__
1380 def __rdiv__(self, other):
1381 if isscalar(other):
1382 return poly1d(other/self.coeffs)
1383 else:
1384 other = poly1d(other)
1385 return polydiv(other, self)
1387 __rtruediv__ = __rdiv__
1389 def __eq__(self, other):
1390 if not isinstance(other, poly1d):
1391 return NotImplemented
1392 if self.coeffs.shape != other.coeffs.shape:
1393 return False
1394 return (self.coeffs == other.coeffs).all()
1396 def __ne__(self, other):
1397 if not isinstance(other, poly1d):
1398 return NotImplemented
1399 return not self.__eq__(other)
1402 def __getitem__(self, val):
1403 ind = self.order - val
1404 if val > self.order:
1405 return self.coeffs.dtype.type(0)
1406 if val < 0:
1407 return self.coeffs.dtype.type(0)
1408 return self.coeffs[ind]
1410 def __setitem__(self, key, val):
1411 ind = self.order - key
1412 if key < 0:
1413 raise ValueError("Does not support negative powers.")
1414 if key > self.order:
1415 zr = NX.zeros(key-self.order, self.coeffs.dtype)
1416 self._coeffs = NX.concatenate((zr, self.coeffs))
1417 ind = 0
1418 self._coeffs[ind] = val
1419 return
1421 def __iter__(self):
1422 return iter(self.coeffs)
1424 def integ(self, m=1, k=0):
1425 """
1426 Return an antiderivative (indefinite integral) of this polynomial.
1428 Refer to `polyint` for full documentation.
1430 See Also
1431 --------
1432 polyint : equivalent function
1434 """
1435 return poly1d(polyint(self.coeffs, m=m, k=k))
1437 def deriv(self, m=1):
1438 """
1439 Return a derivative of this polynomial.
1441 Refer to `polyder` for full documentation.
1443 See Also
1444 --------
1445 polyder : equivalent function
1447 """
1448 return poly1d(polyder(self.coeffs, m=m))
1450# Stuff to do on module import
1452warnings.simplefilter('always', RankWarning)