Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-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 

12import numpy.core.numeric as NX 

13 

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 

22 

23 

24array_function_dispatch = functools.partial( 

25 overrides.array_function_dispatch, module='numpy') 

26 

27 

28@set_module('numpy') 

29class RankWarning(UserWarning): 

30 """ 

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

32 

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

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

35 

36 """ 

37 pass 

38 

39 

40def _poly_dispatcher(seq_of_zeros): 

41 return seq_of_zeros 

42 

43 

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. 

48 

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>`. 

54 

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. 

61 

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. 

66 

67 Returns 

68 ------- 

69 c : ndarray 

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

71 

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

73 where c[0] always equals 1. 

74 

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). 

80 

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. 

87 

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``.) 

96 

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

98 matrix **A** is given by 

99 

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

101 

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

103 

104 References 

105 ---------- 

106 .. [1] M. Sullivan and M. Sullivan, III, "Algebra and Trignometry, 

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

108 

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

110 Academic Press, pg. 182, 1980. 

111 

112 Examples 

113 -------- 

114 Given a sequence of a polynomial's zeros: 

115 

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

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

118 

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

120 

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

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

123 

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

125 

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

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

128 

129 Given a square array object: 

130 

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

132 >>> np.poly(P) 

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

134 

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

136 

137 """ 

138 seq_of_zeros = atleast_1d(seq_of_zeros) 

139 sh = seq_of_zeros.shape 

140 

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.") 

150 

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') 

157 

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() 

163 

164 return a 

165 

166 

167def _roots_dispatcher(p): 

168 return p 

169 

170 

171@array_function_dispatch(_roots_dispatcher) 

172def roots(p): 

173 """ 

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

175 

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>`. 

181 

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:: 

184 

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

186 

187 Parameters 

188 ---------- 

189 p : array_like 

190 Rank-1 array of polynomial coefficients. 

191 

192 Returns 

193 ------- 

194 out : ndarray 

195 An array containing the roots of the polynomial. 

196 

197 Raises 

198 ------ 

199 ValueError 

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

201 

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. 

209 

210 Notes 

211 ----- 

212 The algorithm relies on computing the eigenvalues of the 

213 companion matrix [1]_. 

214 

215 References 

216 ---------- 

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

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

219 

220 Examples 

221 -------- 

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 The defining property of the antiderivative: 

309 

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 

316 

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

318 

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 

329 

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: 

332 

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

334 6.0 

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

336 5.0 

337 >>> P(0) 

338 3.0 

339 

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.") 

352 

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 

366 

367 

368def _polyder_dispatcher(p, m=None): 

369 return (p,) 

370 

371 

372@array_function_dispatch(_polyder_dispatcher) 

373def polyder(p, m=1): 

374 """ 

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

376 

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>`. 

382 

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) 

390 

391 Returns 

392 ------- 

393 der : poly1d 

394 A new polynomial representing the derivative. 

395 

396 See Also 

397 -------- 

398 polyint : Anti-derivative of a polynomial. 

399 poly1d : Class for one-dimensional polynomials. 

400 

401 Examples 

402 -------- 

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

404 

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

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

407 >>> p2 

408 poly1d([3, 2, 1]) 

409 

410 which evaluates to: 

411 

412 >>> p2(2.) 

413 17.0 

414 

415 We can verify this, approximating the derivative with 

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

417 

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

419 17.007000999997857 

420 

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

422 

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]) 

429 

430 """ 

431 m = int(m) 

432 if m < 0: 

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

434 

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 

446 

447 

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

449 return (x, y, w) 

450 

451 

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. 

456 

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>`. 

462 

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`. 

466 

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. 

470 

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. 

505 

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]``. 

511 

512 residuals, rank, singular_values, rcond 

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

514 

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`. 

521 

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

523 

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]`` 

530 

531 

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``. 

537 

538 The warnings can be turned off by 

539 

540 >>> import warnings 

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

542 

543 See Also 

544 -------- 

545 polyval : Compute polynomial values. 

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

547 scipy.interpolate.UnivariateSpline : Computes spline fits. 

548 

549 Notes 

550 ----- 

551 The solution minimizes the squared error 

552 

553 .. math:: 

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

555 

556 in the equations:: 

557 

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] 

562 

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

564 

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. 

572 

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. 

578 

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 

585 

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 

594 

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

596 

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 

604 

605 High-order polynomials may oscillate wildly: 

606 

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 

617 

618 Illustration: 

619 

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() 

626 

627 """ 

628 order = int(deg) + 1 

629 x = NX.asarray(x) + 0.0 

630 y = NX.asarray(y) + 0.0 

631 

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") 

643 

644 # set rcond 

645 if rcond is None: 

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

647 

648 # set up least squares equation for powers of x 

649 lhs = vander(x, order) 

650 rhs = y 

651 

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 

664 

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 

670 

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) 

675 

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 

698 

699 

700def _polyval_dispatcher(p, x): 

701 return (p, x) 

702 

703 

704@array_function_dispatch(_polyval_dispatcher) 

705def polyval(p, x): 

706 """ 

707 Evaluate a polynomial at specific values. 

708 

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>`. 

714 

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

716 

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

718 

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. 

722 

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`. 

732 

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. 

741 

742 See Also 

743 -------- 

744 poly1d: A polynomial class. 

745 

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. 

751 

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

753 

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. 

759 

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]) 

770 

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 

781 

782 

783def _binary_op_dispatcher(a1, a2): 

784 return (a1, a2) 

785 

786 

787@array_function_dispatch(_binary_op_dispatcher) 

788def polyadd(a1, a2): 

789 """ 

790 Find the sum of two polynomials. 

791 

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>`. 

797 

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. 

801 

802 Parameters 

803 ---------- 

804 a1, a2 : array_like or poly1d object 

805 Input polynomials. 

806 

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. 

813 

814 See Also 

815 -------- 

816 poly1d : A one-dimensional polynomial class. 

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

818 

819 Examples 

820 -------- 

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

822 array([9, 6, 6]) 

823 

824 Using poly1d objects: 

825 

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 

836 

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 

853 

854 

855@array_function_dispatch(_binary_op_dispatcher) 

856def polysub(a1, a2): 

857 """ 

858 Difference (subtraction) of two polynomials. 

859 

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>`. 

865 

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. 

869 

870 Parameters 

871 ---------- 

872 a1, a2 : array_like or poly1d 

873 Minuend and subtrahend polynomials, respectively. 

874 

875 Returns 

876 ------- 

877 out : ndarray or poly1d 

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

879 

880 See Also 

881 -------- 

882 polyval, polydiv, polymul, polyadd 

883 

884 Examples 

885 -------- 

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

887 

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

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

890 

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 

907 

908 

909@array_function_dispatch(_binary_op_dispatcher) 

910def polymul(a1, a2): 

911 """ 

912 Find the product of two polynomials. 

913 

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>`. 

919 

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. 

923 

924 Parameters 

925 ---------- 

926 a1, a2 : array_like or poly1d object 

927 Input polynomials. 

928 

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. 

936 

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. 

943 

944 Examples 

945 -------- 

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

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

948 

949 Using poly1d objects: 

950 

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 

962 

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 

970 

971 

972def _polydiv_dispatcher(u, v): 

973 return (u, v) 

974 

975 

976@array_function_dispatch(_polydiv_dispatcher) 

977def polydiv(u, v): 

978 """ 

979 Returns the quotient and remainder of polynomial division. 

980 

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>`. 

986 

987 The input arrays are the coefficients (including any coefficients 

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

989 (divisor) polynomials, respectively. 

990 

991 Parameters 

992 ---------- 

993 u : array_like or poly1d 

994 Dividend polynomial's coefficients. 

995 

996 v : array_like or poly1d 

997 Divisor polynomial's coefficients. 

998 

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. 

1005 

1006 See Also 

1007 -------- 

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

1009 polyval 

1010 

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. 

1017 

1018 Examples 

1019 -------- 

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

1021 

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])) 

1026 

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 

1047 

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]]