Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/numpy/lib/polynomial.py: 21%

1""" 

2Functions to operate on polynomials. 

3 

4""" 

5__all__ = ['poly', 'roots', 'polyint', 'polyder', 'polyadd', 

6 'polysub', 'polymul', 'polydiv', 'polyval', 'poly1d', 

7 'polyfit', 'RankWarning'] 

8 

9import functools 

10import re 

11import warnings 

12 

13from .._utils import set_module 

14import numpy.core.numeric as NX 

15 

16from numpy.core import (isscalar, abs, finfo, atleast_1d, hstack, dot, array, 

17 ones) 

18from numpy.core import overrides 

19from numpy.lib.twodim_base import diag, vander 

20from numpy.lib.function_base import trim_zeros 

21from numpy.lib.type_check import iscomplex, real, imag, mintypecode 

22from numpy.linalg import eigvals, lstsq, inv 

23 

24 

25array_function_dispatch = functools.partial( 

26 overrides.array_function_dispatch, module='numpy') 

27 

28 

29@set_module('numpy') 

30class RankWarning(UserWarning): 

31 """ 

32 Issued by `polyfit` when the Vandermonde matrix is rank deficient. 

33 

34 For more information, a way to suppress the warning, and an example of 

35 `RankWarning` being issued, see `polyfit`. 

36 

37 """ 

38 pass 

39 

40 

41def _poly_dispatcher(seq_of_zeros): 

42 return seq_of_zeros 

43 

44 

45@array_function_dispatch(_poly_dispatcher) 

46def poly(seq_of_zeros): 

47 """ 

48 Find the coefficients of a polynomial with the given sequence of roots. 

49 

50 .. note:: 

51 This forms part of the old polynomial API. Since version 1.4, the 

52 new polynomial API defined in `numpy.polynomial` is preferred. 

53 A summary of the differences can be found in the 

54 :doc:`transition guide </reference/routines.polynomials>`. 

55 

56 Returns the coefficients of the polynomial whose leading coefficient 

57 is one for the given sequence of zeros (multiple roots must be included 

58 in the sequence as many times as their multiplicity; see Examples). 

59 A square matrix (or array, which will be treated as a matrix) can also 

60 be given, in which case the coefficients of the characteristic polynomial 

61 of the matrix are returned. 

62 

63 Parameters 

64 ---------- 

65 seq_of_zeros : array_like, shape (N,) or (N, N) 

66 A sequence of polynomial roots, or a square array or matrix object. 

67 

68 Returns 

69 ------- 

70 c : ndarray 

71 1D array of polynomial coefficients from highest to lowest degree: 

72 

73 ``c[0] * x**(N) + c[1] * x**(N-1) + ... + c[N-1] * x + c[N]`` 

74 where c[0] always equals 1. 

75 

76 Raises 

77 ------ 

78 ValueError 

79 If input is the wrong shape (the input must be a 1-D or square 

80 2-D array). 

81 

82 See Also 

83 -------- 

84 polyval : Compute polynomial values. 

85 roots : Return the roots of a polynomial. 

86 polyfit : Least squares polynomial fit. 

87 poly1d : A one-dimensional polynomial class. 

88 

89 Notes 

90 ----- 

91 Specifying the roots of a polynomial still leaves one degree of 

92 freedom, typically represented by an undetermined leading 

93 coefficient. [1]_ In the case of this function, that coefficient - 

94 the first one in the returned array - is always taken as one. (If 

95 for some reason you have one other point, the only automatic way 

96 presently to leverage that information is to use ``polyfit``.) 

97 

98 The characteristic polynomial, :math:`p_a(t)`, of an `n`-by-`n` 

99 matrix **A** is given by 

100 

101 :math:`p_a(t) = \\mathrm{det}(t\\, \\mathbf{I} - \\mathbf{A})`, 

102 

103 where **I** is the `n`-by-`n` identity matrix. [2]_ 

104 

105 References 

106 ---------- 

107 .. [1] M. Sullivan and M. Sullivan, III, "Algebra and Trigonometry, 

108 Enhanced With Graphing Utilities," Prentice-Hall, pg. 318, 1996. 

109 

110 .. [2] G. Strang, "Linear Algebra and Its Applications, 2nd Edition," 

111 Academic Press, pg. 182, 1980. 

112 

113 Examples 

114 -------- 

115 Given a sequence of a polynomial's zeros: 

116 

117 >>> np.poly((0, 0, 0)) # Multiple root example 

118 array([1., 0., 0., 0.]) 

119 

120 The line above represents z**3 + 0*z**2 + 0*z + 0. 

121 

122 >>> np.poly((-1./2, 0, 1./2)) 

123 array([ 1. , 0. , -0.25, 0. ]) 

124 

125 The line above represents z**3 - z/4 

126 

127 >>> np.poly((np.random.random(1)[0], 0, np.random.random(1)[0])) 

128 array([ 1. , -0.77086955, 0.08618131, 0. ]) # random 

129 

130 Given a square array object: 

131 

132 >>> P = np.array([[0, 1./3], [-1./2, 0]]) 

133 >>> np.poly(P) 

134 array([1. , 0. , 0.16666667]) 

135 

136 Note how in all cases the leading coefficient is always 1. 

137 

138 """ 

139 seq_of_zeros = atleast_1d(seq_of_zeros) 

140 sh = seq_of_zeros.shape 

141 

142 if len(sh) == 2 and sh[0] == sh[1] and sh[0] != 0: 

143 seq_of_zeros = eigvals(seq_of_zeros) 

144 elif len(sh) == 1: 

145 dt = seq_of_zeros.dtype 

146 # Let object arrays slip through, e.g. for arbitrary precision 

147 if dt != object: 

148 seq_of_zeros = seq_of_zeros.astype(mintypecode(dt.char)) 

149 else: 

150 raise ValueError("input must be 1d or non-empty square 2d array.") 

151 

152 if len(seq_of_zeros) == 0: 

153 return 1.0 

154 dt = seq_of_zeros.dtype 

155 a = ones((1,), dtype=dt) 

156 for zero in seq_of_zeros: 

157 a = NX.convolve(a, array([1, -zero], dtype=dt), mode='full') 

158 

159 if issubclass(a.dtype.type, NX.complexfloating): 

160 # if complex roots are all complex conjugates, the roots are real. 

161 roots = NX.asarray(seq_of_zeros, complex) 

162 if NX.all(NX.sort(roots) == NX.sort(roots.conjugate())): 

163 a = a.real.copy() 

164 

165 return a 

166 

167 

168def _roots_dispatcher(p): 

169 return p 

170 

171 

172@array_function_dispatch(_roots_dispatcher) 

173def roots(p): 

174 """ 

175 Return the roots of a polynomial with coefficients given in p. 

176 

177 .. note:: 

178 This forms part of the old polynomial API. Since version 1.4, the 

179 new polynomial API defined in `numpy.polynomial` is preferred. 

180 A summary of the differences can be found in the 

181 :doc:`transition guide </reference/routines.polynomials>`. 

182 

183 The values in the rank-1 array `p` are coefficients of a polynomial. 

184 If the length of `p` is n+1 then the polynomial is described by:: 

185 

186 p[0] * x**n + p[1] * x**(n-1) + ... + p[n-1]*x + p[n] 

187 

188 Parameters 

189 ---------- 

190 p : array_like 

191 Rank-1 array of polynomial coefficients. 

192 

193 Returns 

194 ------- 

195 out : ndarray 

196 An array containing the roots of the polynomial. 

197 

198 Raises 

199 ------ 

200 ValueError 

201 When `p` cannot be converted to a rank-1 array. 

202 

203 See also 

204 -------- 

205 poly : Find the coefficients of a polynomial with a given sequence 

206 of roots. 

207 polyval : Compute polynomial values. 

208 polyfit : Least squares polynomial fit. 

209 poly1d : A one-dimensional polynomial class. 

210 

211 Notes 

212 ----- 

213 The algorithm relies on computing the eigenvalues of the 

214 companion matrix [1]_. 

215 

216 References 

217 ---------- 

218 .. [1] R. A. Horn & C. R. Johnson, *Matrix Analysis*. Cambridge, UK: 

219 Cambridge University Press, 1999, pp. 146-7. 

220 

221 Examples 

222 -------- 

223 >>> coeff = [3.2, 2, 1] 

224 >>> np.roots(coeff) 

225 array([-0.3125+0.46351241j, -0.3125-0.46351241j]) 

226 

227 """ 

228 # If input is scalar, this makes it an array 

229 p = atleast_1d(p) 

230 if p.ndim != 1: 

231 raise ValueError("Input must be a rank-1 array.") 

232 

233 # find non-zero array entries 

234 non_zero = NX.nonzero(NX.ravel(p))[0] 

235 

236 # Return an empty array if polynomial is all zeros 

237 if len(non_zero) == 0: 

238 return NX.array([]) 

239 

240 # find the number of trailing zeros -- this is the number of roots at 0. 

241 trailing_zeros = len(p) - non_zero[-1] - 1 

242 

243 # strip leading and trailing zeros 

244 p = p[int(non_zero[0]):int(non_zero[-1])+1] 

245 

246 # casting: if incoming array isn't floating point, make it floating point. 

247 if not issubclass(p.dtype.type, (NX.floating, NX.complexfloating)): 

248 p = p.astype(float) 

249 

250 N = len(p) 

251 if N > 1: 

252 # build companion matrix and find its eigenvalues (the roots) 

253 A = diag(NX.ones((N-2,), p.dtype), -1) 

254 A[0,:] = -p[1:] / p[0] 

255 roots = eigvals(A) 

256 else: 

257 roots = NX.array([]) 

258 

259 # tack any zeros onto the back of the array 

260 roots = hstack((roots, NX.zeros(trailing_zeros, roots.dtype))) 

261 return roots 

262 

263 

264def _polyint_dispatcher(p, m=None, k=None): 

265 return (p,) 

266 

267 

268@array_function_dispatch(_polyint_dispatcher) 

269def polyint(p, m=1, k=None): 

270 """ 

271 Return an antiderivative (indefinite integral) of a polynomial. 

272 

273 .. note:: 

274 This forms part of the old polynomial API. Since version 1.4, the 

275 new polynomial API defined in `numpy.polynomial` is preferred. 

276 A summary of the differences can be found in the 

277 :doc:`transition guide </reference/routines.polynomials>`. 

278 

279 The returned order `m` antiderivative `P` of polynomial `p` satisfies 

280 :math:`\\frac{d^m}{dx^m}P(x) = p(x)` and is defined up to `m - 1` 

281 integration constants `k`. The constants determine the low-order 

282 polynomial part 

283 

284 .. math:: \\frac{k_{m-1}}{0!} x^0 + \\ldots + \\frac{k_0}{(m-1)!}x^{m-1} 

285 

286 of `P` so that :math:`P^{(j)}(0) = k_{m-j-1}`. 

287 

288 Parameters 

289 ---------- 

290 p : array_like or poly1d 

291 Polynomial to integrate. 

292 A sequence is interpreted as polynomial coefficients, see `poly1d`. 

293 m : int, optional 

294 Order of the antiderivative. (Default: 1) 

295 k : list of `m` scalars or scalar, optional 

296 Integration constants. They are given in the order of integration: 

297 those corresponding to highest-order terms come first. 

298 

299 If ``None`` (default), all constants are assumed to be zero. 

300 If `m = 1`, a single scalar can be given instead of a list. 

301 

302 See Also 

303 -------- 

304 polyder : derivative of a polynomial 

305 poly1d.integ : equivalent method 

306 

307 Examples 

308 -------- 

309 The defining property of the antiderivative: 

310 

311 >>> p = np.poly1d([1,1,1]) 

312 >>> P = np.polyint(p) 

313 >>> P 

314 poly1d([ 0.33333333, 0.5 , 1. , 0. ]) # may vary 

315 >>> np.polyder(P) == p 

316 True 

317 

318 The integration constants default to zero, but can be specified: 

319 

320 >>> P = np.polyint(p, 3) 

321 >>> P(0) 

322 0.0 

323 >>> np.polyder(P)(0) 

324 0.0 

325 >>> np.polyder(P, 2)(0) 

326 0.0 

327 >>> P = np.polyint(p, 3, k=[6,5,3]) 

328 >>> P 

329 poly1d([ 0.01666667, 0.04166667, 0.16666667, 3. , 5. , 3. ]) # may vary 

330 

331 Note that 3 = 6 / 2!, and that the constants are given in the order of 

332 integrations. Constant of the highest-order polynomial term comes first: 

333 

334 >>> np.polyder(P, 2)(0) 

335 6.0 

336 >>> np.polyder(P, 1)(0) 

337 5.0 

338 >>> P(0) 

339 3.0 

340 

341 """ 

342 m = int(m) 

343 if m < 0: 

344 raise ValueError("Order of integral must be positive (see polyder)") 

345 if k is None: 

346 k = NX.zeros(m, float) 

347 k = atleast_1d(k) 

348 if len(k) == 1 and m > 1: 

349 k = k[0]*NX.ones(m, float) 

350 if len(k) < m: 

351 raise ValueError( 

352 "k must be a scalar or a rank-1 array of length 1 or >m.") 

353 

354 truepoly = isinstance(p, poly1d) 

355 p = NX.asarray(p) 

356 if m == 0: 

357 if truepoly: 

358 return poly1d(p) 

359 return p 

360 else: 

361 # Note: this must work also with object and integer arrays 

362 y = NX.concatenate((p.__truediv__(NX.arange(len(p), 0, -1)), [k[0]])) 

363 val = polyint(y, m - 1, k=k[1:]) 

364 if truepoly: 

365 return poly1d(val) 

366 return val 

367 

368 

369def _polyder_dispatcher(p, m=None): 

370 return (p,) 

371 

372 

373@array_function_dispatch(_polyder_dispatcher) 

374def polyder(p, m=1): 

375 """ 

376 Return the derivative of the specified order of a polynomial. 

377 

378 .. note:: 

379 This forms part of the old polynomial API. Since version 1.4, the 

380 new polynomial API defined in `numpy.polynomial` is preferred. 

381 A summary of the differences can be found in the 

382 :doc:`transition guide </reference/routines.polynomials>`. 

383 

384 Parameters 

385 ---------- 

386 p : poly1d or sequence 

387 Polynomial to differentiate. 

388 A sequence is interpreted as polynomial coefficients, see `poly1d`. 

389 m : int, optional 

390 Order of differentiation (default: 1) 

391 

392 Returns 

393 ------- 

394 der : poly1d 

395 A new polynomial representing the derivative. 

396 

397 See Also 

398 -------- 

399 polyint : Anti-derivative of a polynomial. 

400 poly1d : Class for one-dimensional polynomials. 

401 

402 Examples 

403 -------- 

404 The derivative of the polynomial :math:`x^3 + x^2 + x^1 + 1` is: 

405 

406 >>> p = np.poly1d([1,1,1,1]) 

407 >>> p2 = np.polyder(p) 

408 >>> p2 

409 poly1d([3, 2, 1]) 

410 

411 which evaluates to: 

412 

413 >>> p2(2.) 

414 17.0 

415 

416 We can verify this, approximating the derivative with 

417 ``(f(x + h) - f(x))/h``: 

418 

419 >>> (p(2. + 0.001) - p(2.)) / 0.001 

420 17.007000999997857 

421 

422 The fourth-order derivative of a 3rd-order polynomial is zero: 

423 

424 >>> np.polyder(p, 2) 

425 poly1d([6, 2]) 

426 >>> np.polyder(p, 3) 

427 poly1d([6]) 

428 >>> np.polyder(p, 4) 

429 poly1d([0]) 

430 

431 """ 

432 m = int(m) 

433 if m < 0: 

434 raise ValueError("Order of derivative must be positive (see polyint)") 

435 

436 truepoly = isinstance(p, poly1d) 

437 p = NX.asarray(p) 

438 n = len(p) - 1 

439 y = p[:-1] * NX.arange(n, 0, -1) 

440 if m == 0: 

441 val = p 

442 else: 

443 val = polyder(y, m - 1) 

444 if truepoly: 

445 val = poly1d(val) 

446 return val 

447 

448 

449def _polyfit_dispatcher(x, y, deg, rcond=None, full=None, w=None, cov=None): 

450 return (x, y, w) 

451 

452 

453@array_function_dispatch(_polyfit_dispatcher) 

454def polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False): 

455 """ 

456 Least squares polynomial fit. 

457 

458 .. note:: 

459 This forms part of the old polynomial API. Since version 1.4, the 

460 new polynomial API defined in `numpy.polynomial` is preferred. 

461 A summary of the differences can be found in the 

462 :doc:`transition guide </reference/routines.polynomials>`. 

463 

464 Fit a polynomial ``p(x) = p[0] * x**deg + ... + p[deg]`` of degree `deg` 

465 to points `(x, y)`. Returns a vector of coefficients `p` that minimises 

466 the squared error in the order `deg`, `deg-1`, ... `0`. 

467 

468 The `Polynomial.fit <numpy.polynomial.polynomial.Polynomial.fit>` class 

469 method is recommended for new code as it is more stable numerically. See 

470 the documentation of the method for more information. 

471 

472 Parameters 

473 ---------- 

474 x : array_like, shape (M,) 

475 x-coordinates of the M sample points ``(x[i], y[i])``. 

476 y : array_like, shape (M,) or (M, K) 

477 y-coordinates of the sample points. Several data sets of sample 

478 points sharing the same x-coordinates can be fitted at once by 

479 passing in a 2D-array that contains one dataset per column. 

480 deg : int 

481 Degree of the fitting polynomial 

482 rcond : float, optional 

483 Relative condition number of the fit. Singular values smaller than 

484 this relative to the largest singular value will be ignored. The 

485 default value is len(x)*eps, where eps is the relative precision of 

486 the float type, about 2e-16 in most cases. 

487 full : bool, optional 

488 Switch determining nature of return value. When it is False (the 

489 default) just the coefficients are returned, when True diagnostic 

490 information from the singular value decomposition is also returned. 

491 w : array_like, shape (M,), optional 

492 Weights. If not None, the weight ``w[i]`` applies to the unsquared 

493 residual ``y[i] - y_hat[i]`` at ``x[i]``. Ideally the weights are 

494 chosen so that the errors of the products ``w[i]*y[i]`` all have the 

495 same variance. When using inverse-variance weighting, use 

496 ``w[i] = 1/sigma(y[i])``. The default value is None. 

497 cov : bool or str, optional 

498 If given and not `False`, return not just the estimate but also its 

499 covariance matrix. By default, the covariance are scaled by 

500 chi2/dof, where dof = M - (deg + 1), i.e., the weights are presumed 

501 to be unreliable except in a relative sense and everything is scaled 

502 such that the reduced chi2 is unity. This scaling is omitted if 

503 ``cov='unscaled'``, as is relevant for the case that the weights are 

504 w = 1/sigma, with sigma known to be a reliable estimate of the 

505 uncertainty. 

506 

507 Returns 

508 ------- 

509 p : ndarray, shape (deg + 1,) or (deg + 1, K) 

510 Polynomial coefficients, highest power first. If `y` was 2-D, the 

511 coefficients for `k`-th data set are in ``p[:,k]``. 

512 

513 residuals, rank, singular_values, rcond 

514 These values are only returned if ``full == True`` 

515 

516 - residuals -- sum of squared residuals of the least squares fit 

517 - rank -- the effective rank of the scaled Vandermonde 

518 coefficient matrix 

519 - singular_values -- singular values of the scaled Vandermonde 

520 coefficient matrix 

521 - rcond -- value of `rcond`. 

522 

523 For more details, see `numpy.linalg.lstsq`. 

524 

525 V : ndarray, shape (M,M) or (M,M,K) 

526 Present only if ``full == False`` and ``cov == True``. The covariance 

527 matrix of the polynomial coefficient estimates. The diagonal of 

528 this matrix are the variance estimates for each coefficient. If y 

529 is a 2-D array, then the covariance matrix for the `k`-th data set 

530 are in ``V[:,:,k]`` 

531 

532 

533 Warns 

534 ----- 

535 RankWarning 

536 The rank of the coefficient matrix in the least-squares fit is 

537 deficient. The warning is only raised if ``full == False``. 

538 

539 The warnings can be turned off by 

540 

541 >>> import warnings 

542 >>> warnings.simplefilter('ignore', np.RankWarning) 

543 

544 See Also 

545 -------- 

546 polyval : Compute polynomial values. 

547 linalg.lstsq : Computes a least-squares fit. 

548 scipy.interpolate.UnivariateSpline : Computes spline fits. 

549 

550 Notes 

551 ----- 

552 The solution minimizes the squared error 

553 

554 .. math:: 

555 E = \\sum_{j=0}^k |p(x_j) - y_j|^2 

556 

557 in the equations:: 

558 

559 x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0] 

560 x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1] 

561 ... 

562 x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k] 

563 

564 The coefficient matrix of the coefficients `p` is a Vandermonde matrix. 

565 

566 `polyfit` issues a `RankWarning` when the least-squares fit is badly 

567 conditioned. This implies that the best fit is not well-defined due 

568 to numerical error. The results may be improved by lowering the polynomial 

569 degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter 

570 can also be set to a value smaller than its default, but the resulting 

571 fit may be spurious: including contributions from the small singular 

572 values can add numerical noise to the result. 

573 

574 Note that fitting polynomial coefficients is inherently badly conditioned 

575 when the degree of the polynomial is large or the interval of sample points 

576 is badly centered. The quality of the fit should always be checked in these 

577 cases. When polynomial fits are not satisfactory, splines may be a good 

578 alternative. 

579 

580 References 

581 ---------- 

582 .. [1] Wikipedia, "Curve fitting", 

583 https://en.wikipedia.org/wiki/Curve_fitting 

584 .. [2] Wikipedia, "Polynomial interpolation", 

585 https://en.wikipedia.org/wiki/Polynomial_interpolation 

586 

587 Examples 

588 -------- 

589 >>> import warnings 

590 >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) 

591 >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0]) 

592 >>> z = np.polyfit(x, y, 3) 

593 >>> z 

594 array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) # may vary 

595 

596 It is convenient to use `poly1d` objects for dealing with polynomials: 

597 

598 >>> p = np.poly1d(z) 

599 >>> p(0.5) 

600 0.6143849206349179 # may vary 

601 >>> p(3.5) 

602 -0.34732142857143039 # may vary 

603 >>> p(10) 

604 22.579365079365115 # may vary 

605 

606 High-order polynomials may oscillate wildly: 

607 

608 >>> with warnings.catch_warnings(): 

609 ... warnings.simplefilter('ignore', np.RankWarning) 

610 ... p30 = np.poly1d(np.polyfit(x, y, 30)) 

611 ... 

612 >>> p30(4) 

613 -0.80000000000000204 # may vary 

614 >>> p30(5) 

615 -0.99999999999999445 # may vary 

616 >>> p30(4.5) 

617 -0.10547061179440398 # may vary 

618 

619 Illustration: 

620 

621 >>> import matplotlib.pyplot as plt 

622 >>> xp = np.linspace(-2, 6, 100) 

623 >>> _ = plt.plot(x, y, '.', xp, p(xp), '-', xp, p30(xp), '--') 

624 >>> plt.ylim(-2,2) 

625 (-2, 2) 

626 >>> plt.show() 

627 

628 """ 

629 order = int(deg) + 1 

630 x = NX.asarray(x) + 0.0 

631 y = NX.asarray(y) + 0.0 

632 

633 # check arguments. 

634 if deg < 0: 

635 raise ValueError("expected deg >= 0") 

636 if x.ndim != 1: 

637 raise TypeError("expected 1D vector for x") 

638 if x.size == 0: 

639 raise TypeError("expected non-empty vector for x") 

640 if y.ndim < 1 or y.ndim > 2: 

641 raise TypeError("expected 1D or 2D array for y") 

642 if x.shape[0] != y.shape[0]: 

643 raise TypeError("expected x and y to have same length") 

644 

645 # set rcond 

646 if rcond is None: 

647 rcond = len(x)*finfo(x.dtype).eps 

648 

649 # set up least squares equation for powers of x 

650 lhs = vander(x, order) 

651 rhs = y 

652 

653 # apply weighting 

654 if w is not None: 

655 w = NX.asarray(w) + 0.0 

656 if w.ndim != 1: 

657 raise TypeError("expected a 1-d array for weights") 

658 if w.shape[0] != y.shape[0]: 

659 raise TypeError("expected w and y to have the same length") 

660 lhs *= w[:, NX.newaxis] 

661 if rhs.ndim == 2: 

662 rhs *= w[:, NX.newaxis] 

663 else: 

664 rhs *= w 

665 

666 # scale lhs to improve condition number and solve 

667 scale = NX.sqrt((lhs*lhs).sum(axis=0)) 

668 lhs /= scale 

669 c, resids, rank, s = lstsq(lhs, rhs, rcond) 

670 c = (c.T/scale).T # broadcast scale coefficients 

671 

672 # warn on rank reduction, which indicates an ill conditioned matrix 

673 if rank != order and not full: 

674 msg = "Polyfit may be poorly conditioned" 

675 warnings.warn(msg, RankWarning, stacklevel=2) 

676 

677 if full: 

678 return c, resids, rank, s, rcond 

679 elif cov: 

680 Vbase = inv(dot(lhs.T, lhs)) 

681 Vbase /= NX.outer(scale, scale) 

682 if cov == "unscaled": 

683 fac = 1 

684 else: 

685 if len(x) <= order: 

686 raise ValueError("the number of data points must exceed order " 

687 "to scale the covariance matrix") 

688 # note, this used to be: fac = resids / (len(x) - order - 2.0) 

689 # it was deciced that the "- 2" (originally justified by "Bayesian 

690 # uncertainty analysis") is not what the user expects 

691 # (see gh-11196 and gh-11197) 

692 fac = resids / (len(x) - order) 

693 if y.ndim == 1: 

694 return c, Vbase * fac 

695 else: 

696 return c, Vbase[:,:, NX.newaxis] * fac 

697 else: 

698 return c 

699 

700 

701def _polyval_dispatcher(p, x): 

702 return (p, x) 

703 

704 

705@array_function_dispatch(_polyval_dispatcher) 

706def polyval(p, x): 

707 """ 

708 Evaluate a polynomial at specific values. 

709 

710 .. note:: 

711 This forms part of the old polynomial API. Since version 1.4, the 

712 new polynomial API defined in `numpy.polynomial` is preferred. 

713 A summary of the differences can be found in the 

714 :doc:`transition guide </reference/routines.polynomials>`. 

715 

716 If `p` is of length N, this function returns the value: 

717 

718 ``p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]`` 

719 

720 If `x` is a sequence, then ``p(x)`` is returned for each element of ``x``. 

721 If `x` is another polynomial then the composite polynomial ``p(x(t))`` 

722 is returned. 

723 

724 Parameters 

725 ---------- 

726 p : array_like or poly1d object 

727 1D array of polynomial coefficients (including coefficients equal 

728 to zero) from highest degree to the constant term, or an 

729 instance of poly1d. 

730 x : array_like or poly1d object 

731 A number, an array of numbers, or an instance of poly1d, at 

732 which to evaluate `p`. 

733 

734 Returns 

735 ------- 

736 values : ndarray or poly1d 

737 If `x` is a poly1d instance, the result is the composition of the two 

738 polynomials, i.e., `x` is "substituted" in `p` and the simplified 

739 result is returned. In addition, the type of `x` - array_like or 

740 poly1d - governs the type of the output: `x` array_like => `values` 

741 array_like, `x` a poly1d object => `values` is also. 

742 

743 See Also 

744 -------- 

745 poly1d: A polynomial class. 

746 

747 Notes 

748 ----- 

749 Horner's scheme [1]_ is used to evaluate the polynomial. Even so, 

750 for polynomials of high degree the values may be inaccurate due to 

751 rounding errors. Use carefully. 

752 

753 If `x` is a subtype of `ndarray` the return value will be of the same type. 

754 

755 References 

756 ---------- 

757 .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng. 

758 trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand 

759 Reinhold Co., 1985, pg. 720. 

760 

761 Examples 

762 -------- 

763 >>> np.polyval([3,0,1], 5) # 3 * 5**2 + 0 * 5**1 + 1 

764 76 

765 >>> np.polyval([3,0,1], np.poly1d(5)) 

766 poly1d([76]) 

767 >>> np.polyval(np.poly1d([3,0,1]), 5) 

768 76 

769 >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5)) 

770 poly1d([76]) 

771 

772 """ 

773 p = NX.asarray(p) 

774 if isinstance(x, poly1d): 

775 y = 0 

776 else: 

777 x = NX.asanyarray(x) 

778 y = NX.zeros_like(x) 

779 for pv in p: 

780 y = y * x + pv 

781 return y 

782 

783 

784def _binary_op_dispatcher(a1, a2): 

785 return (a1, a2) 

786 

787 

788@array_function_dispatch(_binary_op_dispatcher) 

789def polyadd(a1, a2): 

790 """ 

791 Find the sum of two polynomials. 

792 

793 .. note:: 

794 This forms part of the old polynomial API. Since version 1.4, the 

795 new polynomial API defined in `numpy.polynomial` is preferred. 

796 A summary of the differences can be found in the 

797 :doc:`transition guide </reference/routines.polynomials>`. 

798 

799 Returns the polynomial resulting from the sum of two input polynomials. 

800 Each input must be either a poly1d object or a 1D sequence of polynomial 

801 coefficients, from highest to lowest degree. 

802 

803 Parameters 

804 ---------- 

805 a1, a2 : array_like or poly1d object 

806 Input polynomials. 

807 

808 Returns 

809 ------- 

810 out : ndarray or poly1d object 

811 The sum of the inputs. If either input is a poly1d object, then the 

812 output is also a poly1d object. Otherwise, it is a 1D array of 

813 polynomial coefficients from highest to lowest degree. 

814 

815 See Also 

816 -------- 

817 poly1d : A one-dimensional polynomial class. 

818 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval 

819 

820 Examples 

821 -------- 

822 >>> np.polyadd([1, 2], [9, 5, 4]) 

823 array([9, 6, 6]) 

824 

825 Using poly1d objects: 

826 

827 >>> p1 = np.poly1d([1, 2]) 

828 >>> p2 = np.poly1d([9, 5, 4]) 

829 >>> print(p1) 

830 1 x + 2 

831 >>> print(p2) 

832 2 

833 9 x + 5 x + 4 

834 >>> print(np.polyadd(p1, p2)) 

835 2 

836 9 x + 6 x + 6 

837 

838 """ 

839 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d)) 

840 a1 = atleast_1d(a1) 

841 a2 = atleast_1d(a2) 

842 diff = len(a2) - len(a1) 

843 if diff == 0: 

844 val = a1 + a2 

845 elif diff > 0: 

846 zr = NX.zeros(diff, a1.dtype) 

847 val = NX.concatenate((zr, a1)) + a2 

848 else: 

849 zr = NX.zeros(abs(diff), a2.dtype) 

850 val = a1 + NX.concatenate((zr, a2)) 

851 if truepoly: 

852 val = poly1d(val) 

853 return val 

854 

855 

856@array_function_dispatch(_binary_op_dispatcher) 

857def polysub(a1, a2): 

858 """ 

859 Difference (subtraction) of two polynomials. 

860 

861 .. note:: 

862 This forms part of the old polynomial API. Since version 1.4, the 

863 new polynomial API defined in `numpy.polynomial` is preferred. 

864 A summary of the differences can be found in the 

865 :doc:`transition guide </reference/routines.polynomials>`. 

866 

867 Given two polynomials `a1` and `a2`, returns ``a1 - a2``. 

868 `a1` and `a2` can be either array_like sequences of the polynomials' 

869 coefficients (including coefficients equal to zero), or `poly1d` objects. 

870 

871 Parameters 

872 ---------- 

873 a1, a2 : array_like or poly1d 

874 Minuend and subtrahend polynomials, respectively. 

875 

876 Returns 

877 ------- 

878 out : ndarray or poly1d 

879 Array or `poly1d` object of the difference polynomial's coefficients. 

880 

881 See Also 

882 -------- 

883 polyval, polydiv, polymul, polyadd 

884 

885 Examples 

886 -------- 

887 .. math:: (2 x^2 + 10 x - 2) - (3 x^2 + 10 x -4) = (-x^2 + 2) 

888 

889 >>> np.polysub([2, 10, -2], [3, 10, -4]) 

890 array([-1, 0, 2]) 

891 

892 """ 

893 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d)) 

894 a1 = atleast_1d(a1) 

895 a2 = atleast_1d(a2) 

896 diff = len(a2) - len(a1) 

897 if diff == 0: 

898 val = a1 - a2 

899 elif diff > 0: 

900 zr = NX.zeros(diff, a1.dtype) 

901 val = NX.concatenate((zr, a1)) - a2 

902 else: 

903 zr = NX.zeros(abs(diff), a2.dtype) 

904 val = a1 - NX.concatenate((zr, a2)) 

905 if truepoly: 

906 val = poly1d(val) 

907 return val 

908 

909 

910@array_function_dispatch(_binary_op_dispatcher) 

911def polymul(a1, a2): 

912 """ 

913 Find the product of two polynomials. 

914 

915 .. note:: 

916 This forms part of the old polynomial API. Since version 1.4, the 

917 new polynomial API defined in `numpy.polynomial` is preferred. 

918 A summary of the differences can be found in the 

919 :doc:`transition guide </reference/routines.polynomials>`. 

920 

921 Finds the polynomial resulting from the multiplication of the two input 

922 polynomials. Each input must be either a poly1d object or a 1D sequence 

923 of polynomial coefficients, from highest to lowest degree. 

924 

925 Parameters 

926 ---------- 

927 a1, a2 : array_like or poly1d object 

928 Input polynomials. 

929 

930 Returns 

931 ------- 

932 out : ndarray or poly1d object 

933 The polynomial resulting from the multiplication of the inputs. If 

934 either inputs is a poly1d object, then the output is also a poly1d 

935 object. Otherwise, it is a 1D array of polynomial coefficients from 

936 highest to lowest degree. 

937 

938 See Also 

939 -------- 

940 poly1d : A one-dimensional polynomial class. 

941 poly, polyadd, polyder, polydiv, polyfit, polyint, polysub, polyval 

942 convolve : Array convolution. Same output as polymul, but has parameter 

943 for overlap mode. 

944 

945 Examples 

946 -------- 

947 >>> np.polymul([1, 2, 3], [9, 5, 1]) 

948 array([ 9, 23, 38, 17, 3]) 

949 

950 Using poly1d objects: 

951 

952 >>> p1 = np.poly1d([1, 2, 3]) 

953 >>> p2 = np.poly1d([9, 5, 1]) 

954 >>> print(p1) 

955 2 

956 1 x + 2 x + 3 

957 >>> print(p2) 

958 2 

959 9 x + 5 x + 1 

960 >>> print(np.polymul(p1, p2)) 

961 4 3 2 

962 9 x + 23 x + 38 x + 17 x + 3 

963 

964 """ 

965 truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d)) 

966 a1, a2 = poly1d(a1), poly1d(a2) 

967 val = NX.convolve(a1, a2) 

968 if truepoly: 

969 val = poly1d(val) 

970 return val 

971 

972 

973def _polydiv_dispatcher(u, v): 

974 return (u, v) 

975 

976 

977@array_function_dispatch(_polydiv_dispatcher) 

978def polydiv(u, v): 

979 """ 

980 Returns the quotient and remainder of polynomial division. 

981 

982 .. note:: 

983 This forms part of the old polynomial API. Since version 1.4, the 

984 new polynomial API defined in `numpy.polynomial` is preferred. 

985 A summary of the differences can be found in the 

986 :doc:`transition guide </reference/routines.polynomials>`. 

987 

988 The input arrays are the coefficients (including any coefficients 

989 equal to zero) of the "numerator" (dividend) and "denominator" 

990 (divisor) polynomials, respectively. 

991 

992 Parameters 

993 ---------- 

994 u : array_like or poly1d 

995 Dividend polynomial's coefficients. 

996 

997 v : array_like or poly1d 

998 Divisor polynomial's coefficients. 

999 

1000 Returns 

1001 ------- 

1002 q : ndarray 

1003 Coefficients, including those equal to zero, of the quotient. 

1004 r : ndarray 

1005 Coefficients, including those equal to zero, of the remainder. 

1006 

1007 See Also 

1008 -------- 

1009 poly, polyadd, polyder, polydiv, polyfit, polyint, polymul, polysub 

1010 polyval 

1011 

1012 Notes 

1013 ----- 

1014 Both `u` and `v` must be 0-d or 1-d (ndim = 0 or 1), but `u.ndim` need 

1015 not equal `v.ndim`. In other words, all four possible combinations - 

1016 ``u.ndim = v.ndim = 0``, ``u.ndim = v.ndim = 1``, 

1017 ``u.ndim = 1, v.ndim = 0``, and ``u.ndim = 0, v.ndim = 1`` - work. 

1018 

1019 Examples 

1020 -------- 

1021 .. math:: \\frac{3x^2 + 5x + 2}{2x + 1} = 1.5x + 1.75, remainder 0.25 

1022 

1023 >>> x = np.array([3.0, 5.0, 2.0]) 

1024 >>> y = np.array([2.0, 1.0]) 

1025 >>> np.polydiv(x, y) 

1026 (array([1.5 , 1.75]), array([0.25])) 

1027 

1028 """ 

1029 truepoly = (isinstance(u, poly1d) or isinstance(v, poly1d)) 

1030 u = atleast_1d(u) + 0.0 

1031 v = atleast_1d(v) + 0.0 

1032 # w has the common type 

1033 w = u[0] + v[0] 

1034 m = len(u) - 1 

1035 n = len(v) - 1 

1036 scale = 1. / v[0] 

1037 q = NX.zeros((max(m - n + 1, 1),), w.dtype) 

1038 r = u.astype(w.dtype) 

1039 for k in range(0, m-n+1): 

1040 d = scale * r[k] 

1041 q[k] = d 

1042 r[k:k+n+1] -= d*v 

1043 while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1): 

1044 r = r[1:] 

1045 if truepoly: 

1046 return poly1d(q), poly1d(r) 

1047 return q, r 

1048 

1049_poly_mat = re.compile(r"\*\*([0-9]*)") 

1050def _raise_power(astr, wrap=70): 

1051 n = 0 

1052 line1 = '' 

1053 line2 = '' 

1054 output = ' ' 

1055 while True: 

1056 mat = _poly_mat.search(astr, n) 

1057 if mat is None: 

1058 break 

1059 span = mat.span() 

1060 power = mat.groups()[0] 

1061 partstr = astr[n:span[0]]