codekingpro/portable-devtools
115k
1"""
2Functions to operate on polynomials.
3
4"""
5__all__ = ['poly', 'roots', 'polyint', 'polyder', 'polyadd',
6 'polysub', 'polymul', 'polydiv', 'polyval', 'poly1d',
7 'polyfit']
8
9import functools
10import re
11import warnings
12
13import numpy._core.numeric as NX
14from numpy._core import (
15 abs,
16 array,
17 atleast_1d,
18 dot,
19 finfo,
20 hstack,
21 isscalar,
22 ones,
23 overrides,
24)
25from numpy._utils import set_module
26from numpy.exceptions import RankWarning
27from numpy.lib._function_base_impl import trim_zeros
28from numpy.lib._twodim_base_impl import diag, vander
29from numpy.lib._type_check_impl import imag, iscomplex, mintypecode, real
30from numpy.linalg import eigvals, inv, lstsq
31
32array_function_dispatch = functools.partial(
33 overrides.array_function_dispatch, module='numpy')
34
35
36def _poly_dispatcher(seq_of_zeros):
37 return seq_of_zeros
38
39
40@array_function_dispatch(_poly_dispatcher)
41def poly(seq_of_zeros):
42 """
43 Find the coefficients of a polynomial with the given sequence of roots.
44
45 .. note::
46 This forms part of the old polynomial API. Since version 1.4, the
47 new polynomial API defined in `numpy.polynomial` is preferred.
48 A summary of the differences can be found in the
49 :doc:`transition guide </reference/routines.polynomials>`.
50
51 Returns the coefficients of the polynomial whose leading coefficient
52 is one for the given sequence of zeros (multiple roots must be included
53 in the sequence as many times as their multiplicity; see Examples).
54 A square matrix (or array, which will be treated as a matrix) can also
55 be given, in which case the coefficients of the characteristic polynomial
56 of the matrix are returned.
57
58 Parameters
59 ----------
60 seq_of_zeros : array_like, shape (N,) or (N, N)
61 A sequence of polynomial roots, or a square array or matrix object.
62
63 Returns
64 -------
65 c : ndarray
66 1D array of polynomial coefficients from highest to lowest degree:
67
68 ``c[0] * x**(N) + c[1] * x**(N-1) + ... + c[N-1] * x + c[N]``
69 where c[0] always equals 1.
70
71 Raises
72 ------
73 ValueError
74 If input is the wrong shape (the input must be a 1-D or square
75 2-D array).
76
77 See Also
78 --------
79 polyval : Compute polynomial values.
80 roots : Return the roots of a polynomial.
81 polyfit : Least squares polynomial fit.
82 poly1d : A one-dimensional polynomial class.
83
84 Notes
85 -----
86 Specifying the roots of a polynomial still leaves one degree of
87 freedom, typically represented by an undetermined leading
88 coefficient. [1]_ In the case of this function, that coefficient -
89 the first one in the returned array - is always taken as one. (If
90 for some reason you have one other point, the only automatic way
91 presently to leverage that information is to use ``polyfit``.)
92
93 The characteristic polynomial, :math:`p_a(t)`, of an `n`-by-`n`
94 matrix **A** is given by
95
96 :math:`p_a(t) = \\mathrm{det}(t\\, \\mathbf{I} - \\mathbf{A})`,
97
98 where **I** is the `n`-by-`n` identity matrix. [2]_
99
100 References
101 ----------
102 .. [1] M. Sullivan and M. Sullivan, III, "Algebra and Trigonometry,
103 Enhanced With Graphing Utilities," Prentice-Hall, pg. 318, 1996.
104
105 .. [2] G. Strang, "Linear Algebra and Its Applications, 2nd Edition,"
106 Academic Press, pg. 182, 1980.
107
108 Examples
109 --------
110
111 Given a sequence of a polynomial's zeros:
112
113 >>> import numpy as np
114
115 >>> np.poly((0, 0, 0)) # Multiple root example
116 array([1., 0., 0., 0.])
117
118 The line above represents z**3 + 0*z**2 + 0*z + 0.
119
120 >>> np.poly((-1./2, 0, 1./2))
121 array([ 1. , 0. , -0.25, 0. ])
122
123 The line above represents z**3 - z/4
124
125 >>> np.poly((np.random.random(1)[0], 0, np.random.random(1)[0]))
126 array([ 1. , -0.77086955, 0.08618131, 0. ]) # random
127
128 Given a square array object:
129
130 >>> P = np.array([[0, 1./3], [-1./2, 0]])
131 >>> np.poly(P)
132 array([1. , 0. , 0.16666667])
133
134 Note how in all cases the leading coefficient is always 1.
135
136 """
137 seq_of_zeros = atleast_1d(seq_of_zeros)
138 sh = seq_of_zeros.shape
139
140 if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0:
141 seq_of_zeros = eigvals(seq_of_zeros)
142 elif len(sh) == 1:
143 dt = seq_of_zeros.dtype
144 # Let object arrays slip through, e.g. for arbitrary precision
145 if dt != object:
146 seq_of_zeros = seq_of_zeros.astype(mintypecode(dt.char))
147 else:
148 raise ValueError("input must be 1d or non-empty square 2d array.")
149
150 if len(seq_of_zeros) == 0:
151 return 1.0
152 dt = seq_of_zeros.dtype
153 a = ones((1,), dtype=dt)
154 for zero in seq_of_zeros:
155 a = NX.convolve(a, array([1, -zero], dtype=dt), mode='full')
156
157 if issubclass(a.dtype.type, NX.complexfloating):
158 # if complex roots are all complex conjugates, the roots are real.
159 roots = NX.asarray(seq_of_zeros, complex)
160 if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())):
161 a = a.real.copy()
162
163 return a
164
165
166def _roots_dispatcher(p):
167 return p
168
169
170@array_function_dispatch(_roots_dispatcher)
171def roots(p):
172 """
173 Return the roots of a polynomial with coefficients given in p.
174
175 .. note::
176 This forms part of the old polynomial API. Since version 1.4, the
177 new polynomial API defined in `numpy.polynomial` is preferred.
178 A summary of the differences can be found in the
179 :doc:`transition guide </reference/routines.polynomials>`.
180
181 The values in the rank-1 array `p` are coefficients of a polynomial.
182 If the length of `p` is n+1 then the polynomial is described by::
183
184 p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n]
185
186 Parameters
187 ----------
188 p : array_like
189 Rank-1 array of polynomial coefficients.
190
191 Returns
192 -------
193 out : ndarray
194 An array containing the roots of the polynomial.
195
196 Raises
197 ------
198 ValueError
199 When `p` cannot be converted to a rank-1 array.
200
201 See also
202 --------
203 poly : Find the coefficients of a polynomial with a given sequence
204 of roots.
205 polyval : Compute polynomial values.
206 polyfit : Least squares polynomial fit.
207 poly1d : A one-dimensional polynomial class.
208
209 Notes
210 -----
211 The algorithm relies on computing the eigenvalues of the
212 companion matrix [1]_.
213
214 References
215 ----------
216 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK:
217 Cambridge University Press, 1999, pp. 146-7.
218
219 Examples
220 --------
221 >>> import numpy as np
222 >>> coeff = [3.2, 2, 1]
223 >>> np.roots(coeff)
224 array([-0.3125+0.46351241j, -0.3125-0.46351241j])
225
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.")
231
232 # find non-zero array entries
233 non_zero = NX.nonzero(NX.ravel(p))[0]
234
235 # Return an empty array if polynomial is all zeros
236 if len(non_zero) == 0:
237 return NX.array([])
238
239 # find the number of trailing zeros -- this is the number of roots at 0.
240 trailing_zeros = len(p) - non_zero[-1] - 1
241
242 # strip leading and trailing zeros
243 p = p[int(non_zero[0]):int(non_zero[-1]) + 1]
244
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)
248
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([])
257
258 # tack any zeros onto the back of the array
259 roots = hstack((roots, NX.zeros(trailing_zeros, roots.dtype)))
260 return roots
261
262
263def _polyint_dispatcher(p, m=None, k=None):
264 return (p,)
265
266
267@array_function_dispatch(_polyint_dispatcher)
268def polyint(p, m=1, k=None):
269 """
270 Return an antiderivative (indefinite integral) of a polynomial.
271
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>`.
277
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
282
283 .. math:: \\frac{k_{m-1}}{0!} x^0 + \\ldots + \\frac{k_0}{(m-1)!}x^{m-1}
284
285 of `P` so that :math:`P^{(j)}(0) = k_{m-j-1}`.
286
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.
297
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.
300
301 See Also
302 --------
303 polyder : derivative of a polynomial
304 poly1d.integ : equivalent method
305
306 Examples
307 --------
308
309 The defining property of the antiderivative:
310
311 >>> import numpy as np
312
313 >>> p = np.poly1d([1,1,1])
314 >>> P = np.polyint(p)
315 >>> P
316 poly1d([ 0.33333333, 0.5 , 1. , 0. ]) # may vary
317 >>> np.polyder(P) == p
318 True
319
320 The integration constants default to zero, but can be specified:
321
322 >>> P = np.polyint(p, 3)
323 >>> P(0)
324 0.0
325 >>> np.polyder(P)(0)
326 0.0
327 >>> np.polyder(P, 2)(0)
328 0.0
329 >>> P = np.polyint(p, 3, k=[6,5,3])
330 >>> P
331 poly1d([ 0.01666667, 0.04166667, 0.16666667, 3. , 5. , 3. ]) # may vary
332
333 Note that 3 = 6 / 2!, and that the constants are given in the order of
334 integrations. Constant of the highest-order polynomial term comes first:
335
336 >>> np.polyder(P, 2)(0)
337 6.0
338 >>> np.polyder(P, 1)(0)
339 5.0
340 >>> P(0)
341 3.0
342
343 """
344 m = int(m)
345 if m < 0:
346 raise ValueError("Order of integral must be positive (see polyder)")
347 if k is None:
348 k = NX.zeros(m, float)
349 k = atleast_1d(k)
350 if len(k) == 1 and m > 1:
351 k = k[0] * NX.ones(m, float)
352 if len(k) < m:
353 raise ValueError(
354 "k must be a scalar or a rank-1 array of length 1 or >m.")
355
356 truepoly = isinstance(p, poly1d)
357 p = NX.asarray(p)
358 if m == 0:
359 if truepoly:
360 return poly1d(p)
361 return p
362 else:
363 # Note: this must work also with object and integer arrays
364 y = NX.concatenate((p.__truediv__(NX.arange(len(p), 0, -1)), [k[0]]))
365 val = polyint(y, m - 1, k=k[1:])
366 if truepoly:
367 return poly1d(val)
368 return val
369
370
371def _polyder_dispatcher(p, m=None):
372 return (p,)
373
374
375@array_function_dispatch(_polyder_dispatcher)
376def polyder(p, m=1):
377 """
378 Return the derivative of the specified order of a polynomial.
379
380 .. note::
381 This forms part of the old polynomial API. Since version 1.4, the
382 new polynomial API defined in `numpy.polynomial` is preferred.
383 A summary of the differences can be found in the
384 :doc:`transition guide </reference/routines.polynomials>`.
385
386 Parameters
387 ----------
388 p : poly1d or sequence
389 Polynomial to differentiate.
390 A sequence is interpreted as polynomial coefficients, see `poly1d`.
391 m : int, optional
392 Order of differentiation (default: 1)
393
394 Returns
395 -------
396 der : poly1d
397 A new polynomial representing the derivative.
398
399 See Also
400 --------
401 polyint : Anti-derivative of a polynomial.
402 poly1d : Class for one-dimensional polynomials.
403
404 Examples
405 --------
406
407 The derivative of the polynomial :math:`x^3 + x^2 + x^1 + 1` is:
408
409 >>> import numpy as np
410
411 >>> p = np.poly1d([1,1,1,1])
412 >>> p2 = np.polyder(p)
413 >>> p2
414 poly1d([3, 2, 1])
415
416 which evaluates to:
417
418 >>> p2(2.)
419 17.0
420
421 We can verify this, approximating the derivative with
422 ``(f(x + h) - f(x))/h``:
423
424 >>> (p(2. + 0.001) - p(2.)) / 0.001
425 17.007000999997857
426
427 The fourth-order derivative of a 3rd-order polynomial is zero:
428
429 >>> np.polyder(p, 2)
430 poly1d([6, 2])
431 >>> np.polyder(p, 3)
432 poly1d([6])
433 >>> np.polyder(p, 4)
434 poly1d([0])
435
436 """
437 m = int(m)
438 if m < 0:
439 raise ValueError("Order of derivative must be positive (see polyint)")
440
441 truepoly = isinstance(p, poly1d)
442 p = NX.asarray(p)
443 n = len(p) - 1
444 y = p[:-1] * NX.arange(n, 0, -1)
445 if m == 0:
446 val = p
447 else:
448 val = polyder(y, m - 1)
449 if truepoly:
450 val = poly1d(val)
451 return val
452
453
454def _polyfit_dispatcher(x, y, deg, rcond=None, full=None, w=None, cov=None):
455 return (x, y, w)
456
457
458@array_function_dispatch(_polyfit_dispatcher)
459def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
460 """
461 Least squares polynomial fit.
462
463 .. note::
464 This forms part of the old polynomial API. Since version 1.4, the
465 new polynomial API defined in `numpy.polynomial` is preferred.
466 A summary of the differences can be found in the
467 :doc:`transition guide </reference/routines.polynomials>`.
468
469 Fit a polynomial ``p[0] * x**deg + ... + p[deg]`` of degree `deg`
470 to points `(x, y)`. Returns a vector of coefficients `p` that minimises
471 the squared error in the order `deg`, `deg-1`, ... `0`.
472
473 The `Polynomial.fit <numpy.polynomial.polynomial.Polynomial.fit>` class
474 method is recommended for new code as it is more stable numerically. See
475 the documentation of the method for more information.
476
477 Parameters
478 ----------
479 x : array_like, shape (M,)
480 x-coordinates of the M sample points ``(x[i], y[i])``.
481 y : array_like, shape (M,) or (M, K)
482 y-coordinates of the sample points. Several data sets of sample
483 points sharing the same x-coordinates can be fitted at once by
484 passing in a 2D-array that contains one dataset per column.
485 deg : int
486 Degree of the fitting polynomial
487 rcond : float, optional
488 Relative condition number of the fit. Singular values smaller than
489 this relative to the largest singular value will be ignored. The
490 default value is len(x)*eps, where eps is the relative precision of
491 the float type, about 2e-16 in most cases.
492 full : bool, optional
493 Switch determining nature of return value. When it is False (the
494 default) just the coefficients are returned, when True diagnostic
495 information from the singular value decomposition is also returned.
496 w : array_like, shape (M,), optional
497 Weights. If not None, the weight ``w[i]`` applies to the unsquared
498 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are
499 chosen so that the errors of the products ``w[i]*y[i]`` all have the
500 same variance. When using inverse-variance weighting, use
501 ``w[i] = 1/sigma(y[i])``. The default value is None.
502 cov : bool or str, optional
503 If given and not `False`, return not just the estimate but also its
504 covariance matrix. By default, the covariance are scaled by
505 chi2/dof, where dof = M - (deg + 1), i.e., the weights are presumed
506 to be unreliable except in a relative sense and everything is scaled
507 such that the reduced chi2 is unity. This scaling is omitted if
508 ``cov='unscaled'``, as is relevant for the case that the weights are
509 w = 1/sigma, with sigma known to be a reliable estimate of the
510 uncertainty.
511
512 Returns
513 -------
514 p : ndarray, shape (deg + 1,) or (deg + 1, K)
515 Polynomial coefficients, highest power first. If `y` was 2-D, the
516 coefficients for `k`-th data set are in ``p[:,k]``.
517
518 residuals, rank, singular_values, rcond
519 These values are only returned if ``full == True``
520
521 - residuals -- sum of squared residuals of the least squares fit
522 - rank -- the effective rank of the scaled Vandermonde
523 coefficient matrix
524 - singular_values -- singular values of the scaled Vandermonde
525 coefficient matrix
526 - rcond -- value of `rcond`.
527
528 For more details, see `numpy.linalg.lstsq`.
529
530 V : ndarray, shape (deg + 1, deg + 1) or (deg + 1, deg + 1, K)
531 Present only if ``full == False`` and ``cov == True``. The covariance
532 matrix of the polynomial coefficient estimates. The diagonal of
533 this matrix are the variance estimates for each coefficient. If y
534 is a 2-D array, then the covariance matrix for the `k`-th data set
535 are in ``V[:,:,k]``
536
537
538 Warns
539 -----
540 RankWarning
541 The rank of the coefficient matrix in the least-squares fit is
542 deficient. The warning is only raised if ``full == False``.
543
544 The warnings can be turned off by
545
546 >>> import warnings
547 >>> warnings.simplefilter('ignore', np.exceptions.RankWarning)
548
549 See Also
550 --------
551 polyval : Compute polynomial values.
552 linalg.lstsq : Computes a least-squares fit.
553 scipy.interpolate.UnivariateSpline : Computes spline fits.
554
555 Notes
556 -----
557 The solution minimizes the squared error
558
559 .. math::
560 E = \\sum_{j=0}^k |p(x_j) - y_j|^2
561
562 in the equations::
563
564 x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0]
565 x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1]
566 ...
567 x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k]
568
569 The coefficient matrix of the coefficients `p` is a Vandermonde matrix.
570
571 `polyfit` issues a `~exceptions.RankWarning` when the least-squares fit is
572 badly conditioned. This implies that the best fit is not well-defined due
573 to numerical error. The results may be improved by lowering the polynomial
574 degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter
575 can also be set to a value smaller than its default, but the resulting
576 fit may be spurious: including contributions from the small singular
577 values can add numerical noise to the result.
578
579 Note that fitting polynomial coefficients is inherently badly conditioned
580 when the degree of the polynomial is large or the interval of sample points
581 is badly centered. The quality of the fit should always be checked in these
582 cases. When polynomial fits are not satisfactory, splines may be a good
583 alternative.
584
585 References
586 ----------
587 .. [1] Wikipedia, "Curve fitting",
588 https://en.wikipedia.org/wiki/Curve_fitting
589 .. [2] Wikipedia, "Polynomial interpolation",
590 https://en.wikipedia.org/wiki/Polynomial_interpolation
591
592 Examples
593 --------
594 >>> import numpy as np
595 >>> import warnings
596 >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0])
597 >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
598 >>> z = np.polyfit(x, y, 3)
599 >>> z
600 array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) # may vary
601
602 It is convenient to use `poly1d` objects for dealing with polynomials:
603
604 >>> p = np.poly1d(z)
605 >>> p(0.5)
606 0.6143849206349179 # may vary
607 >>> p(3.5)
608 -0.34732142857143039 # may vary
609 >>> p(10)
610 22.579365079365115 # may vary
611
612 High-order polynomials may oscillate wildly:
613
614 >>> with warnings.catch_warnings():
615 ... warnings.simplefilter('ignore', np.exceptions.RankWarning)
616 ... p30 = np.poly1d(np.polyfit(x, y, 30))
617 ...
618 >>> p30(4)
619 -0.80000000000000204 # may vary
620 >>> p30(5)
621 -0.99999999999999445 # may vary
622 >>> p30(4.5)
623 -0.10547061179440398 # may vary
624
625 Illustration:
626
627 >>> import matplotlib.pyplot as plt
628 >>> xp = np.linspace(-2, 6, 100)
629 >>> _ = plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--')
630 >>> plt.ylim(-2,2)
631 (-2, 2)
632 >>> plt.show()
633
634 """
635 order = int(deg) + 1
636 x = NX.asarray(x) + 0.0
637 y = NX.asarray(y) + 0.0
638
639 # check arguments.
640 if deg < 0:
641 raise ValueError("expected deg >= 0")
642 if x.ndim != 1:
643 raise TypeError("expected 1D vector for x")
644 if x.size == 0:
645 raise TypeError("expected non-empty vector for x")
646 if y.ndim < 1 or y.ndim > 2:
647 raise TypeError("expected 1D or 2D array for y")
648 if x.shape[0] != y.shape[0]:
649 raise TypeError("expected x and y to have same length")
650
651 # set rcond
652 if rcond is None:
653 rcond = len(x) * finfo(x.dtype).eps
654
655 # set up least squares equation for powers of x
656 lhs = vander(x, order)
657 rhs = y
658
659 # apply weighting
660 if w is not None:
661 w = NX.asarray(w) + 0.0
662 if w.ndim != 1:
663 raise TypeError("expected a 1-d array for weights")
664 if w.shape[0] != y.shape[0]:
665 raise TypeError("expected w and y to have the same length")
666 lhs *= w[:, NX.newaxis]
667 if rhs.ndim == 2:
668 rhs *= w[:, NX.newaxis]
669 else:
670 rhs *= w
671
672 # scale lhs to improve condition number and solve
673 scale = NX.sqrt((lhs * lhs).sum(axis=0))
674 lhs /= scale
675 c, resids, rank, s = lstsq(lhs, rhs, rcond)
676 c = (c.T / scale).T # broadcast scale coefficients
677
678 # warn on rank reduction, which indicates an ill conditioned matrix
679 if rank != order and not full:
680 msg = "Polyfit may be poorly conditioned"
681 warnings.warn(msg, RankWarning, stacklevel=2)
682
683 if full:
684 return c, resids, rank, s, rcond
685 elif cov:
686 Vbase = inv(dot(lhs.T, lhs))
687 Vbase /= NX.outer(scale, scale)
688 if cov == "unscaled":
689 fac = 1
690 else:
691 if len(x) <= order:
692 raise ValueError("the number of data points must exceed order "
693 "to scale the covariance matrix")
694 # note, this used to be: fac = resids / (len(x) - order - 2.0)
695 # it was decided that the "- 2" (originally justified by "Bayesian
696 # uncertainty analysis") is not what the user expects
697 # (see gh-11196 and gh-11197)
698 fac = resids / (len(x) - order)
699 if y.ndim == 1:
700 return c, Vbase * fac
701 else:
702 return c, Vbase[:, :, NX.newaxis] * fac
703 else:
704 return c
705
706
707def _polyval_dispatcher(p, x):
708 return (p, x)
709
710
711@array_function_dispatch(_polyval_dispatcher)
712def polyval(p, x):
713 """
714 Evaluate a polynomial at specific values.
715
716 .. note::
717 This forms part of the old polynomial API. Since version 1.4, the
718 new polynomial API defined in `numpy.polynomial` is preferred.
719 A summary of the differences can be found in the
720 :doc:`transition guide </reference/routines.polynomials>`.
721
722 If `p` is of length N, this function returns the value::
723
724 p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]
725
726 If `x` is a sequence, then ``p(x)`` is returned for each element of ``x``.
727 If `x` is another polynomial then the composite polynomial ``p(x(t))``
728 is returned.
729
730 Parameters
731 ----------
732 p : array_like or poly1d object
733 1D array of polynomial coefficients (including coefficients equal
734 to zero) from highest degree to the constant term, or an
735 instance of poly1d.
736 x : array_like or poly1d object
737 A number, an array of numbers, or an instance of poly1d, at
738 which to evaluate `p`.
739
740 Returns
741 -------
742 values : ndarray or poly1d
743 If `x` is a poly1d instance, the result is the composition of the two
744 polynomials, i.e., `x` is "substituted" in `p` and the simplified
745 result is returned. In addition, the type of `x` - array_like or
746 poly1d - governs the type of the output: `x` array_like => `values`
747 array_like, `x` a poly1d object => `values` is also.
748
749 See Also
750 --------
751 poly1d: A polynomial class.
752
753 Notes
754 -----
755 Horner's scheme [1]_ is used to evaluate the polynomial. Even so,
756 for polynomials of high degree the values may be inaccurate due to
757 rounding errors. Use carefully.
758
759 If `x` is a subtype of `ndarray` the return value will be of the same type.
760
761 References
762 ----------
763 .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng.
764 trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand
765 Reinhold Co., 1985, pg. 720.
766
767 Examples
768 --------
769 >>> import numpy as np
770 >>> np.polyval([3,0,1], 5) # 3 * 5**2 + 0 * 5**1 + 1
771 76
772 >>> np.polyval([3,0,1], np.poly1d(5))
773 poly1d([76])
774 >>> np.polyval(np.poly1d([3,0,1]), 5)
775 76
776 >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5))
777 poly1d([76])
778
779 """
780 p = NX.asarray(p)
781 if isinstance(x, poly1d):
782 y = 0
783 else:
784 x = NX.asanyarray(x)
785 y = NX.zeros_like(x)
786 for pv in p:
787 y = y * x + pv
788 return y
789
790
791def _binary_op_dispatcher(a1, a2):
792 return (a1, a2)
793
794
795@array_function_dispatch(_binary_op_dispatcher)
796def polyadd(a1, a2):
797 """
798 Find the sum of two polynomials.
799
800 .. note::
801 This forms part of the old polynomial API. Since version 1.4, the
802 new polynomial API defined in `numpy.polynomial` is preferred.
803 A summary of the differences can be found in the
804 :doc:`transition guide </reference/routines.polynomials>`.
805
806 Returns the polynomial resulting from the sum of two input polynomials.
807 Each input must be either a poly1d object or a 1D sequence of polynomial
808 coefficients, from highest to lowest degree.
809
810 Parameters
811 ----------
812 a1, a2 : array_like or poly1d object
813 Input polynomials.
814
815 Returns
816 -------
817 out : ndarray or poly1d object
818 The sum of the inputs. If either input is a poly1d object, then the
819 output is also a poly1d object. Otherwise, it is a 1D array of
820 polynomial coefficients from highest to lowest degree.
821
822 See Also
823 --------
824 poly1d : A one-dimensional polynomial class.
825 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
826
827 Examples
828 --------
829 >>> import numpy as np
830 >>> np.polyadd([1, 2], [9, 5, 4])
831 array([9, 6, 6])
832
833 Using poly1d objects:
834
835 >>> p1 = np.poly1d([1, 2])
836 >>> p2 = np.poly1d([9, 5, 4])
837 >>> print(p1)
838 1 x + 2
839 >>> print(p2)
840 2
841 9 x + 5 x + 4
842 >>> print(np.polyadd(p1, p2))
843 2
844 9 x + 6 x + 6
845
846 """
847 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
848 a1 = atleast_1d(a1)
849 a2 = atleast_1d(a2)
850 diff = len(a2) - len(a1)
851 if diff == 0:
852 val = a1 + a2
853 elif diff > 0:
854 zr = NX.zeros(diff, a1.dtype)
855 val = NX.concatenate((zr, a1)) + a2
856 else:
857 zr = NX.zeros(abs(diff), a2.dtype)
858 val = a1 + NX.concatenate((zr, a2))
859 if truepoly:
860 val = poly1d(val)
861 return val
862
863
864@array_function_dispatch(_binary_op_dispatcher)
865def polysub(a1, a2):
866 """
867 Difference (subtraction) of two polynomials.
868
869 .. note::
870 This forms part of the old polynomial API. Since version 1.4, the
871 new polynomial API defined in `numpy.polynomial` is preferred.
872 A summary of the differences can be found in the
873 :doc:`transition guide </reference/routines.polynomials>`.
874
875 Given two polynomials `a1` and `a2`, returns ``a1 - a2``.
876 `a1` and `a2` can be either array_like sequences of the polynomials'
877 coefficients (including coefficients equal to zero), or `poly1d` objects.
878
879 Parameters
880 ----------
881 a1, a2 : array_like or poly1d
882 Minuend and subtrahend polynomials, respectively.
883
884 Returns
885 -------
886 out : ndarray or poly1d
887 Array or `poly1d` object of the difference polynomial's coefficients.
888
889 See Also
890 --------
891 polyval, polydiv, polymul, polyadd
892
893 Examples
894 --------
895
896 .. math:: (2 x^2 + 10 x - 2) - (3 x^2 + 10 x -4) = (-x^2 + 2)
897
898 >>> import numpy as np
899
900 >>> np.polysub([2, 10, -2], [3, 10, -4])
901 array([-1, 0, 2])
902
903 """
904 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
905 a1 = atleast_1d(a1)
906 a2 = atleast_1d(a2)
907 diff = len(a2) - len(a1)
908 if diff == 0:
909 val = a1 - a2
910 elif diff > 0:
911 zr = NX.zeros(diff, a1.dtype)
912 val = NX.concatenate((zr, a1)) - a2
913 else:
914 zr = NX.zeros(abs(diff), a2.dtype)
915 val = a1 - NX.concatenate((zr, a2))
916 if truepoly:
917 val = poly1d(val)
918 return val
919
920
921@array_function_dispatch(_binary_op_dispatcher)
922def polymul(a1, a2):
923 """
924 Find the product of two polynomials.
925
926 .. note::
927 This forms part of the old polynomial API. Since version 1.4, the
928 new polynomial API defined in `numpy.polynomial` is preferred.
929 A summary of the differences can be found in the
930 :doc:`transition guide </reference/routines.polynomials>`.
931
932 Finds the polynomial resulting from the multiplication of the two input
933 polynomials. Each input must be either a poly1d object or a 1D sequence
934 of polynomial coefficients, from highest to lowest degree.
935
936 Parameters
937 ----------
938 a1, a2 : array_like or poly1d object
939 Input polynomials.
940
941 Returns
942 -------
943 out : ndarray or poly1d object
944 The polynomial resulting from the multiplication of the inputs. If
945 either inputs is a poly1d object, then the output is also a poly1d
946 object. Otherwise, it is a 1D array of polynomial coefficients from
947 highest to lowest degree.
948
949 See Also
950 --------
951 poly1d : A one-dimensional polynomial class.
952 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval
953 convolve : Array convolution. Same output as polymul, but has parameter
954 for overlap mode.
955
956 Examples
957 --------
958 >>> import numpy as np
959 >>> np.polymul([1, 2, 3], [9, 5, 1])
960 array([ 9, 23, 38, 17, 3])
961
962 Using poly1d objects:
963
964 >>> p1 = np.poly1d([1, 2, 3])
965 >>> p2 = np.poly1d([9, 5, 1])
966 >>> print(p1)
967 2
968 1 x + 2 x + 3
969 >>> print(p2)
970 2
971 9 x + 5 x + 1
972 >>> print(np.polymul(p1, p2))
973 4 3 2
974 9 x + 23 x + 38 x + 17 x + 3
975
976 """
977 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
978 a1, a2 = poly1d(a1), poly1d(a2)
979 val = NX.convolve(a1, a2)
980 if truepoly:
981 val = poly1d(val)
982 return val
983
984
985def _polydiv_dispatcher(u, v):
986 return (u, v)
987
988
989@array_function_dispatch(_polydiv_dispatcher)
990def polydiv(u, v):
991 """
992 Returns the quotient and remainder of polynomial division.
993
994 .. note::
995 This forms part of the old polynomial API. Since version 1.4, the
996 new polynomial API defined in `numpy.polynomial` is preferred.
997 A summary of the differences can be found in the
998 :doc:`transition guide </reference/routines.polynomials>`.
999
1000 The input arrays are the coefficients (including any coefficients
1001 equal to zero) of the "numerator" (dividend) and "denominator"
1002 (divisor) polynomials, respectively.
1003
1004 Parameters
1005 ----------
1006 u : array_like or poly1d
1007 Dividend polynomial's coefficients.
1008
1009 v : array_like or poly1d
1010 Divisor polynomial's coefficients.
1011
1012 Returns
1013 -------
1014 q : ndarray
1015 Coefficients, including those equal to zero, of the quotient.
1016 r : ndarray
1017 Coefficients, including those equal to zero, of the remainder.
1018
1019 See Also
1020 --------
1021 poly, polyadd, polyder, polydiv, polyfit, polyint, polymul, polysub
1022 polyval
1023
1024 Notes
1025 -----
1026 Both `u` and `v` must be 0-d or 1-d (ndim = 0 or 1), but `u.ndim` need
1027 not equal `v.ndim`. In other words, all four possible combinations -
1028 ``u.ndim = v.ndim = 0``, ``u.ndim = v.ndim = 1``,
1029 ``u.ndim = 1, v.ndim = 0``, and ``u.ndim = 0, v.ndim = 1`` - work.
1030
1031 Examples
1032 --------
1033
1034 .. math:: \\frac{3x^2 + 5x + 2}{2x + 1} = 1.5x + 1.75, remainder 0.25
1035
1036 >>> import numpy as np
1037
1038 >>> x = np.array([3.0, 5.0, 2.0])
1039 >>> y = np.array([2.0, 1.0])
1040 >>> np.polydiv(x, y)
1041 (array([1.5 , 1.75]), array([0.25]))
1042
1043 """
1044 truepoly = (isinstance(u, poly1d) or isinstance(v, poly1d))
1045 u = atleast_1d(u) + 0.0
1046 v = atleast_1d(v) + 0.0
1047 # w has the common type
1048 w = u[0] + v[0]
1049 m = len(u) - 1
1050 n = len(v) - 1
1051 scale = 1. / v[0]
1052 q = NX.zeros((max(m - n + 1, 1),), w.dtype)
1053 r = u.astype(w.dtype)
1054 for k in range(m - n + 1):
1055 d = scale * r[k]
1056 q[k] = d
1057 r[k:k + n + 1] -= d * v
1058 while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
1059 r = r[1:]
1060 if truepoly:
1061 return poly1d(q), poly1d(r)
1062 return q, r
1063
1064
1065_poly_mat = re.compile(r"\*\*([0-9]*)")
1066def _raise_power(astr, wrap=70):
1067 n = 0
1068 line1 = ''
1069 line2 = ''
1070 output = ' '
1071 while True:
1072 mat = _poly_mat.search(astr, n)
1073 if mat is None:
1074 break
1075 span = mat.span()
1076 power = mat.groups()[0]
1077 partstr = astr[n:span[0]]
1078 n = span[1]
1079 toadd2 = partstr + ' ' * (len(power) - 1)
1080 toadd1 = ' ' * (len(partstr) - 1) + power
1081 if ((len(line2) + len(toadd2) > wrap) or
1082 (len(line1) + len(toadd1) > wrap)):
1083 output += line1 + "\n" + line2 + "\n "
1084 line1 = toadd1
1085 line2 = toadd2
1086 else:
1087 line2 += partstr + ' ' * (len(power) - 1)
1088 line1 += ' ' * (len(partstr) - 1) + power
1089 output += line1 + "\n" + line2
1090 return output + astr[n:]
1091
1092
1093@set_module('numpy')
1094class poly1d:
1095 """
1096 A one-dimensional polynomial class.
1097
1098 .. note::
1099 This forms part of the old polynomial API. Since version 1.4, the
1100 new polynomial API defined in `numpy.polynomial` is preferred.
1101 A summary of the differences can be found in the
1102 :doc:`transition guide </reference/routines.polynomials>`.
1103
1104 A convenience class, used to encapsulate "natural" operations on
1105 polynomials so that said operations may take on their customary
1106 form in code (see Examples).
1107
1108 Parameters
1109 ----------
1110 c_or_r : array_like
1111 The polynomial's coefficients, in decreasing powers, or if
1112 the value of the second parameter is True, the polynomial's
1113 roots (values where the polynomial evaluates to 0). For example,
1114 ``poly1d([1, 2, 3])`` returns an object that represents
1115 :math:`x^2 + 2x + 3`, whereas ``poly1d([1, 2, 3], True)`` returns
1116 one that represents :math:`(x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x -6`.
1117 r : bool, optional
1118 If True, `c_or_r` specifies the polynomial's roots; the default
1119 is False.
1120 variable : str, optional
1121 Changes the variable used when printing `p` from `x` to `variable`
1122 (see Examples).
1123
1124 Examples
1125 --------
1126 >>> import numpy as np
1127
1128 Construct the polynomial :math:`x^2 + 2x + 3`:
1129
1130 >>> import numpy as np
1131
1132 >>> p = np.poly1d([1, 2, 3])
1133 >>> print(np.poly1d(p))
1134 2
1135 1 x + 2 x + 3
1136
1137 Evaluate the polynomial at :math:`x = 0.5`:
1138
1139 >>> p(0.5)
1140 4.25
1141
1142 Find the roots:
1143
1144 >>> p.r
1145 array([-1.+1.41421356j, -1.-1.41421356j])
1146 >>> p(p.r)
1147 array([ -4.44089210e-16+0.j, -4.44089210e-16+0.j]) # may vary
1148
1149 These numbers in the previous line represent (0, 0) to machine precision
1150
1151 Show the coefficients:
1152
1153 >>> p.c
1154 array([1, 2, 3])
1155
1156 Display the order (the leading zero-coefficients are removed):
1157
1158 >>> p.order
1159 2
1160
1161 Show the coefficient of the k-th power in the polynomial
1162 (which is equivalent to ``p.c[-(i+1)]``):
1163
1164 >>> p[1]
1165 2
1166
1167 Polynomials can be added, subtracted, multiplied, and divided
1168 (returns quotient and remainder):
1169
1170 >>> p * p
1171 poly1d([ 1, 4, 10, 12, 9])
1172
1173 >>> (p**3 + 4) / p
1174 (poly1d([ 1., 4., 10., 12., 9.]), poly1d([4.]))
1175
1176 ``asarray(p)`` gives the coefficient array, so polynomials can be
1177 used in all functions that accept arrays:
1178
1179 >>> p**2 # square of polynomial
1180 poly1d([ 1, 4, 10, 12, 9])
1181
1182 >>> np.square(p) # square of individual coefficients
1183 array([1, 4, 9])
1184
1185 The variable used in the string representation of `p` can be modified,
1186 using the `variable` parameter:
1187
1188 >>> p = np.poly1d([1,2,3], variable='z')
1189 >>> print(p)
1190 2
1191 1 z + 2 z + 3
1192
1193 Construct a polynomial from its roots:
1194
1195 >>> np.poly1d([1, 2], True)
1196 poly1d([ 1., -3., 2.])
1197
1198 This is the same polynomial as obtained by:
1199
1200 >>> np.poly1d([1, -1]) * np.poly1d([1, -2])