NumPy - polynomials

https://numpy.org/doc/stable/reference/generated/numpy.poly1d.html

https://numpy.org/doc/stable/reference/routines.polynomials.package.html

https://numpy.org/doc/stable/reference/routines.polynomials.classes.html

POLY1D (OLD)

# numpy.poly1d(c_or_r, r=False, variable=None)
# c_or_r : array_like, coefficients in DECREASING powers.
# r : bool, optional, the default is False.
#     If True, c_or_r specifies the polynomial’s roots.
# variable : str, optional
#     Default variable is 'x' for printing.

p = np.poly1d([1, 2, 3])
print(p)
#    2
# 1 x + 2 x + 3
p.r   # or p.roots, array([-1.+1.41421356j, -1.-1.41421356j])
p(0.5)   # 4.25, value for x = 0.5
p.c   # or p.coefficients, array([1, 2, 3])
p[1]   # 2, coefficient
p.order   # 2
p * p   # poly1d([ 1,  4, 10, 12,  9])
p ** 3 + 4   # poly1d([ 1,  6, 21, 44, 63, 54, 31])
p.deriv()   # poly1d([2, 2]), derivative
p.integ()   # poly1d([0.33333333, 1., 3., 0.]), indefinite integral

# Construct a polynomial from its roots, (x-1)*(x-2)

np.poly1d([1, 2], r=True)   # poly1d([ 1., -3.,  2.])
np.poly1d([1, -1]) * np.poly1d([1, -2])   # poly1d([ 1, -3,  2])

POLYNOMIAL (NEW)

The sub-package 'numpy.polynomial' provides 'convenience classes' 
for each of six different kinds of polynomials.
+------------+------------------+
| Name       | Provides         |
+------------+------------------+
| Polynomial | Power series     |
| Chebyshev  | Chebyshev series |
| Legendre   | Legendre series  |
| Laguerre   | Laguerre series  |
| Hermite    | Hermite series   |
| HermiteE   | HermiteE series  |
+------------+------------------+
P.Polynomial, a0 + a1*x + a2*x^2, coefficients [a0, a1, a2]
P.Chebyshev, b0 + b1*T1(x) + b2*T2(x), coefficients [b0, b1, b2]
P.Hermite, c0 + c1*H1(x) + c2*H2(x), coefficients [c0, c1, c2]

import numpy as np
from numpy import polynomial as P

one_x = P.Polynomial.identity()
# Polynomial([0., 1.], domain=[-1.,  1.], window=[-1.,  1.])
# 'domain' and 'window' are important for fitting.

P.Polynomial.basis(3)   # x^3
# Polynomial([0., 0., 0., 1.], domain=[-1.,  1.], window=[-1.,  1.])

p1 = P.Polynomial([1, 2, 3])   # 1 + 2 x + 3 x^3
# Polynomial([1., 2., 3.], domain=[-1,  1], window=[-1,  1])

p1.coef   # array([1., 2., 3.])
p1.domain   # array([-1,  1])
p1.window   # array([-1,  1])
p1.degree()   # 2, the degree of the series

p2 = P.Polynomial.fromroots([-1,1])   # -1 + x^2
# Polynomial([-1.,  0.,  1.], domain=[-1.,  1.], window=[-1.,  1.])

r2 = p2.roots()   # array([-1.,  1.])

p1(1)   # 6.0, evaluation
p1(np.array([1.5, 2.5, 3.5]))   # array([10.75, 24.75, 44.75])

p1(p2)   # substitution
# Polynomial([ 2.,  0., -4.,  0.,  3.], domain=[-1.,  1.], window=[-1.,  1.])

p1 + p2
# Polynomial([0., 2., 4.], domain=[-1.,  1.], window=[-1.,  1.])

p1 - p2
# Polynomial([2., 2., 2.], domain=[-1.,  1.], window=[-1.,  1.])

p1 * p2
# Polynomial([-1., -2., -2.,  2.,  3.], domain=[-1.,  1.], window=[-1.,  1.])

p1 // p2   # floor division
# Polynomial([3.], domain=[-1.,  1.], window=[-1.,  1.])

p1 ** 2
# Polynomial([ 1.,  4., 10., 12.,  9.], domain=[-1.,  1.], window=[-1.,  1.])

p1 % P.Polynomial.identity()   # remainder
# Polynomial([1.], domain=[-1.,  1.], window=[-1.,  1.])

p1_copy = p1.copy()   # return a copy
assert p1_copy == p1

p1.deriv()   # differentiate
# Polynomial([2., 6.], domain=[-1.,  1.], window=[-1.,  1.])

p1.cutdeg(1)   # truncate series to the given degree
# Polynomial([1., 2.], domain=[-1.,  1.], window=[-1.,  1.])

p1.integ()   # integrate, the default lower bound is 0 and the integration constant is 0
# Polynomial([0., 1., 1., 1.], domain=[-1.,  1.], window=[-1.,  1.])

p1.integ(lbnd=-1, k=3)
# Polynomial([4., 1., 1., 1.], domain=[-1.,  1.], window=[-1.,  1.])

leg2 = P.Legendre([0, 0, 1])   # -0.5 + 1.5 x^2
#leg2 = P.Legendre.basis(2)
# Legendre([0., 0., 1.], domain=[-1,  1], window=[-1,  1])

leg2(0)   # -0.5

p3 = P.Polynomial.cast(leg2)   # convert series to series of this class
# Polynomial([-0.5,  0. ,  1.5], domain=[-1.,  1.], window=[-1.,  1.])
# Warning: the loss of numerical precision is possible.

p1.convert(kind=P.Chebyshev)
# Chebyshev([2.5, 2. , 1.5], domain=[-1.,  1.], window=[-1.,  1.])

# Polynomials that differ in domain, window, or class can’t be mixed in arithmetic.
# But different types can be used for substitution.
p1(leg2)
# Legendre([1.6, 0., 2.85714286, 0., 1.54285714], domain=[-1.,  1.], window=[-1.,  1.])