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>

6 files changed, 420 insertions, 69 deletions

diff --git a/doc/release/upcoming_changes/16154.new_feature.rst b/doc/release/upcoming_changes/16154.new_feature.rst
new file mode 100644
index 000000000..99d4b1b04
--- /dev/null
+++ b/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.
diff --git a/doc/source/reference/routines.polynomials.classes.rst b/doc/source/reference/routines.polynomials.classes.rst
index 2ce29d9d0..05e1c5476 100644
--- a/doc/source/reference/routines.polynomials.classes.rst
+++ b/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
diff --git a/doc/source/reference/routines.polynomials.rst b/doc/source/reference/routines.polynomials.rst
index 4aea963c0..6ad692e70 100644
--- a/doc/source/reference/routines.polynomials.rst
+++ b/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
diff --git a/numpy/polynomial/_polybase.py b/numpy/polynomial/_polybase.py
index 155d72805..6382732dc 100644
--- a/numpy/polynomial/_polybase.py
+++ b/numpy/polynomial/_polybase.py
@@ -37,6 +37,12 @@ class ABCPolyBase(abc.ABC):
     window : (2,) array_like, optional
         Window, see domain for its use. The default value is the
         derived class window.
+    symbol : str, optional
+        Symbol used to represent the independent variable in string 
+        representations of the polynomial expression, e.g. for printing.
+        The symbol must be a valid Python identifier. Default value is 'x'.
+
+        .. versionadded:: 1.24
 
     Attributes
     ----------
@@ -46,6 +52,8 @@ class ABCPolyBase(abc.ABC):
         Domain that is mapped to window.
     window : (2,) ndarray
         Window that domain is mapped to.
+    symbol : str
+        Symbol representing the independent variable.
 
     Class Attributes
     ----------------
@@ -100,6 +108,10 @@ class ABCPolyBase(abc.ABC):
     _use_unicode = not os.name == 'nt'
 
     @property
+    def symbol(self):
+        return self._symbol
+
+    @property
     @abc.abstractmethod
     def domain(self):
         pass
@@ -284,10 +296,12 @@ class ABCPolyBase(abc.ABC):
                 raise TypeError("Domains differ")
             elif not np.all(self.window == other.window):
                 raise TypeError("Windows differ")
+            elif self.symbol != other.symbol:
+                raise ValueError("Polynomial symbols differ")
             return other.coef
         return other
 
-    def __init__(self, coef, domain=None, window=None):
+    def __init__(self, coef, domain=None, window=None, symbol='x'):
         [coef] = pu.as_series([coef], trim=False)
         self.coef = coef
 
@@ -303,12 +317,27 @@ class ABCPolyBase(abc.ABC):
                 raise ValueError("Window has wrong number of elements.")
             self.window = window
 
+        # Validation for symbol
+        try:
+            if not symbol.isidentifier():
+                raise ValueError(
+                    "Symbol string must be a valid Python identifier"
+                )
+        # If a user passes in something other than a string, the above
+        # results in an AttributeError. Catch this and raise a more
+        # informative exception
+        except AttributeError:
+            raise TypeError("Symbol must be a non-empty string")
+
+        self._symbol = symbol
+
     def __repr__(self):
         coef = repr(self.coef)[6:-1]
         domain = repr(self.domain)[6:-1]
         window = repr(self.window)[6:-1]
         name = self.__class__.__name__
-        return f"{name}({coef}, domain={domain}, window={window})"
+        return (f"{name}({coef}, domain={domain}, window={window}, "
+                f"symbol='{self.symbol}')")
 
     def __format__(self, fmt_str):
         if fmt_str == '':
@@ -353,7 +382,7 @@ class ABCPolyBase(abc.ABC):
             except TypeError:
                 next_term = f"+ {coef}"
             # Polynomial term
-            next_term += term_method(power, "x")
+            next_term += term_method(power, self.symbol)
             # Length of the current line with next term added
             line_len = len(out.split('\n')[-1]) + len(next_term)
             # If not the last term in the polynomial, it will be two
@@ -412,18 +441,18 @@ class ABCPolyBase(abc.ABC):
         # get the scaled argument string to the basis functions
         off, scale = self.mapparms()
         if off == 0 and scale == 1:
-            term = 'x'
+            term = self.symbol
             needs_parens = False
         elif scale == 1:
-            term = f"{self._repr_latex_scalar(off)} + x"
+            term = f"{self._repr_latex_scalar(off)} + {self.symbol}"
             needs_parens = True
         elif off == 0:
-            term = f"{self._repr_latex_scalar(scale)}x"
+            term = f"{self._repr_latex_scalar(scale)}{self.symbol}"
             needs_parens = True
         else:
             term = (
                 f"{self._repr_latex_scalar(off)} + "
-                f"{self._repr_latex_scalar(scale)}x"
+                f"{self._repr_latex_scalar(scale)}{self.symbol}"
             )
             needs_parens = True
 
@@ -459,7 +488,7 @@ class ABCPolyBase(abc.ABC):
             # in case somehow there are no coefficients at all
             body = '0'
 
-        return rf"$x \mapsto {body}$"
+        return rf"${self.symbol} \mapsto {body}$"
 
 
 
@@ -470,6 +499,7 @@ class ABCPolyBase(abc.ABC):
         ret['coef'] = self.coef.copy()
         ret['domain'] = self.domain.copy()
         ret['window'] = self.window.copy()
+        ret['symbol'] = self.symbol.copy()
         return ret
 
     def __setstate__(self, dict):
@@ -491,7 +521,9 @@ class ABCPolyBase(abc.ABC):
     # Numeric properties.
 
     def __neg__(self):
-        return self.__class__(-self.coef, self.domain, self.window)
+        return self.__class__(
+            -self.coef, self.domain, self.window, self.symbol
+        )
 
     def __pos__(self):
         return self
@@ -502,7 +534,7 @@ class ABCPolyBase(abc.ABC):
             coef = self._add(self.coef, othercoef)
         except Exception:
             return NotImplemented
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def __sub__(self, other):
         othercoef = self._get_coefficients(other)
@@ -510,7 +542,7 @@ class ABCPolyBase(abc.ABC):
             coef = self._sub(self.coef, othercoef)
         except Exception:
             return NotImplemented
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def __mul__(self, other):
         othercoef = self._get_coefficients(other)
@@ -518,7 +550,7 @@ class ABCPolyBase(abc.ABC):
             coef = self._mul(self.coef, othercoef)
         except Exception:
             return NotImplemented
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def __truediv__(self, other):
         # there is no true divide if the rhs is not a Number, although it
@@ -551,13 +583,13 @@ class ABCPolyBase(abc.ABC):
             raise
         except Exception:
             return NotImplemented
-        quo = self.__class__(quo, self.domain, self.window)
-        rem = self.__class__(rem, self.domain, self.window)
+        quo = self.__class__(quo, self.domain, self.window, self.symbol)
+        rem = self.__class__(rem, self.domain, self.window, self.symbol)
         return quo, rem
 
     def __pow__(self, other):
         coef = self._pow(self.coef, other, maxpower=self.maxpower)
-        res = self.__class__(coef, self.domain, self.window)
+        res = self.__class__(coef, self.domain, self.window, self.symbol)
         return res
 
     def __radd__(self, other):
@@ -565,21 +597,21 @@ class ABCPolyBase(abc.ABC):
             coef = self._add(other, self.coef)
         except Exception:
             return NotImplemented
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def __rsub__(self, other):
         try:
             coef = self._sub(other, self.coef)
         except Exception:
             return NotImplemented
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def __rmul__(self, other):
         try:
             coef = self._mul(other, self.coef)
         except Exception:
             return NotImplemented
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def __rdiv__(self, other):
         # set to __floordiv__ /.
@@ -609,8 +641,8 @@ class ABCPolyBase(abc.ABC):
             raise
         except Exception:
             return NotImplemented
-        quo = self.__class__(quo, self.domain, self.window)
-        rem = self.__class__(rem, self.domain, self.window)
+        quo = self.__class__(quo, self.domain, self.window, self.symbol)
+        rem = self.__class__(rem, self.domain, self.window, self.symbol)
         return quo, rem
 
     def __eq__(self, other):
@@ -618,7 +650,8 @@ class ABCPolyBase(abc.ABC):
                np.all(self.domain == other.domain) and
                np.all(self.window == other.window) and
                (self.coef.shape == other.coef.shape) and
-               np.all(self.coef == other.coef))
+               np.all(self.coef == other.coef) and
+               (self.symbol == other.symbol))
         return res
 
     def __ne__(self, other):
@@ -637,7 +670,7 @@ class ABCPolyBase(abc.ABC):
             Copy of self.
 
         """
-        return self.__class__(self.coef, self.domain, self.window)
+        return self.__class__(self.coef, self.domain, self.window, self.symbol)
 
     def degree(self):
         """The degree of the series.
@@ -698,7 +731,7 @@ class ABCPolyBase(abc.ABC):
 
         """
         coef = pu.trimcoef(self.coef, tol)
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def truncate(self, size):
         """Truncate series to length `size`.
@@ -727,7 +760,7 @@ class ABCPolyBase(abc.ABC):
             coef = self.coef
         else:
             coef = self.coef[:isize]
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def convert(self, domain=None, kind=None, window=None):
         """Convert series to a different kind and/or domain and/or window.
@@ -764,7 +797,7 @@ class ABCPolyBase(abc.ABC):
             domain = kind.domain
         if window is None:
             window = kind.window
-        return self(kind.identity(domain, window=window))
+        return self(kind.identity(domain, window=window, symbol=self.symbol))
 
     def mapparms(self):
         """Return the mapping parameters.
@@ -826,7 +859,7 @@ class ABCPolyBase(abc.ABC):
         else:
             lbnd = off + scl*lbnd
         coef = self._int(self.coef, m, k, lbnd, 1./scl)
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def deriv(self, m=1):
         """Differentiate.
@@ -848,7 +881,7 @@ class ABCPolyBase(abc.ABC):
         """
         off, scl = self.mapparms()
         coef = self._der(self.coef, m, scl)
-        return self.__class__(coef, self.domain, self.window)
+        return self.__class__(coef, self.domain, self.window, self.symbol)
 
     def roots(self):
         """Return the roots of the series polynomial.
@@ -899,7 +932,7 @@ class ABCPolyBase(abc.ABC):
 
     @classmethod
     def fit(cls, x, y, deg, domain=None, rcond=None, full=False, w=None,
-        window=None):
+        window=None, symbol='x'):
         """Least squares fit to data.
 
         Return a series instance that is the least squares fit to the data
@@ -948,6 +981,8 @@ class ABCPolyBase(abc.ABC):
             value is the default class domain
 
             .. versionadded:: 1.6.0
+        symbol : str, optional
+            Symbol representing the independent variable. Default is 'x'.
 
         Returns
         -------
@@ -980,13 +1015,15 @@ class ABCPolyBase(abc.ABC):
         res = cls._fit(xnew, y, deg, w=w, rcond=rcond, full=full)
         if full:
             [coef, status] = res
-            return cls(coef, domain=domain, window=window), status
+            return (
+                cls(coef, domain=domain, window=window, symbol=symbol), status
+            )
         else:
             coef = res
-            return cls(coef, domain=domain, window=window)
+            return cls(coef, domain=domain, window=window, symbol=symbol)
 
     @classmethod
-    def fromroots(cls, roots, domain=[], window=None):
+    def fromroots(cls, roots, domain=[], window=None, symbol='x'):
         """Return series instance that has the specified roots.
 
         Returns a series representing the product
@@ -1004,6 +1041,8 @@ class ABCPolyBase(abc.ABC):
         window : {None, array_like}, optional
             Window for the returned series. If None the class window is
             used. The default is None.
+        symbol : str, optional
+            Symbol representing the independent variable. Default is 'x'.
 
         Returns
         -------
@@ -1024,10 +1063,10 @@ class ABCPolyBase(abc.ABC):
         off, scl = pu.mapparms(domain, window)
         rnew = off + scl*roots
         coef = cls._fromroots(rnew) / scl**deg
-        return cls(coef, domain=domain, window=window)
+        return cls(coef, domain=domain, window=window, symbol=symbol)
 
     @classmethod
-    def identity(cls, domain=None, window=None):
+    def identity(cls, domain=None, window=None, symbol='x'):
         """Identity function.
 
         If ``p`` is the returned series, then ``p(x) == x`` for all
@@ -1044,6 +1083,8 @@ class ABCPolyBase(abc.ABC):
             ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
             the window. If None is given then the class window is used. The
             default is None.
+        symbol : str, optional
+            Symbol representing the independent variable. Default is 'x'.
 
         Returns
         -------
@@ -1057,10 +1098,10 @@ class ABCPolyBase(abc.ABC):
             window = cls.window
         off, scl = pu.mapparms(window, domain)
         coef = cls._line(off, scl)
-        return cls(coef, domain, window)
+        return cls(coef, domain, window, symbol)
 
     @classmethod
-    def basis(cls, deg, domain=None, window=None):
+    def basis(cls, deg, domain=None, window=None, symbol='x'):
         """Series basis polynomial of degree `deg`.
 
         Returns the series representing the basis polynomial of degree `deg`.
@@ -1080,6 +1121,8 @@ class ABCPolyBase(abc.ABC):
             ``[beg, end]``, where ``beg`` and ``end`` are the endpoints of
             the window. If None is given then the class window is used. The
             default is None.
+        symbol : str, optional
+            Symbol representing the independent variable. Default is 'x'.
 
         Returns
         -------
@@ -1096,7 +1139,7 @@ class ABCPolyBase(abc.ABC):
 
         if ideg != deg or ideg < 0:
             raise ValueError("deg must be non-negative integer")
-        return cls([0]*ideg + [1], domain, window)
+        return cls([0]*ideg + [1], domain, window, symbol)
 
     @classmethod
     def cast(cls, series, domain=None, window=None):
diff --git a/numpy/polynomial/tests/test_printing.py b/numpy/polynomial/tests/test_printing.py
index 4e9902a69..0c4316223 100644
--- a/numpy/polynomial/tests/test_printing.py
+++ b/numpy/polynomial/tests/test_printing.py
@@ -309,35 +309,66 @@ class TestFormat:
             format(p, '.2f')
 
 
+@pytest.mark.parametrize(('poly', 'tgt'), (
+    (poly.Polynomial, '1.0 + 2.0·z¹ + 3.0·z²'),
+    (poly.Chebyshev, '1.0 + 2.0·T₁(z) + 3.0·T₂(z)'),
+    (poly.Hermite, '1.0 + 2.0·H₁(z) + 3.0·H₂(z)'),
+    (poly.HermiteE, '1.0 + 2.0·He₁(z) + 3.0·He₂(z)'),
+    (poly.Laguerre, '1.0 + 2.0·L₁(z) + 3.0·L₂(z)'),
+    (poly.Legendre, '1.0 + 2.0·P₁(z) + 3.0·P₂(z)'),
+))
+def test_symbol(poly, tgt):
+    p = poly([1, 2, 3], symbol='z')
+    assert_equal(f"{p:unicode}", tgt)
+
+
 class TestRepr:
     def test_polynomial_str(self):
         res = repr(poly.Polynomial([0, 1]))
-        tgt = 'Polynomial([0., 1.], domain=[-1,  1], window=[-1,  1])'
+        tgt = (
+            "Polynomial([0., 1.], domain=[-1,  1], window=[-1,  1], "
+            "symbol='x')"
+        )
         assert_equal(res, tgt)
 
     def test_chebyshev_str(self):
         res = repr(poly.Chebyshev([0, 1]))
-        tgt = 'Chebyshev([0., 1.], domain=[-1,  1], window=[-1,  1])'
+        tgt = (
+            "Chebyshev([0., 1.], domain=[-1,  1], window=[-1,  1], "
+            "symbol='x')"
+        )
         assert_equal(res, tgt)
 
     def test_legendre_repr(self):
         res = repr(poly.Legendre([0, 1]))
-        tgt = 'Legendre([0., 1.], domain=[-1,  1], window=[-1,  1])'
+        tgt = (
+            "Legendre([0., 1.], domain=[-1,  1], window=[-1,  1], "
+            "symbol='x')"
+        )
         assert_equal(res, tgt)
 
     def test_hermite_repr(self):
         res = repr(poly.Hermite([0, 1]))
-        tgt = 'Hermite([0., 1.], domain=[-1,  1], window=[-1,  1])'
+        tgt = (
+            "Hermite([0., 1.], domain=[-1,  1], window=[-1,  1], "
+            "symbol='x')"
+        )
         assert_equal(res, tgt)
 
     def test_hermiteE_repr(self):
         res = repr(poly.HermiteE([0, 1]))
-        tgt = 'HermiteE([0., 1.], domain=[-1,  1], window=[-1,  1])'
+        tgt = (
+            "HermiteE([0., 1.], domain=[-1,  1], window=[-1,  1], "
+            "symbol='x')"
+        )
         assert_equal(res, tgt)
 
     def test_laguerre_repr(self):
         res = repr(poly.Laguerre([0, 1]))
-        tgt = 'Laguerre([0., 1.], domain=[0, 1], window=[0, 1])'
+        tgt = (
+            "Laguerre([0., 1.], domain=[0, 1], window=[0, 1], "
+            "symbol='x')"
+        )
         assert_equal(res, tgt)
 
 
@@ -388,3 +419,39 @@ class TestLatexRepr:
         p = poly.HermiteE([1, 2, 3])
         assert_equal(self.as_latex(p),
             r'$x \mapsto 1.0\,{He}_{0}(x) + 2.0\,{He}_{1}(x) + 3.0\,{He}_{2}(x)$')
+
+    def test_symbol_basic(self):
+        # default input
+        p = poly.Polynomial([1, 2, 3], symbol='z')
+        assert_equal(self.as_latex(p),
+            r'$z \mapsto 1.0 + 2.0\,z + 3.0\,z^{2}$')
+
+        # translated input
+        p = poly.Polynomial([1, 2, 3], domain=[-2, 0], symbol='z')
+        assert_equal(
+            self.as_latex(p),
+            (
+                r'$z \mapsto 1.0 + 2.0\,\left(1.0 + z\right) + 3.0\,'
+                r'\left(1.0 + z\right)^{2}$'
+            ),
+        )
+
+        # scaled input
+        p = poly.Polynomial([1, 2, 3], domain=[-0.5, 0.5], symbol='z')
+        assert_equal(
+            self.as_latex(p),
+            (
+                r'$z \mapsto 1.0 + 2.0\,\left(2.0z\right) + 3.0\,'
+                r'\left(2.0z\right)^{2}$'
+            ),
+        )
+
+        # affine input
+        p = poly.Polynomial([1, 2, 3], domain=[-1, 0], symbol='z')
+        assert_equal(
+            self.as_latex(p),
+            (
+                r'$z \mapsto 1.0 + 2.0\,\left(1.0 + 2.0z\right) + 3.0\,'
+                r'\left(1.0 + 2.0z\right)^{2}$'
+            ),
+        )
diff --git a/numpy/polynomial/tests/test_symbol.py b/numpy/polynomial/tests/test_symbol.py
new file mode 100644
index 000000000..4ea6035ef
--- /dev/null
+++ b/numpy/polynomial/tests/test_symbol.py
@@ -0,0 +1,216 @@
+"""
+Tests related to the ``symbol`` attribute of the ABCPolyBase class.
+"""
+
+import pytest
+import numpy.polynomial as poly
+from numpy.core import array
+from numpy.testing import assert_equal, assert_raises, assert_
+
+
+class TestInit:
+    """
+    Test polynomial creation with symbol kwarg.
+    """
+    c = [1, 2, 3]
+
+    def test_default_symbol(self):
+        p = poly.Polynomial(self.c)
+        assert_equal(p.symbol, 'x')
+
+    @pytest.mark.parametrize(('bad_input', 'exception'), (
+        ('', ValueError),
+        ('3', ValueError),
+        (None, TypeError),
+        (1, TypeError),
+    ))
+    def test_symbol_bad_input(self, bad_input, exception):
+        with pytest.raises(exception):
+            p = poly.Polynomial(self.c, symbol=bad_input)
+
+    @pytest.mark.parametrize('symbol', (
+        'x',
+        'x_1',
+        'A',
+        'xyz',
+        'β',
+    ))
+    def test_valid_symbols(self, symbol):
+        """
+        Values for symbol that should pass input validation.
+        """
+        p = poly.Polynomial(self.c, symbol=symbol)
+        assert_equal(p.symbol, symbol)
+
+    def test_property(self):
+        """
+        'symbol' attribute is read only.
+        """
+        p = poly.Polynomial(self.c, symbol='x')
+        with pytest.raises(AttributeError):
+            p.symbol = 'z'
+
+    def test_change_symbol(self):
+        p = poly.Polynomial(self.c, symbol='y')
+        # Create new polynomial from p with different symbol
+        pt = poly.Polynomial(p.coef, symbol='t')
+        assert_equal(pt.symbol, 't')
+
+
+class TestUnaryOperators:
+    p = poly.Polynomial([1, 2, 3], symbol='z')
+
+    def test_neg(self):
+        n = -self.p
+        assert_equal(n.symbol, 'z')
+
+    def test_scalarmul(self):
+        out = self.p * 10
+        assert_equal(out.symbol, 'z')
+
+    def test_rscalarmul(self):
+        out = 10 * self.p
+        assert_equal(out.symbol, 'z')
+
+    def test_pow(self):
+        out = self.p ** 3
+        assert_equal(out.symbol, 'z')
+
+
+@pytest.mark.parametrize(
+    'rhs',
+    (
+        poly.Polynomial([4, 5, 6], symbol='z'),
+        array([4, 5, 6]),
+    ),
+)
+class TestBinaryOperatorsSameSymbol:
+    """
+    Ensure symbol is preserved for numeric operations on polynomials with
+    the same symbol
+    """
+    p = poly.Polynomial([1, 2, 3], symbol='z')
+
+    def test_add(self, rhs):
+        out = self.p + rhs
+        assert_equal(out.symbol, 'z')
+
+    def test_sub(self, rhs):
+        out = self.p - rhs
+        assert_equal(out.symbol, 'z')
+
+    def test_polymul(self, rhs):
+        out = self.p * rhs
+        assert_equal(out.symbol, 'z')
+
+    def test_divmod(self, rhs):
+        for out in divmod(self.p, rhs):
+            assert_equal(out.symbol, 'z')
+
+    def test_radd(self, rhs):
+        out = rhs + self.p
+        assert_equal(out.symbol, 'z')
+
+    def test_rsub(self, rhs):
+        out = rhs - self.p
+        assert_equal(out.symbol, 'z')
+
+    def test_rmul(self, rhs):
+        out = rhs * self.p
+        assert_equal(out.symbol, 'z')
+
+    def test_rdivmod(self, rhs):
+        for out in divmod(rhs, self.p):
+            assert_equal(out.symbol, 'z')
+
+
+class TestBinaryOperatorsDifferentSymbol:
+    p = poly.Polynomial([1, 2, 3], symbol='x')
+    other = poly.Polynomial([4, 5, 6], symbol='y')
+    ops = (p.__add__, p.__sub__, p.__mul__, p.__floordiv__, p.__mod__)
+
+    @pytest.mark.parametrize('f', ops)
+    def test_binops_fails(self, f):
+        assert_raises(ValueError, f, self.other)
+
+
+class TestEquality:
+    p = poly.Polynomial([1, 2, 3], symbol='x')
+
+    def test_eq(self):
+        other = poly.Polynomial([1, 2, 3], symbol='x')
+        assert_(self.p == other)
+
+    def test_neq(self):
+        other = poly.Polynomial([1, 2, 3], symbol='y')
+        assert_(not self.p == other)
+
+
+class TestExtraMethods:
+    """
+    Test other methods for manipulating/creating polynomial objects.
+    """
+    p = poly.Polynomial([1, 2, 3, 0], symbol='z')
+
+    def test_copy(self):
+        other = self.p.copy()
+        assert_equal(other.symbol, 'z')
+
+    def test_trim(self):
+        other = self.p.trim()
+        assert_equal(other.symbol, 'z')
+
+    def test_truncate(self):
+        other = self.p.truncate(2)
+        assert_equal(other.symbol, 'z')
+
+    @pytest.mark.parametrize('kwarg', (
+        {'domain': [-10, 10]},
+        {'window': [-10, 10]},
+        {'kind': poly.Chebyshev},
+    ))
+    def test_convert(self, kwarg):
+        other = self.p.convert(**kwarg)
+        assert_equal(other.symbol, 'z')
+
+    def test_integ(self):
+        other = self.p.integ()
+        assert_equal(other.symbol, 'z')
+
+    def test_deriv(self):
+        other = self.p.deriv()
+        assert_equal(other.symbol, 'z')
+
+
+def test_composition():
+    p = poly.Polynomial([3, 2, 1], symbol="t")
+    q = poly.Polynomial([5, 1, 0, -1], symbol="λ_1")
+    r = p(q)
+    assert r.symbol == "λ_1"
+
+
+#
+# Class methods that result in new polynomial class instances
+#
+
+
+def test_fit():
+    x, y = (range(10),)*2
+    p = poly.Polynomial.fit(x, y, deg=1, symbol='z')
+    assert_equal(p.symbol, 'z')
+
+
+def test_froomroots():
+    roots = [-2, 2]
+    p = poly.Polynomial.fromroots(roots, symbol='z')
+    assert_equal(p.symbol, 'z')
+
+
+def test_identity():
+    p = poly.Polynomial.identity(domain=[-1, 1], window=[5, 20], symbol='z')
+    assert_equal(p.symbol, 'z')
+
+
+def test_basis():
+    p = poly.Polynomial.basis(3, symbol='z')
+    assert_equal(p.symbol, 'z')