ENH: Add support for symbol to polynomial package (#16154)

Adds a symbol attribute to the polynomials from the np.polynomial package to allow the user to control/modify the symbol used to represent the independent variable for a polynomial expression. This attribute corresponds to the variable attribute of the poly1d class from the old np.lib.polynomial module.

Marked as draft for now as it depends on #15666 - all _str* and _repr* methods of ABCPolyBase and derived classes would need to be modified (and tested) to support this change.

Co-authored-by: Warren Weckesser <warren.weckesser@gmail.com>

doc/release/upcoming_changes/16154.new_feature.rst

@@ -0,0 +1,25 @@
+New attribute ``symbol`` added to polynomial classes
+----------------------------------------------------
+
+The polynomial classes in the ``numpy.polynomial`` package have a new
+``symbol`` attribute which is used to represent the indeterminate
+of the polynomial.
+This can be used to change the value of the variable when printing::
+
+    >>> P_y = np.polynomial.Polynomial([1, 0, -1], symbol="y")
+    >>> print(P_y)
+    1.0 + 0.0·y¹ - 1.0·y²
+
+Note that the polynomial classes only support 1D polynomials, so operations
+that involve polynomials with different symbols are disallowed when the
+result would be multivariate::
+
+    >>> P = np.polynomial.Polynomial([1, -1])  # default symbol is "x"
+    >>> P_z = np.polynomial.Polynomial([1, 1], symbol="z")
+    >>> P * P_z
+    Traceback (most recent call last)
+       ...
+    ValueError: Polynomial symbols differ
+
+The symbol can be any valid Python identifier. The default is ``symbol=x``,
+consistent with existing behavior.

doc/source/reference/routines.polynomials.classes.rst

@@ -52,7 +52,7 @@ the conventional Polynomial class because of its familiarity::
    >>> from numpy.polynomial import Polynomial as P
    >>> p = P([1,2,3])
    >>> p
-   Polynomial([1., 2., 3.], domain=[-1,  1], window=[-1,  1])
+   Polynomial([1., 2., 3.], domain=[-1,  1], window=[-1,  1], symbol='x')
 
 Note that there are three parts to the long version of the printout. The
 first is the coefficients, the second is the domain, and the third is the
@@ -92,19 +92,19 @@ we ignore them and run through the basic algebraic and arithmetic operations.
 Addition and Subtraction::
 
    >>> p + p
-   Polynomial([2., 4., 6.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([2., 4., 6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
    >>> p - p
-   Polynomial([0.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([0.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Multiplication::
 
    >>> p * p
-   Polynomial([ 1.,   4.,  10.,  12.,   9.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([ 1.,   4.,  10.,  12.,   9.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Powers::
 
    >>> p**2
-   Polynomial([ 1.,   4., 10., 12.,  9.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([ 1.,   4., 10., 12.,  9.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Division:
 
@@ -115,20 +115,20 @@ versions the '/' will only work for division by scalars. At some point it
 will be deprecated::
 
    >>> p // P([-1, 1])
-   Polynomial([5.,  3.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([5.,  3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Remainder::
 
    >>> p % P([-1, 1])
-   Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Divmod::
 
    >>> quo, rem = divmod(p, P([-1, 1]))
    >>> quo
-   Polynomial([5.,  3.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([5.,  3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
    >>> rem
-   Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Evaluation::
 
@@ -149,7 +149,7 @@ the polynomials are regarded as functions this is composition of
 functions::
 
    >>> p(p)
-   Polynomial([ 6., 16., 36., 36., 27.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([ 6., 16., 36., 36., 27.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Roots::
 
@@ -163,11 +163,11 @@ tuples, lists, arrays, and scalars are automatically cast in the arithmetic
 operations::
 
    >>> p + [1, 2, 3]
-   Polynomial([2., 4., 6.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([2., 4., 6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
    >>> [1, 2, 3] * p
-   Polynomial([ 1.,  4., 10., 12.,  9.], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([ 1.,  4., 10., 12.,  9.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
    >>> p / 2
-   Polynomial([0.5, 1. , 1.5], domain=[-1.,  1.], window=[-1.,  1.])
+   Polynomial([0.5, 1. , 1.5], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Polynomials that differ in domain, window, or class can't be mixed in
 arithmetic::
@@ -195,7 +195,7 @@ conversion of Polynomial classes among themselves is done for type, domain,
 and window casting::
 
     >>> p(T([0, 1]))
-    Chebyshev([2.5, 2. , 1.5], domain=[-1.,  1.], window=[-1.,  1.])
+    Chebyshev([2.5, 2. , 1.5], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Which gives the polynomial `p` in Chebyshev form. This works because
 :math:`T_1(x) = x` and substituting :math:`x` for :math:`x` doesn't change
@@ -215,18 +215,18 @@ Polynomial instances can be integrated and differentiated.::
     >>> from numpy.polynomial import Polynomial as P
     >>> p = P([2, 6])
     >>> p.integ()
-    Polynomial([0., 2., 3.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([0., 2., 3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
     >>> p.integ(2)
-    Polynomial([0., 0., 1., 1.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([0., 0., 1., 1.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 The first example integrates `p` once, the second example integrates it
 twice. By default, the lower bound of the integration and the integration
 constant are 0, but both can be specified.::
 
     >>> p.integ(lbnd=-1)
-    Polynomial([-1.,  2.,  3.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([-1.,  2.,  3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
     >>> p.integ(lbnd=-1, k=1)
-    Polynomial([0., 2., 3.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([0., 2., 3.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 In the first case the lower bound of the integration is set to -1 and the
 integration constant is 0. In the second the constant of integration is set
@@ -235,9 +235,9 @@ number of times the polynomial is differentiated::
 
     >>> p = P([1, 2, 3])
     >>> p.deriv(1)
-    Polynomial([2., 6.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([2., 6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
     >>> p.deriv(2)
-    Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([6.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 
 Other Polynomial Constructors
@@ -253,25 +253,25 @@ are demonstrated below::
     >>> from numpy.polynomial import Chebyshev as T
     >>> p = P.fromroots([1, 2, 3])
     >>> p
-    Polynomial([-6., 11., -6.,  1.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([-6., 11., -6.,  1.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
     >>> p.convert(kind=T)
-    Chebyshev([-9.  , 11.75, -3.  ,  0.25], domain=[-1.,  1.], window=[-1.,  1.])
+    Chebyshev([-9.  , 11.75, -3.  ,  0.25], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 The convert method can also convert domain and window::
 
     >>> p.convert(kind=T, domain=[0, 1])
-    Chebyshev([-2.4375 ,  2.96875, -0.5625 ,  0.03125], domain=[0.,  1.], window=[-1.,  1.])
+    Chebyshev([-2.4375 ,  2.96875, -0.5625 ,  0.03125], domain=[0.,  1.], window=[-1.,  1.], symbol='x')
     >>> p.convert(kind=P, domain=[0, 1])
-    Polynomial([-1.875,  2.875, -1.125,  0.125], domain=[0.,  1.], window=[-1.,  1.])
+    Polynomial([-1.875,  2.875, -1.125,  0.125], domain=[0.,  1.], window=[-1.,  1.], symbol='x')
 
 In numpy versions >= 1.7.0 the `basis` and `cast` class methods are also
 available. The cast method works like the convert method while the basis
 method returns the basis polynomial of given degree::
 
     >>> P.basis(3)
-    Polynomial([0., 0., 0., 1.], domain=[-1.,  1.], window=[-1.,  1.])
+    Polynomial([0., 0., 0., 1.], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
     >>> T.cast(p)
-    Chebyshev([-9.  , 11.75, -3. ,  0.25], domain=[-1.,  1.], window=[-1.,  1.])
+    Chebyshev([-9.  , 11.75, -3. ,  0.25], domain=[-1.,  1.], window=[-1.,  1.], symbol='x')
 
 Conversions between types can be useful, but it is *not* recommended
 for routine use. The loss of numerical precision in passing from a

doc/source/reference/routines.polynomials.rst

@@ -97,7 +97,7 @@ can't be mixed in arithmetic::
 
     >>> p1 = np.polynomial.Polynomial([1, 2, 3])
     >>> p1
-    Polynomial([1., 2., 3.], domain=[-1,  1], window=[-1,  1])
+    Polynomial([1., 2., 3.], domain=[-1,  1], window=[-1,  1], symbol='x')
     >>> p2 = np.polynomial.Polynomial([1, 2, 3], domain=[-2, 2])
     >>> p1 == p2
     False