1""" 2================================================= 3Power Series (:mod:`numpy.polynomial.polynomial`) 4================================================= 5 6This module provides a number of objects (mostly functions) useful for 7dealing with polynomials, including a `Polynomial` class that 8encapsulates the usual arithmetic operations. (General information 9on how this module represents and works with polynomial objects is in 10the docstring for its "parent" sub-package, `numpy.polynomial`). 11 12Classes 13------- 14.. autosummary:: 15 :toctree: generated/ 16 17 Polynomial 18 19Constants 20--------- 21.. autosummary:: 22 :toctree: generated/ 23 24 polydomain 25 polyzero 26 polyone 27 polyx 28 29Arithmetic 30---------- 31.. autosummary:: 32 :toctree: generated/ 33 34 polyadd 35 polysub 36 polymulx 37 polymul 38 polydiv 39 polypow 40 polyval 41 polyval2d 42 polyval3d 43 polygrid2d 44 polygrid3d 45 46Calculus 47-------- 48.. autosummary:: 49 :toctree: generated/ 50 51 polyder 52 polyint 53 54Misc Functions 55-------------- 56.. autosummary:: 57 :toctree: generated/ 58 59 polyfromroots 60 polyroots 61 polyvalfromroots 62 polyvander 63 polyvander2d 64 polyvander3d 65 polycompanion 66 polyfit 67 polytrim 68 polyline 69 70See Also 71-------- 72`numpy.polynomial` 73 74""" 75__all__ = [ 76 'polyzero', 'polyone', 'polyx', 'polydomain', 'polyline', 'polyadd', 77 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow', 'polyval', 78 'polyvalfromroots', 'polyder', 'polyint', 'polyfromroots', 'polyvander', 79 'polyfit', 'polytrim', 'polyroots', 'Polynomial', 'polyval2d', 'polyval3d', 80 'polygrid2d', 'polygrid3d', 'polyvander2d', 'polyvander3d'] 81 82import numpy as np 83import numpy.linalg as la 84from numpy.core.multiarray import normalize_axis_index 85 86from . import polyutils as pu 87from ._polybase import ABCPolyBase 88 89polytrim = pu.trimcoef 90 91# 92# These are constant arrays are of integer type so as to be compatible 93# with the widest range of other types, such as Decimal. 94# 95 96# Polynomial default domain. 97polydomain = np.array([-1, 1]) 98 99# Polynomial coefficients representing zero. 100polyzero = np.array([0]) 101 102# Polynomial coefficients representing one. 103polyone = np.array([1]) 104 105# Polynomial coefficients representing the identity x. 106polyx = np.array([0, 1]) 107 108# 109# Polynomial series functions 110# 111 112 113def polyline(off, scl): 114 """ 115 Returns an array representing a linear polynomial. 116 117 Parameters 118 ---------- 119 off, scl : scalars 120 The "y-intercept" and "slope" of the line, respectively. 121 122 Returns 123 ------- 124 y : ndarray 125 This module's representation of the linear polynomial ``off + 126 scl*x``. 127 128 See Also 129 -------- 130 numpy.polynomial.chebyshev.chebline 131 numpy.polynomial.legendre.legline 132 numpy.polynomial.laguerre.lagline 133 numpy.polynomial.hermite.hermline 134 numpy.polynomial.hermite_e.hermeline 135 136 Examples 137 -------- 138 >>> from numpy.polynomial import polynomial as P 139 >>> P.polyline(1,-1) 140 array([ 1, -1]) 141 >>> P.polyval(1, P.polyline(1,-1)) # should be 0 142 0.0 143 144 """ 145 if scl != 0: 146 return np.array([off, scl]) 147 else: 148 return np.array([off]) 149 150 151def polyfromroots(roots): 152 """ 153 Generate a monic polynomial with given roots. 154 155 Return the coefficients of the polynomial 156 157 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), 158 159 where the `r_n` are the roots specified in `roots`. If a zero has 160 multiplicity n, then it must appear in `roots` n times. For instance, 161 if 2 is a root of multiplicity three and 3 is a root of multiplicity 2, 162 then `roots` looks something like [2, 2, 2, 3, 3]. The roots can appear 163 in any order. 164 165 If the returned coefficients are `c`, then 166 167 .. math:: p(x) = c_0 + c_1 * x + ... + x^n 168 169 The coefficient of the last term is 1 for monic polynomials in this 170 form. 171 172 Parameters 173 ---------- 174 roots : array_like 175 Sequence containing the roots. 176 177 Returns 178 ------- 179 out : ndarray 180 1-D array of the polynomial's coefficients If all the roots are 181 real, then `out` is also real, otherwise it is complex. (see 182 Examples below). 183 184 See Also 185 -------- 186 numpy.polynomial.chebyshev.chebfromroots 187 numpy.polynomial.legendre.legfromroots 188 numpy.polynomial.laguerre.lagfromroots 189 numpy.polynomial.hermite.hermfromroots 190 numpy.polynomial.hermite_e.hermefromroots 191 192 Notes 193 ----- 194 The coefficients are determined by multiplying together linear factors 195 of the form `(x - r_i)`, i.e. 196 197 .. math:: p(x) = (x - r_0) (x - r_1) ... (x - r_n) 198 199 where ``n == len(roots) - 1``; note that this implies that `1` is always 200 returned for :math:`a_n`. 201 202 Examples 203 -------- 204 >>> from numpy.polynomial import polynomial as P 205 >>> P.polyfromroots((-1,0,1)) # x(x - 1)(x + 1) = x^3 - x 206 array([ 0., -1., 0., 1.]) 207 >>> j = complex(0,1) 208 >>> P.polyfromroots((-j,j)) # complex returned, though values are real 209 array([1.+0.j, 0.+0.j, 1.+0.j]) 210 211 """ 212 return pu._fromroots(polyline, polymul, roots) 213 214 215def polyadd(c1, c2): 216 """ 217 Add one polynomial to another. 218 219 Returns the sum of two polynomials `c1` + `c2`. The arguments are 220 sequences of coefficients from lowest order term to highest, i.e., 221 [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``. 222 223 Parameters 224 ---------- 225 c1, c2 : array_like 226 1-D arrays of polynomial coefficients ordered from low to high. 227 228 Returns 229 ------- 230 out : ndarray 231 The coefficient array representing their sum. 232 233 See Also 234 -------- 235 polysub, polymulx, polymul, polydiv, polypow 236 237 Examples 238 -------- 239 >>> from numpy.polynomial import polynomial as P 240 >>> c1 = (1,2,3) 241 >>> c2 = (3,2,1) 242 >>> sum = P.polyadd(c1,c2); sum 243 array([4., 4., 4.]) 244 >>> P.polyval(2, sum) # 4 + 4(2) + 4(2**2) 245 28.0 246 247 """ 248 return pu._add(c1, c2) 249 250 251def polysub(c1, c2): 252 """ 253 Subtract one polynomial from another. 254 255 Returns the difference of two polynomials `c1` - `c2`. The arguments 256 are sequences of coefficients from lowest order term to highest, i.e., 257 [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``. 258 259 Parameters 260 ---------- 261 c1, c2 : array_like 262 1-D arrays of polynomial coefficients ordered from low to 263 high. 264 265 Returns 266 ------- 267 out : ndarray 268 Of coefficients representing their difference. 269 270 See Also 271 -------- 272 polyadd, polymulx, polymul, polydiv, polypow 273 274 Examples 275 -------- 276 >>> from numpy.polynomial import polynomial as P 277 >>> c1 = (1,2,3) 278 >>> c2 = (3,2,1) 279 >>> P.polysub(c1,c2) 280 array([-2., 0., 2.]) 281 >>> P.polysub(c2,c1) # -P.polysub(c1,c2) 282 array([ 2., 0., -2.]) 283 284 """ 285 return pu._sub(c1, c2) 286 287 288def polymulx(c): 289 """Multiply a polynomial by x. 290 291 Multiply the polynomial `c` by x, where x is the independent 292 variable. 293 294 295 Parameters 296 ---------- 297 c : array_like 298 1-D array of polynomial coefficients ordered from low to 299 high. 300 301 Returns 302 ------- 303 out : ndarray 304 Array representing the result of the multiplication. 305 306 See Also 307 -------- 308 polyadd, polysub, polymul, polydiv, polypow 309 310 Notes 311 ----- 312 313 .. versionadded:: 1.5.0 314 315 """ 316 # c is a trimmed copy 317 [c] = pu.as_series([c]) 318 # The zero series needs special treatment 319 if len(c) == 1 and c[0] == 0: 320 return c 321 322 prd = np.empty(len(c) + 1, dtype=c.dtype) 323 prd[0] = c[0]*0 324 prd[1:] = c 325 return prd 326 327 328def polymul(c1, c2): 329 """ 330 Multiply one polynomial by another. 331 332 Returns the product of two polynomials `c1` * `c2`. The arguments are 333 sequences of coefficients, from lowest order term to highest, e.g., 334 [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2.`` 335 336 Parameters 337 ---------- 338 c1, c2 : array_like 339 1-D arrays of coefficients representing a polynomial, relative to the 340 "standard" basis, and ordered from lowest order term to highest. 341 342 Returns 343 ------- 344 out : ndarray 345 Of the coefficients of their product. 346 347 See Also 348 -------- 349 polyadd, polysub, polymulx, polydiv, polypow 350 351 Examples 352 -------- 353 >>> from numpy.polynomial import polynomial as P 354 >>> c1 = (1,2,3) 355 >>> c2 = (3,2,1) 356 >>> P.polymul(c1,c2) 357 array([ 3., 8., 14., 8., 3.]) 358 359 """ 360 # c1, c2 are trimmed copies 361 [c1, c2] = pu.as_series([c1, c2]) 362 ret = np.convolve(c1, c2) 363 return pu.trimseq(ret) 364 365 366def polydiv(c1, c2): 367 """ 368 Divide one polynomial by another. 369 370 Returns the quotient-with-remainder of two polynomials `c1` / `c2`. 371 The arguments are sequences of coefficients, from lowest order term 372 to highest, e.g., [1,2,3] represents ``1 + 2*x + 3*x**2``. 373 374 Parameters 375 ---------- 376 c1, c2 : array_like 377 1-D arrays of polynomial coefficients ordered from low to high. 378 379 Returns 380 ------- 381 [quo, rem] : ndarrays 382 Of coefficient series representing the quotient and remainder. 383 384 See Also 385 -------- 386 polyadd, polysub, polymulx, polymul, polypow 387 388 Examples 389 -------- 390 >>> from numpy.polynomial import polynomial as P 391 >>> c1 = (1,2,3) 392 >>> c2 = (3,2,1) 393 >>> P.polydiv(c1,c2) 394 (array([3.]), array([-8., -4.])) 395 >>> P.polydiv(c2,c1) 396 (array([ 0.33333333]), array([ 2.66666667, 1.33333333])) # may vary 397 398 """ 399 # c1, c2 are trimmed copies 400 [c1, c2] = pu.as_series([c1, c2]) 401 if c2[-1] == 0: 402 raise ZeroDivisionError() 403 404 # note: this is more efficient than `pu._div(polymul, c1, c2)` 405 lc1 = len(c1) 406 lc2 = len(c2) 407 if lc1 < lc2: 408 return c1[:1]*0, c1 409 elif lc2 == 1: 410 return c1/c2[-1], c1[:1]*0 411 else: 412 dlen = lc1 - lc2 413 scl = c2[-1] 414 c2 = c2[:-1]/scl 415 i = dlen 416 j = lc1 - 1 417 while i >= 0: 418 c1[i:j] -= c2*c1[j] 419 i -= 1 420 j -= 1 421 return c1[j+1:]/scl, pu.trimseq(c1[:j+1]) 422 423 424def polypow(c, pow, maxpower=None): 425 """Raise a polynomial to a power. 426 427 Returns the polynomial `c` raised to the power `pow`. The argument 428 `c` is a sequence of coefficients ordered from low to high. i.e., 429 [1,2,3] is the series ``1 + 2*x + 3*x**2.`` 430 431 Parameters 432 ---------- 433 c : array_like 434 1-D array of array of series coefficients ordered from low to 435 high degree. 436 pow : integer 437 Power to which the series will be raised 438 maxpower : integer, optional 439 Maximum power allowed. This is mainly to limit growth of the series 440 to unmanageable size. Default is 16 441 442 Returns 443 ------- 444 coef : ndarray 445 Power series of power. 446 447 See Also 448 -------- 449 polyadd, polysub, polymulx, polymul, polydiv 450 451 Examples 452 -------- 453 >>> from numpy.polynomial import polynomial as P 454 >>> P.polypow([1,2,3], 2) 455 array([ 1., 4., 10., 12., 9.]) 456 457 """ 458 # note: this is more efficient than `pu._pow(polymul, c1, c2)`, as it 459 # avoids calling `as_series` repeatedly 460 return pu._pow(np.convolve, c, pow, maxpower) 461 462 463def polyder(c, m=1, scl=1, axis=0): 464 """ 465 Differentiate a polynomial. 466 467 Returns the polynomial coefficients `c` differentiated `m` times along 468 `axis`. At each iteration the result is multiplied by `scl` (the 469 scaling factor is for use in a linear change of variable). The 470 argument `c` is an array of coefficients from low to high degree along 471 each axis, e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2`` 472 while [[1,2],[1,2]] represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is 473 ``x`` and axis=1 is ``y``. 474 475 Parameters 476 ---------- 477 c : array_like 478 Array of polynomial coefficients. If c is multidimensional the 479 different axis correspond to different variables with the degree 480 in each axis given by the corresponding index. 481 m : int, optional 482 Number of derivatives taken, must be non-negative. (Default: 1) 483 scl : scalar, optional 484 Each differentiation is multiplied by `scl`. The end result is 485 multiplication by ``scl**m``. This is for use in a linear change 486 of variable. (Default: 1) 487 axis : int, optional 488 Axis over which the derivative is taken. (Default: 0). 489 490 .. versionadded:: 1.7.0 491 492 Returns 493 ------- 494 der : ndarray 495 Polynomial coefficients of the derivative. 496 497 See Also 498 -------- 499 polyint 500 501 Examples 502 -------- 503 >>> from numpy.polynomial import polynomial as P 504 >>> c = (1,2,3,4) # 1 + 2x + 3x**2 + 4x**3 505 >>> P.polyder(c) # (d/dx)(c) = 2 + 6x + 12x**2 506 array([ 2., 6., 12.]) 507 >>> P.polyder(c,3) # (d**3/dx**3)(c) = 24 508 array([24.]) 509 >>> P.polyder(c,scl=-1) # (d/d(-x))(c) = -2 - 6x - 12x**2 510 array([ -2., -6., -12.]) 511 >>> P.polyder(c,2,-1) # (d**2/d(-x)**2)(c) = 6 + 24x 512 array([ 6., 24.]) 513 514 """ 515 c = np.array(c, ndmin=1, copy=True) 516 if c.dtype.char in '?bBhHiIlLqQpP': 517 # astype fails with NA 518 c = c + 0.0 519 cdt = c.dtype 520 cnt = pu._deprecate_as_int(m, "the order of derivation") 521 iaxis = pu._deprecate_as_int(axis, "the axis") 522 if cnt < 0: 523 raise ValueError("The order of derivation must be non-negative") 524 iaxis = normalize_axis_index(iaxis, c.ndim) 525 526 if cnt == 0: 527 return c 528 529 c = np.moveaxis(c, iaxis, 0) 530 n = len(c) 531 if cnt >= n: 532 c = c[:1]*0 533 else: 534 for i in range(cnt): 535 n = n - 1 536 c *= scl 537 der = np.empty((n,) + c.shape[1:], dtype=cdt) 538 for j in range(n, 0, -1): 539 der[j - 1] = j*c[j] 540 c = der 541 c = np.moveaxis(c, 0, iaxis) 542 return c 543 544 545def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0): 546 """ 547 Integrate a polynomial. 548 549 Returns the polynomial coefficients `c` integrated `m` times from 550 `lbnd` along `axis`. At each iteration the resulting series is 551 **multiplied** by `scl` and an integration constant, `k`, is added. 552 The scaling factor is for use in a linear change of variable. ("Buyer 553 beware": note that, depending on what one is doing, one may want `scl` 554 to be the reciprocal of what one might expect; for more information, 555 see the Notes section below.) The argument `c` is an array of 556 coefficients, from low to high degree along each axis, e.g., [1,2,3] 557 represents the polynomial ``1 + 2*x + 3*x**2`` while [[1,2],[1,2]] 558 represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is ``x`` and axis=1 is 559 ``y``. 560 561 Parameters 562 ---------- 563 c : array_like 564 1-D array of polynomial coefficients, ordered from low to high. 565 m : int, optional 566 Order of integration, must be positive. (Default: 1) 567 k : {[], list, scalar}, optional 568 Integration constant(s). The value of the first integral at zero 569 is the first value in the list, the value of the second integral 570 at zero is the second value, etc. If ``k == []`` (the default), 571 all constants are set to zero. If ``m == 1``, a single scalar can 572 be given instead of a list. 573 lbnd : scalar, optional 574 The lower bound of the integral. (Default: 0) 575 scl : scalar, optional 576 Following each integration the result is *multiplied* by `scl` 577 before the integration constant is added. (Default: 1) 578 axis : int, optional 579 Axis over which the integral is taken. (Default: 0). 580 581 .. versionadded:: 1.7.0 582 583 Returns 584 ------- 585 S : ndarray 586 Coefficient array of the integral. 587 588 Raises 589 ------ 590 ValueError 591 If ``m < 1``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or 592 ``np.ndim(scl) != 0``. 593 594 See Also 595 -------- 596 polyder 597 598 Notes 599 ----- 600 Note that the result of each integration is *multiplied* by `scl`. Why 601 is this important to note? Say one is making a linear change of 602 variable :math:`u = ax + b` in an integral relative to `x`. Then 603 :math:`dx = du/a`, so one will need to set `scl` equal to 604 :math:`1/a` - perhaps not what one would have first thought. 605 606 Examples 607 -------- 608 >>> from numpy.polynomial import polynomial as P 609 >>> c = (1,2,3) 610 >>> P.polyint(c) # should return array([0, 1, 1, 1]) 611 array([0., 1., 1., 1.]) 612 >>> P.polyint(c,3) # should return array([0, 0, 0, 1/6, 1/12, 1/20]) 613 array([ 0. , 0. , 0. , 0.16666667, 0.08333333, # may vary 614 0.05 ]) 615 >>> P.polyint(c,k=3) # should return array([3, 1, 1, 1]) 616 array([3., 1., 1., 1.]) 617 >>> P.polyint(c,lbnd=-2) # should return array([6, 1, 1, 1]) 618 array([6., 1., 1., 1.]) 619 >>> P.polyint(c,scl=-2) # should return array([0, -2, -2, -2]) 620 array([ 0., -2., -2., -2.]) 621 622 """ 623 c = np.array(c, ndmin=1, copy=True) 624 if c.dtype.char in '?bBhHiIlLqQpP': 625 # astype doesn't preserve mask attribute. 626 c = c + 0.0 627 cdt = c.dtype 628 if not np.iterable(k): 629 k = [k] 630 cnt = pu._deprecate_as_int(m, "the order of integration") 631 iaxis = pu._deprecate_as_int(axis, "the axis") 632 if cnt < 0: 633 raise ValueError("The order of integration must be non-negative") 634 if len(k) > cnt: 635 raise ValueError("Too many integration constants") 636 if np.ndim(lbnd) != 0: 637 raise ValueError("lbnd must be a scalar.") 638 if np.ndim(scl) != 0: 639 raise ValueError("scl must be a scalar.") 640 iaxis = normalize_axis_index(iaxis, c.ndim) 641 642 if cnt == 0: 643 return c 644 645 k = list(k) + [0]*(cnt - len(k)) 646 c = np.moveaxis(c, iaxis, 0) 647 for i in range(cnt): 648 n = len(c) 649 c *= scl 650 if n == 1 and np.all(c[0] == 0): 651 c[0] += k[i] 652 else: 653 tmp = np.empty((n + 1,) + c.shape[1:], dtype=cdt) 654 tmp[0] = c[0]*0 655 tmp[1] = c[0] 656 for j in range(1, n): 657 tmp[j + 1] = c[j]/(j + 1) 658 tmp[0] += k[i] - polyval(lbnd, tmp) 659 c = tmp 660 c = np.moveaxis(c, 0, iaxis) 661 return c 662 663 664def polyval(x, c, tensor=True): 665 """ 666 Evaluate a polynomial at points x. 667 668 If `c` is of length `n + 1`, this function returns the value 669 670 .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n 671 672 The parameter `x` is converted to an array only if it is a tuple or a 673 list, otherwise it is treated as a scalar. In either case, either `x` 674 or its elements must support multiplication and addition both with 675 themselves and with the elements of `c`. 676 677 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If 678 `c` is multidimensional, then the shape of the result depends on the 679 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + 680 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that 681 scalars have shape (,). 682 683 Trailing zeros in the coefficients will be used in the evaluation, so 684 they should be avoided if efficiency is a concern. 685 686 Parameters 687 ---------- 688 x : array_like, compatible object 689 If `x` is a list or tuple, it is converted to an ndarray, otherwise 690 it is left unchanged and treated as a scalar. In either case, `x` 691 or its elements must support addition and multiplication with 692 with themselves and with the elements of `c`. 693 c : array_like 694 Array of coefficients ordered so that the coefficients for terms of 695 degree n are contained in c[n]. If `c` is multidimensional the 696 remaining indices enumerate multiple polynomials. In the two 697 dimensional case the coefficients may be thought of as stored in 698 the columns of `c`. 699 tensor : boolean, optional 700 If True, the shape of the coefficient array is extended with ones 701 on the right, one for each dimension of `x`. Scalars have dimension 0 702 for this action. The result is that every column of coefficients in 703 `c` is evaluated for every element of `x`. If False, `x` is broadcast 704 over the columns of `c` for the evaluation. This keyword is useful 705 when `c` is multidimensional. The default value is True. 706 707 .. versionadded:: 1.7.0 708 709 Returns 710 ------- 711 values : ndarray, compatible object 712 The shape of the returned array is described above. 713 714 See Also 715 -------- 716 polyval2d, polygrid2d, polyval3d, polygrid3d 717 718 Notes 719 ----- 720 The evaluation uses Horner's method. 721 722 Examples 723 -------- 724 >>> from numpy.polynomial.polynomial import polyval 725 >>> polyval(1, [1,2,3]) 726 6.0 727 >>> a = np.arange(4).reshape(2,2) 728 >>> a 729 array([[0, 1], 730 [2, 3]]) 731 >>> polyval(a, [1,2,3]) 732 array([[ 1., 6.], 733 [17., 34.]]) 734 >>> coef = np.arange(4).reshape(2,2) # multidimensional coefficients 735 >>> coef 736 array([[0, 1], 737 [2, 3]]) 738 >>> polyval([1,2], coef, tensor=True) 739 array([[2., 4.], 740 [4., 7.]]) 741 >>> polyval([1,2], coef, tensor=False) 742 array([2., 7.]) 743 744 """ 745 c = np.array(c, ndmin=1, copy=False) 746 if c.dtype.char in '?bBhHiIlLqQpP': 747 # astype fails with NA 748 c = c + 0.0 749 if isinstance(x, (tuple, list)): 750 x = np.asarray(x) 751 if isinstance(x, np.ndarray) and tensor: 752 c = c.reshape(c.shape + (1,)*x.ndim) 753 754 c0 = c[-1] + x*0 755 for i in range(2, len(c) + 1): 756 c0 = c[-i] + c0*x 757 return c0 758 759 760def polyvalfromroots(x, r, tensor=True): 761 """ 762 Evaluate a polynomial specified by its roots at points x. 763 764 If `r` is of length `N`, this function returns the value 765 766 .. math:: p(x) = \\prod_{n=1}^{N} (x - r_n) 767 768 The parameter `x` is converted to an array only if it is a tuple or a 769 list, otherwise it is treated as a scalar. In either case, either `x` 770 or its elements must support multiplication and addition both with 771 themselves and with the elements of `r`. 772 773 If `r` is a 1-D array, then `p(x)` will have the same shape as `x`. If `r` 774 is multidimensional, then the shape of the result depends on the value of 775 `tensor`. If `tensor is ``True`` the shape will be r.shape[1:] + x.shape; 776 that is, each polynomial is evaluated at every value of `x`. If `tensor` is 777 ``False``, the shape will be r.shape[1:]; that is, each polynomial is 778 evaluated only for the corresponding broadcast value of `x`. Note that 779 scalars have shape (,). 780 781 .. versionadded:: 1.12 782 783 Parameters 784 ---------- 785 x : array_like, compatible object 786 If `x` is a list or tuple, it is converted to an ndarray, otherwise 787 it is left unchanged and treated as a scalar. In either case, `x` 788 or its elements must support addition and multiplication with 789 with themselves and with the elements of `r`. 790 r : array_like 791 Array of roots. If `r` is multidimensional the first index is the 792 root index, while the remaining indices enumerate multiple 793 polynomials. For instance, in the two dimensional case the roots 794 of each polynomial may be thought of as stored in the columns of `r`. 795 tensor : boolean, optional 796 If True, the shape of the roots array is extended with ones on the 797 right, one for each dimension of `x`. Scalars have dimension 0 for this 798 action. The result is that every column of coefficients in `r` is 799 evaluated for every element of `x`. If False, `x` is broadcast over the 800 columns of `r` for the evaluation. This keyword is useful when `r` is 801 multidimensional. The default value is True. 802 803 Returns 804 ------- 805 values : ndarray, compatible object 806 The shape of the returned array is described above. 807 808 See Also 809 -------- 810 polyroots, polyfromroots, polyval 811 812 Examples 813 -------- 814 >>> from numpy.polynomial.polynomial import polyvalfromroots 815 >>> polyvalfromroots(1, [1,2,3]) 816 0.0 817 >>> a = np.arange(4).reshape(2,2) 818 >>> a 819 array([[0, 1], 820 [2, 3]]) 821 >>> polyvalfromroots(a, [-1, 0, 1]) 822 array([[-0., 0.], 823 [ 6., 24.]]) 824 >>> r = np.arange(-2, 2).reshape(2,2) # multidimensional coefficients 825 >>> r # each column of r defines one polynomial 826 array([[-2, -1], 827 [ 0, 1]]) 828 >>> b = [-2, 1] 829 >>> polyvalfromroots(b, r, tensor=True) 830 array([[-0., 3.], 831 [ 3., 0.]]) 832 >>> polyvalfromroots(b, r, tensor=False) 833 array([-0., 0.]) 834 """ 835 r = np.array(r, ndmin=1, copy=False) 836 if r.dtype.char in '?bBhHiIlLqQpP': 837 r = r.astype(np.double) 838 if isinstance(x, (tuple, list)): 839 x = np.asarray(x) 840 if isinstance(x, np.ndarray): 841 if tensor: 842 r = r.reshape(r.shape + (1,)*x.ndim) 843 elif x.ndim >= r.ndim: 844 raise ValueError("x.ndim must be < r.ndim when tensor == False") 845 return np.prod(x - r, axis=0) 846 847 848def polyval2d(x, y, c): 849 """ 850 Evaluate a 2-D polynomial at points (x, y). 851 852 This function returns the value 853 854 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * x^i * y^j 855 856 The parameters `x` and `y` are converted to arrays only if they are 857 tuples or a lists, otherwise they are treated as a scalars and they 858 must have the same shape after conversion. In either case, either `x` 859 and `y` or their elements must support multiplication and addition both 860 with themselves and with the elements of `c`. 861 862 If `c` has fewer than two dimensions, ones are implicitly appended to 863 its shape to make it 2-D. The shape of the result will be c.shape[2:] + 864 x.shape. 865 866 Parameters 867 ---------- 868 x, y : array_like, compatible objects 869 The two dimensional series is evaluated at the points `(x, y)`, 870 where `x` and `y` must have the same shape. If `x` or `y` is a list 871 or tuple, it is first converted to an ndarray, otherwise it is left 872 unchanged and, if it isn't an ndarray, it is treated as a scalar. 873 c : array_like 874 Array of coefficients ordered so that the coefficient of the term 875 of multi-degree i,j is contained in `c[i,j]`. If `c` has 876 dimension greater than two the remaining indices enumerate multiple 877 sets of coefficients. 878 879 Returns 880 ------- 881 values : ndarray, compatible object 882 The values of the two dimensional polynomial at points formed with 883 pairs of corresponding values from `x` and `y`. 884 885 See Also 886 -------- 887 polyval, polygrid2d, polyval3d, polygrid3d 888 889 Notes 890 ----- 891 892 .. versionadded:: 1.7.0 893 894 """ 895 return pu._valnd(polyval, c, x, y) 896 897 898def polygrid2d(x, y, c): 899 """ 900 Evaluate a 2-D polynomial on the Cartesian product of x and y. 901 902 This function returns the values: 903 904 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * a^i * b^j 905 906 where the points `(a, b)` consist of all pairs formed by taking 907 `a` from `x` and `b` from `y`. The resulting points form a grid with 908 `x` in the first dimension and `y` in the second. 909 910 The parameters `x` and `y` are converted to arrays only if they are 911 tuples or a lists, otherwise they are treated as a scalars. In either 912 case, either `x` and `y` or their elements must support multiplication 913 and addition both with themselves and with the elements of `c`. 914 915 If `c` has fewer than two dimensions, ones are implicitly appended to 916 its shape to make it 2-D. The shape of the result will be c.shape[2:] + 917 x.shape + y.shape. 918 919 Parameters 920 ---------- 921 x, y : array_like, compatible objects 922 The two dimensional series is evaluated at the points in the 923 Cartesian product of `x` and `y`. If `x` or `y` is a list or 924 tuple, it is first converted to an ndarray, otherwise it is left 925 unchanged and, if it isn't an ndarray, it is treated as a scalar. 926 c : array_like 927 Array of coefficients ordered so that the coefficients for terms of 928 degree i,j are contained in ``c[i,j]``. If `c` has dimension 929 greater than two the remaining indices enumerate multiple sets of 930 coefficients. 931 932 Returns 933 ------- 934 values : ndarray, compatible object 935 The values of the two dimensional polynomial at points in the Cartesian 936 product of `x` and `y`. 937 938 See Also 939 -------- 940 polyval, polyval2d, polyval3d, polygrid3d 941 942 Notes 943 ----- 944 945 .. versionadded:: 1.7.0 946 947 """ 948 return pu._gridnd(polyval, c, x, y) 949 950 951def polyval3d(x, y, z, c): 952 """ 953 Evaluate a 3-D polynomial at points (x, y, z). 954 955 This function returns the values: 956 957 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * x^i * y^j * z^k 958 959 The parameters `x`, `y`, and `z` are converted to arrays only if 960 they are tuples or a lists, otherwise they are treated as a scalars and 961 they must have the same shape after conversion. In either case, either 962 `x`, `y`, and `z` or their elements must support multiplication and 963 addition both with themselves and with the elements of `c`. 964 965 If `c` has fewer than 3 dimensions, ones are implicitly appended to its 966 shape to make it 3-D. The shape of the result will be c.shape[3:] + 967 x.shape. 968 969 Parameters 970 ---------- 971 x, y, z : array_like, compatible object 972 The three dimensional series is evaluated at the points 973 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If 974 any of `x`, `y`, or `z` is a list or tuple, it is first converted 975 to an ndarray, otherwise it is left unchanged and if it isn't an 976 ndarray it is treated as a scalar. 977 c : array_like 978 Array of coefficients ordered so that the coefficient of the term of 979 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension 980 greater than 3 the remaining indices enumerate multiple sets of 981 coefficients. 982 983 Returns 984 ------- 985 values : ndarray, compatible object 986 The values of the multidimensional polynomial on points formed with 987 triples of corresponding values from `x`, `y`, and `z`. 988 989 See Also 990 -------- 991 polyval, polyval2d, polygrid2d, polygrid3d 992 993 Notes 994 ----- 995 996 .. versionadded:: 1.7.0 997 998 """ 999 return pu._valnd(polyval, c, x, y, z) 1000 1001 1002def polygrid3d(x, y, z, c): 1003 """ 1004 Evaluate a 3-D polynomial on the Cartesian product of x, y and z. 1005 1006 This function returns the values: 1007 1008 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * a^i * b^j * c^k 1009 1010 where the points `(a, b, c)` consist of all triples formed by taking 1011 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form 1012 a grid with `x` in the first dimension, `y` in the second, and `z` in 1013 the third. 1014 1015 The parameters `x`, `y`, and `z` are converted to arrays only if they 1016 are tuples or a lists, otherwise they are treated as a scalars. In 1017 either case, either `x`, `y`, and `z` or their elements must support 1018 multiplication and addition both with themselves and with the elements 1019 of `c`. 1020 1021 If `c` has fewer than three dimensions, ones are implicitly appended to 1022 its shape to make it 3-D. The shape of the result will be c.shape[3:] + 1023 x.shape + y.shape + z.shape. 1024 1025 Parameters 1026 ---------- 1027 x, y, z : array_like, compatible objects 1028 The three dimensional series is evaluated at the points in the 1029 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a 1030 list or tuple, it is first converted to an ndarray, otherwise it is 1031 left unchanged and, if it isn't an ndarray, it is treated as a 1032 scalar. 1033 c : array_like 1034 Array of coefficients ordered so that the coefficients for terms of 1035 degree i,j are contained in ``c[i,j]``. If `c` has dimension 1036 greater than two the remaining indices enumerate multiple sets of 1037 coefficients. 1038 1039 Returns 1040 ------- 1041 values : ndarray, compatible object 1042 The values of the two dimensional polynomial at points in the Cartesian 1043 product of `x` and `y`. 1044 1045 See Also 1046 -------- 1047 polyval, polyval2d, polygrid2d, polyval3d 1048 1049 Notes 1050 ----- 1051 1052 .. versionadded:: 1.7.0 1053 1054 """ 1055 return pu._gridnd(polyval, c, x, y, z) 1056 1057 1058def polyvander(x, deg): 1059 """Vandermonde matrix of given degree. 1060 1061 Returns the Vandermonde matrix of degree `deg` and sample points 1062 `x`. The Vandermonde matrix is defined by 1063 1064 .. math:: V[..., i] = x^i, 1065 1066 where `0 <= i <= deg`. The leading indices of `V` index the elements of 1067 `x` and the last index is the power of `x`. 1068 1069 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the 1070 matrix ``V = polyvander(x, n)``, then ``np.dot(V, c)`` and 1071 ``polyval(x, c)`` are the same up to roundoff. This equivalence is 1072 useful both for least squares fitting and for the evaluation of a large 1073 number of polynomials of the same degree and sample points. 1074 1075 Parameters 1076 ---------- 1077 x : array_like 1078 Array of points. The dtype is converted to float64 or complex128 1079 depending on whether any of the elements are complex. If `x` is 1080 scalar it is converted to a 1-D array. 1081 deg : int 1082 Degree of the resulting matrix. 1083 1084 Returns 1085 ------- 1086 vander : ndarray. 1087 The Vandermonde matrix. The shape of the returned matrix is 1088 ``x.shape + (deg + 1,)``, where the last index is the power of `x`. 1089 The dtype will be the same as the converted `x`. 1090 1091 See Also 1092 -------- 1093 polyvander2d, polyvander3d 1094 1095 """ 1096 ideg = pu._deprecate_as_int(deg, "deg") 1097 if ideg < 0: 1098 raise ValueError("deg must be non-negative") 1099 1100 x = np.array(x, copy=False, ndmin=1) + 0.0 1101 dims = (ideg + 1,) + x.shape 1102 dtyp = x.dtype 1103 v = np.empty(dims, dtype=dtyp) 1104 v[0] = x*0 + 1 1105 if ideg > 0: 1106 v[1] = x 1107 for i in range(2, ideg + 1): 1108 v[i] = v[i-1]*x 1109 return np.moveaxis(v, 0, -1) 1110 1111 1112def polyvander2d(x, y, deg): 1113 """Pseudo-Vandermonde matrix of given degrees. 1114 1115 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample 1116 points `(x, y)`. The pseudo-Vandermonde matrix is defined by 1117 1118 .. math:: V[..., (deg[1] + 1)*i + j] = x^i * y^j, 1119 1120 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of 1121 `V` index the points `(x, y)` and the last index encodes the powers of 1122 `x` and `y`. 1123 1124 If ``V = polyvander2d(x, y, [xdeg, ydeg])``, then the columns of `V` 1125 correspond to the elements of a 2-D coefficient array `c` of shape 1126 (xdeg + 1, ydeg + 1) in the order 1127 1128 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ... 1129 1130 and ``np.dot(V, c.flat)`` and ``polyval2d(x, y, c)`` will be the same 1131 up to roundoff. This equivalence is useful both for least squares 1132 fitting and for the evaluation of a large number of 2-D polynomials 1133 of the same degrees and sample points. 1134 1135 Parameters 1136 ---------- 1137 x, y : array_like 1138 Arrays of point coordinates, all of the same shape. The dtypes 1139 will be converted to either float64 or complex128 depending on 1140 whether any of the elements are complex. Scalars are converted to 1141 1-D arrays. 1142 deg : list of ints 1143 List of maximum degrees of the form [x_deg, y_deg]. 1144 1145 Returns 1146 ------- 1147 vander2d : ndarray 1148 The shape of the returned matrix is ``x.shape + (order,)``, where 1149 :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same 1150 as the converted `x` and `y`. 1151 1152 See Also 1153 -------- 1154 polyvander, polyvander3d, polyval2d, polyval3d 1155 1156 """ 1157 return pu._vander_nd_flat((polyvander, polyvander), (x, y), deg) 1158 1159 1160def polyvander3d(x, y, z, deg): 1161 """Pseudo-Vandermonde matrix of given degrees. 1162 1163 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample 1164 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`, 1165 then The pseudo-Vandermonde matrix is defined by 1166 1167 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = x^i * y^j * z^k, 1168 1169 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading 1170 indices of `V` index the points `(x, y, z)` and the last index encodes 1171 the powers of `x`, `y`, and `z`. 1172 1173 If ``V = polyvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns 1174 of `V` correspond to the elements of a 3-D coefficient array `c` of 1175 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order 1176 1177 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},... 1178 1179 and ``np.dot(V, c.flat)`` and ``polyval3d(x, y, z, c)`` will be the 1180 same up to roundoff. This equivalence is useful both for least squares 1181 fitting and for the evaluation of a large number of 3-D polynomials 1182 of the same degrees and sample points. 1183 1184 Parameters 1185 ---------- 1186 x, y, z : array_like 1187 Arrays of point coordinates, all of the same shape. The dtypes will 1188 be converted to either float64 or complex128 depending on whether 1189 any of the elements are complex. Scalars are converted to 1-D 1190 arrays. 1191 deg : list of ints 1192 List of maximum degrees of the form [x_deg, y_deg, z_deg]. 1193 1194 Returns 1195 ------- 1196 vander3d : ndarray 1197 The shape of the returned matrix is ``x.shape + (order,)``, where 1198 :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will 1199 be the same as the converted `x`, `y`, and `z`. 1200 1201 See Also 1202 -------- 1203 polyvander, polyvander3d, polyval2d, polyval3d 1204 1205 Notes 1206 ----- 1207 1208 .. versionadded:: 1.7.0 1209 1210 """ 1211 return pu._vander_nd_flat((polyvander, polyvander, polyvander), (x, y, z), deg) 1212 1213 1214def polyfit(x, y, deg, rcond=None, full=False, w=None): 1215 """ 1216 Least-squares fit of a polynomial to data. 1217 1218 Return the coefficients of a polynomial of degree `deg` that is the 1219 least squares fit to the data values `y` given at points `x`. If `y` is 1220 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple 1221 fits are done, one for each column of `y`, and the resulting 1222 coefficients are stored in the corresponding columns of a 2-D return. 1223 The fitted polynomial(s) are in the form 1224 1225 .. math:: p(x) = c_0 + c_1 * x + ... + c_n * x^n, 1226 1227 where `n` is `deg`. 1228 1229 Parameters 1230 ---------- 1231 x : array_like, shape (`M`,) 1232 x-coordinates of the `M` sample (data) points ``(x[i], y[i])``. 1233 y : array_like, shape (`M`,) or (`M`, `K`) 1234 y-coordinates of the sample points. Several sets of sample points 1235 sharing the same x-coordinates can be (independently) fit with one 1236 call to `polyfit` by passing in for `y` a 2-D array that contains 1237 one data set per column. 1238 deg : int or 1-D array_like 1239 Degree(s) of the fitting polynomials. If `deg` is a single integer 1240 all terms up to and including the `deg`'th term are included in the 1241 fit. For NumPy versions >= 1.11.0 a list of integers specifying the 1242 degrees of the terms to include may be used instead. 1243 rcond : float, optional 1244 Relative condition number of the fit. Singular values smaller 1245 than `rcond`, relative to the largest singular value, will be 1246 ignored. The default value is ``len(x)*eps``, where `eps` is the 1247 relative precision of the platform's float type, about 2e-16 in 1248 most cases. 1249 full : bool, optional 1250 Switch determining the nature of the return value. When ``False`` 1251 (the default) just the coefficients are returned; when ``True``, 1252 diagnostic information from the singular value decomposition (used 1253 to solve the fit's matrix equation) is also returned. 1254 w : array_like, shape (`M`,), optional 1255 Weights. If not None, the contribution of each point 1256 ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the 1257 weights are chosen so that the errors of the products ``w[i]*y[i]`` 1258 all have the same variance. The default value is None. 1259 1260 .. versionadded:: 1.5.0 1261 1262 Returns 1263 ------- 1264 coef : ndarray, shape (`deg` + 1,) or (`deg` + 1, `K`) 1265 Polynomial coefficients ordered from low to high. If `y` was 2-D, 1266 the coefficients in column `k` of `coef` represent the polynomial 1267 fit to the data in `y`'s `k`-th column. 1268 1269 [residuals, rank, singular_values, rcond] : list 1270 These values are only returned if `full` = True 1271 1272 resid -- sum of squared residuals of the least squares fit 1273 rank -- the numerical rank of the scaled Vandermonde matrix 1274 sv -- singular values of the scaled Vandermonde matrix 1275 rcond -- value of `rcond`. 1276 1277 For more details, see `numpy.linalg.lstsq`. 1278 1279 Raises 1280 ------ 1281 RankWarning 1282 Raised if the matrix in the least-squares fit is rank deficient. 1283 The warning is only raised if `full` == False. The warnings can 1284 be turned off by: 1285 1286 >>> import warnings 1287 >>> warnings.simplefilter('ignore', np.RankWarning) 1288 1289 See Also 1290 -------- 1291 numpy.polynomial.chebyshev.chebfit 1292 numpy.polynomial.legendre.legfit 1293 numpy.polynomial.laguerre.lagfit 1294 numpy.polynomial.hermite.hermfit 1295 numpy.polynomial.hermite_e.hermefit 1296 polyval : Evaluates a polynomial. 1297 polyvander : Vandermonde matrix for powers. 1298 numpy.linalg.lstsq : Computes a least-squares fit from the matrix. 1299 scipy.interpolate.UnivariateSpline : Computes spline fits. 1300 1301 Notes 1302 ----- 1303 The solution is the coefficients of the polynomial `p` that minimizes 1304 the sum of the weighted squared errors 1305 1306 .. math :: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2, 1307 1308 where the :math:`w_j` are the weights. This problem is solved by 1309 setting up the (typically) over-determined matrix equation: 1310 1311 .. math :: V(x) * c = w * y, 1312 1313 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the 1314 coefficients to be solved for, `w` are the weights, and `y` are the 1315 observed values. This equation is then solved using the singular value 1316 decomposition of `V`. 1317 1318 If some of the singular values of `V` are so small that they are 1319 neglected (and `full` == ``False``), a `RankWarning` will be raised. 1320 This means that the coefficient values may be poorly determined. 1321 Fitting to a lower order polynomial will usually get rid of the warning 1322 (but may not be what you want, of course; if you have independent 1323 reason(s) for choosing the degree which isn't working, you may have to: 1324 a) reconsider those reasons, and/or b) reconsider the quality of your 1325 data). The `rcond` parameter can also be set to a value smaller than 1326 its default, but the resulting fit may be spurious and have large 1327 contributions from roundoff error. 1328 1329 Polynomial fits using double precision tend to "fail" at about 1330 (polynomial) degree 20. Fits using Chebyshev or Legendre series are 1331 generally better conditioned, but much can still depend on the 1332 distribution of the sample points and the smoothness of the data. If 1333 the quality of the fit is inadequate, splines may be a good 1334 alternative. 1335 1336 Examples 1337 -------- 1338 >>> np.random.seed(123) 1339 >>> from numpy.polynomial import polynomial as P 1340 >>> x = np.linspace(-1,1,51) # x "data": [-1, -0.96, ..., 0.96, 1] 1341 >>> y = x**3 - x + np.random.randn(len(x)) # x^3 - x + N(0,1) "noise" 1342 >>> c, stats = P.polyfit(x,y,3,full=True) 1343 >>> np.random.seed(123) 1344 >>> c # c[0], c[2] should be approx. 0, c[1] approx. -1, c[3] approx. 1 1345 array([ 0.01909725, -1.30598256, -0.00577963, 1.02644286]) # may vary 1346 >>> stats # note the large SSR, explaining the rather poor results 1347 [array([ 38.06116253]), 4, array([ 1.38446749, 1.32119158, 0.50443316, # may vary 1348 0.28853036]), 1.1324274851176597e-014] 1349 1350 Same thing without the added noise 1351 1352 >>> y = x**3 - x 1353 >>> c, stats = P.polyfit(x,y,3,full=True) 1354 >>> c # c[0], c[2] should be "very close to 0", c[1] ~= -1, c[3] ~= 1 1355 array([-6.36925336e-18, -1.00000000e+00, -4.08053781e-16, 1.00000000e+00]) 1356 >>> stats # note the minuscule SSR 1357 [array([ 7.46346754e-31]), 4, array([ 1.38446749, 1.32119158, # may vary 1358 0.50443316, 0.28853036]), 1.1324274851176597e-014] 1359 1360 """ 1361 return pu._fit(polyvander, x, y, deg, rcond, full, w) 1362 1363 1364def polycompanion(c): 1365 """ 1366 Return the companion matrix of c. 1367 1368 The companion matrix for power series cannot be made symmetric by 1369 scaling the basis, so this function differs from those for the 1370 orthogonal polynomials. 1371 1372 Parameters 1373 ---------- 1374 c : array_like 1375 1-D array of polynomial coefficients ordered from low to high 1376 degree. 1377 1378 Returns 1379 ------- 1380 mat : ndarray 1381 Companion matrix of dimensions (deg, deg). 1382 1383 Notes 1384 ----- 1385 1386 .. versionadded:: 1.7.0 1387 1388 """ 1389 # c is a trimmed copy 1390 [c] = pu.as_series([c]) 1391 if len(c) < 2: 1392 raise ValueError('Series must have maximum degree of at least 1.') 1393 if len(c) == 2: 1394 return np.array([[-c[0]/c[1]]]) 1395 1396 n = len(c) - 1 1397 mat = np.zeros((n, n), dtype=c.dtype) 1398 bot = mat.reshape(-1)[n::n+1] 1399 bot[...] = 1 1400 mat[:, -1] -= c[:-1]/c[-1] 1401 return mat 1402 1403 1404def polyroots(c): 1405 """ 1406 Compute the roots of a polynomial. 1407 1408 Return the roots (a.k.a. "zeros") of the polynomial 1409 1410 .. math:: p(x) = \\sum_i c[i] * x^i. 1411 1412 Parameters 1413 ---------- 1414 c : 1-D array_like 1415 1-D array of polynomial coefficients. 1416 1417 Returns 1418 ------- 1419 out : ndarray 1420 Array of the roots of the polynomial. If all the roots are real, 1421 then `out` is also real, otherwise it is complex. 1422 1423 See Also 1424 -------- 1425 numpy.polynomial.chebyshev.chebroots 1426 numpy.polynomial.legendre.legroots 1427 numpy.polynomial.laguerre.lagroots 1428 numpy.polynomial.hermite.hermroots 1429 numpy.polynomial.hermite_e.hermeroots 1430 1431 Notes 1432 ----- 1433 The root estimates are obtained as the eigenvalues of the companion 1434 matrix, Roots far from the origin of the complex plane may have large 1435 errors due to the numerical instability of the power series for such 1436 values. Roots with multiplicity greater than 1 will also show larger 1437 errors as the value of the series near such points is relatively 1438 insensitive to errors in the roots. Isolated roots near the origin can 1439 be improved by a few iterations of Newton's method. 1440 1441 Examples 1442 -------- 1443 >>> import numpy.polynomial.polynomial as poly 1444 >>> poly.polyroots(poly.polyfromroots((-1,0,1))) 1445 array([-1., 0., 1.]) 1446 >>> poly.polyroots(poly.polyfromroots((-1,0,1))).dtype 1447 dtype('float64') 1448 >>> j = complex(0,1) 1449 >>> poly.polyroots(poly.polyfromroots((-j,0,j))) 1450 array([ 0.00000000e+00+0.j, 0.00000000e+00+1.j, 2.77555756e-17-1.j]) # may vary 1451 1452 """ 1453 # c is a trimmed copy 1454 [c] = pu.as_series([c]) 1455 if len(c) < 2: 1456 return np.array([], dtype=c.dtype) 1457 if len(c) == 2: 1458 return np.array([-c[0]/c[1]]) 1459 1460 # rotated companion matrix reduces error 1461 m = polycompanion(c)[::-1,::-1] 1462 r = la.eigvals(m) 1463 r.sort() 1464 return r 1465 1466 1467# 1468# polynomial class 1469# 1470 1471class Polynomial(ABCPolyBase): 1472 """A power series class. 1473 1474 The Polynomial class provides the standard Python numerical methods 1475 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the 1476 attributes and methods listed in the `ABCPolyBase` documentation. 1477 1478 Parameters 1479 ---------- 1480 coef : array_like 1481 Polynomial coefficients in order of increasing degree, i.e., 1482 ``(1, 2, 3)`` give ``1 + 2*x + 3*x**2``. 1483 domain : (2,) array_like, optional 1484 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped 1485 to the interval ``[window[0], window[1]]`` by shifting and scaling. 1486 The default value is [-1, 1]. 1487 window : (2,) array_like, optional 1488 Window, see `domain` for its use. The default value is [-1, 1]. 1489 1490 .. versionadded:: 1.6.0 1491 1492 """ 1493 # Virtual Functions 1494 _add = staticmethod(polyadd) 1495 _sub = staticmethod(polysub) 1496 _mul = staticmethod(polymul) 1497 _div = staticmethod(polydiv) 1498 _pow = staticmethod(polypow) 1499 _val = staticmethod(polyval) 1500 _int = staticmethod(polyint) 1501 _der = staticmethod(polyder) 1502 _fit = staticmethod(polyfit) 1503 _line = staticmethod(polyline) 1504 _roots = staticmethod(polyroots) 1505 _fromroots = staticmethod(polyfromroots) 1506 1507 # Virtual properties 1508 domain = np.array(polydomain) 1509 window = np.array(polydomain) 1510 basis_name = None 1511 1512 @classmethod 1513 def _str_term_unicode(cls, i, arg_str): 1514 return f"·{arg_str}{i.translate(cls._superscript_mapping)}" 1515 1516 @staticmethod 1517 def _str_term_ascii(i, arg_str): 1518 return f" {arg_str}**{i}" 1519 1520 @staticmethod 1521 def _repr_latex_term(i, arg_str, needs_parens): 1522 if needs_parens: 1523 arg_str = rf"\left({arg_str}\right)" 1524 if i == 0: 1525 return '1' 1526 elif i == 1: 1527 return arg_str 1528 else: 1529 return f"{arg_str}^{{{i}}}" 1530