A MATLAB Tool for the Stable Generation of Matrix Polynomial Evaluation Schemes with Two-Product Savings

Computing numerical approximations of matrix functions frequently relies on the efficient evaluation of high-degree matrix polynomials. Although computational bounds are historically governed by the Paterson--Stockmeyer (PS) method, recent theoretical developments have demonstrated the viability of evaluation schemes that eliminate two matrix products ($2M$). Existing literature documents stable instances of this $2M$ reduction only for isolated cases, such as specific degrees of Taylor approximations for the matrix exponential and the matrix logarithm. However, a generalized approach for arbitrary polynomials remains unestablished. To address this limitation, this work presents a software-driven procedure that extends these computational savings to polynomials of degrees $m \in \{18, 21, 24, 26, 27, 28\}$ and all $m \ge 30$, requiring primarily a non-zero leading coefficient. Since the underlying evaluation coefficients must be determined by solving systems of nonlinear equations (SNEs), selecting a numerically stable solution set is critical. We introduce an automated verification routine designed to filter and validate robust coefficient sets for floating-point execution. The primary contribution is a MATLAB implementation leveraging variable precision arithmetic to handle the underlying SNEs, verify stability, and project precision bounds. Numerical experiments involving various matrix functions verify that the developed implementation preserves or, in some instances, enhances the numerical accuracy of the classic PS method, while systematically achieving the theoretical reduction of $2M$.