  In lattice QCD, the calculation of physical quantities from disconnected
quark loop calculations have large variance due to the use of Monte Carlo
methods for the estimation of the trace of the inverse lattice Dirac operator.
In this work, we build upon our POLY and HFPOLY variance reduction methods by
using high-degree polynomials. Previously, the GMRES polynomials used were only
stable for low-degree polynomials, but through application of a new, stable
form of the GMRES polynomial, we have achieved higher polynomial degrees than
previously used. While the variance is not dependent on the trace correction
term within the methods, the evaluation of this term will be necessary for
forming the vacuum expectation value estimates. This requires computing the
trace of high-degree polynomials, which can be evaluated stochastically through
our new Multipolynomial Monte Carlo method. With these new high-degree noise
subtraction polynomials, we obtained a variance reduction for the scalar
operator of nearly an order of magnitude over that of no subtraction on a $24^3
\times 32$ quenched lattice at $\beta = 6.0$ and $\kappa = 0.1570 \approx
\kappa_{crit}$. Additionally, we observe that for sufficiently high polynomial
degrees, POLY and HFPOLY approach the same level of effectiveness. We also
explore the viability of using double polynomials for variance reduction as a
means of reducing the required orthogonalization and memory costs associated
with forming high-degree GMRES polynomials.
