Add ability to evaluate polynomials at square matrices

galois/_fields/_functions.py

@@ -231,6 +231,17 @@ def _poly_evaluate(cls, coeffs, x):
 
  return results
 
+ def _poly_evaluate_matrix(cls, coeffs, X):
+ field = cls
+ assert X.ndim == 2 and X.shape[0] == X.shape[1]
+ I = field.Identity(X.shape[0])
+
+ results = coeffs[0]*I
+ for j in range(1, coeffs.size):
+ results = coeffs[j]*I + results @ X
+
+ return results
+
  def _poly_divmod(cls, a, b):
  assert isinstance(a, cls) and isinstance(b, cls)
  assert 1 <= a.ndim <= 2 and b.ndim == 1

galois/_fields/_main.py

@@ -1826,8 +1826,9 @@ def characteristic_poly(self) -> "Poly":
  An :math:`n \times n` matrix :math:`\mathbf{A}` has characteristic polynomial
  :math:`p_A(x) = \textrm{det}(x\mathbf{I} - \mathbf{A})`.
 
- Special properties of the characteristic polynomial are the constant coefficient is
- :math:`\textrm{det}(-\mathbf{A})` and the :math:`x^{n-1}` coefficient is :math:`-\textrm{Tr}(\mathbf{A})`.
+ Special properties of the characteristic polynomial are: the constant coefficient is
+ :math:`\textrm{det}(-\mathbf{A})`, the :math:`x^{n-1}` coefficient is :math:`-\textrm{Tr}(\mathbf{A})`,
+ and :math:`p_A(\mathbf{A}) = \mathbf{0}`.
 
  References
  ----------
@@ -1847,6 +1848,8 @@ def characteristic_poly(self) -> "Poly":
  poly.coeffs[-1] == np.linalg.det(-A)
  # The x^n-1 coefficient is -Tr(A)
  poly.coeffs[1] == -np.trace(A)
+ # The characteristic polynomial annihilates the matrix A
+ poly(A, elementwise=False)
  """
  if not self.ndim == 2:
  raise ValueError(f"The array must be 2-D to compute its characteristic poly, not have shape {self.shape}.")
@@ -3380,9 +3383,9 @@ def __hash__(self):
  t = tuple([self.field.order,] + self.nonzero_degrees.tolist() + self.nonzero_coeffs.tolist())
  return hash(t)
 
- def __call__(self, x: FieldArray, field: Optional[FieldClass] = None) -> FieldArray:
+ def __call__(self, x: FieldArray, field: Optional[FieldClass] = None, elementwise: bool = True) -> FieldArray:
  """
- Evaluate the polynomial.
+ Evaluates the polynomial at :math:`x`.
 
  Parameters
  ----------
@@ -3391,21 +3394,52 @@ def __call__(self, x: FieldArray, field: Optional[FieldClass] = None) -> FieldAr
  field : galois.FieldClass, optional
  The Galois field to evaluate the polynomial over. The default is `None` which represents
  the polynomial's current field, i.e. :obj:`field`.
+ elementwise : bool, optional
+ Indicates to evaluate arrays elementwise. The default is `True`. If `False`, the polynomial
+ indeterminate is evaluated at the square matrix :math:`X`.
 
  Returns
  -------
  galois.FieldArray
- The result of the polynomial evaluation of the same shape as `x`.
+ The result of the polynomial evaluation of the same shape as :math:`x`.
+
+ Examples
+ --------
+ .. ipython:: python
+
+ GF = galois.GF(2**8)
+ p = galois.Poly([37, 123, 0, 201], field=GF); p
+
+ Evaluate the polynomial elementwise at :math:`x`.
+
+ .. ipython:: python
+
+ x = GF.Random(4); x
+ p(x)
+ GF(37)*x**3 + GF(123)*x**2 + GF(201)
+
+ Evaluate the polynomial at the matrix :math:`X`.
+
+ .. ipython:: python
+
+ X = GF.Random((2,2)); X
+ p(X, elementwise=False)
+ GF(37)*np.linalg.matrix_power(X,3) + GF(123)*np.linalg.matrix_power(X,2) + GF(201)*GF.Identity(2)
  """
- if field is None:
- field = self.field
- coeffs = self.coeffs
+ if not isinstance(field, (type(None), FieldClass)):
+ raise TypeError(f"Argument `field` must be a Galois field array class, not {type(field)}.")
+
+ field = self.field if field is None else field
+ coeffs = field(self.coeffs)
+ x = field(x)
+
+ if elementwise:
+ return field._poly_evaluate(coeffs, x)
  else:
- assert isinstance(field, FieldClass)
- coeffs = field(self.coeffs)
- if not isinstance(x, field):
- x = field(x)
- return field._poly_evaluate(coeffs, x)
+ if not (x.ndim == 2 and x.shape[0] == x.shape[1]):
+ raise ValueError(f"Argument `x` must be a square matrix when evaluating the polynomial not elementwise, not have shape {x.shape}.")
+ return field._poly_evaluate_matrix(coeffs, x)
+
  def __len__(self) -> int:
  """
  Returns the length of the coefficient array.

tests/polys/test_evaluate_matrix.py

@@ -0,0 +1,63 @@
+"""
+A pytest module to test polynomial evaluation of matrices.
+
+Sage:
+ F = GF(3**5, name="x", repr="int")
+ R = PolynomialRing(F, "x")
+ X = []
+ n = 2
+ for i in range(n):
+ row = []
+ for j in range(n):
+ row.append(F.random_element())
+ X.append(row)
+ print(X)
+ X = matrix(X)
+ for _ in range(3):
+ p = R.random_element(5)
+ print(p)
+ print(list(p(X)))
+"""
+import pytest
+import numpy as np
+
+import galois
+
+
+def test_exceptions():
+ GF = galois.GF(2**8)
+ p = galois.Poly([3,0,2,1], field=GF)
+ x = GF([101, 102, 103, 104])
+
+ with pytest.raises(TypeError):
+ p(x, field="invalid type")
+ with pytest.raises(ValueError):
+ p(x, elementwise=False)
+
+
+def test_binary_extension():
+ GF = galois.GF(2**8)
+ X = GF([[66, 232], [44, 46]])
+
+ p = galois.Poly.String("19*x^5 + 10*x^4 + 145*x^3 + 143*x^2 + 133*x + 87", field=GF)
+ assert np.array_equal(p(X, elementwise=False), GF([(178, 163), (190, 189)]))
+
+ p = galois.Poly.String("216*x^5 + 48*x^4 + 250*x^3 + 181*x^2 + 182*x + 216", field=GF)
+ assert np.array_equal(p(X, elementwise=False), GF([(243, 59), (113, 254)]))
+
+ p = galois.Poly.String("121*x^5 + 165*x^4 + 184*x^3 + 198*x^2 + 248*x + 156", field=GF)
+ assert np.array_equal(p(X, elementwise=False), GF([(24, 145), (66, 188)]))
+
+
+def test_prime_extension():
+ GF = galois.GF(3**5)
+ X = GF([[170, 221], [175, 156]])
+
+ p = galois.Poly.String("141*x^5 + 126*x^4 + 43*x^3 + 46*x^2 + 95*x + 106", field=GF)
+ assert np.array_equal(p(X, elementwise=False), GF([(50, 179), (49, 133)]))
+
+ p = galois.Poly.String("91*x^5 + 206*x^4 + 143*x^3 + 34*x^2 + 211*x + 162", field=GF)
+ assert np.array_equal(p(X, elementwise=False), GF([(89, 156), (230, 139)]))
+
+ p = galois.Poly.String("171*x^5 + 214*x^4 + 82*x^3 + x^2 + 97*x + 109", field=GF)
+ assert np.array_equal(p(X, elementwise=False), GF([(170, 230), (202, 104)]))