Use a better orthogonalization procedure for basis finding in quadrature finder

[inducer]

The procedure in this paper could replace this code:

modepy/modepy/quadrature/construction.py

Lines 171 to 188 in 8f06f1c

	def orthogonalize_basis(
	integrate: Callable[[Integrand], np.inexact],
	basis: Sequence[Integrand],
	) -> Sequence[LinearCombinationIntegrand]:
	r"""
	Let :math:`\Omega\subset\mathbb C^n` be a domain. (Domains
	over the reals are allowable as well.) Returns linear combinations
	of functions in *basis* that is orthogonal under the (complex-valued)
	:math:`L^2` inner product induced by *integrate*.
	:arg integrate: Computes an integral of the passed integrand over
	:math:`\Omega`. Must use integration nodes compatible with
	*basis*.
	:arg basis: A sequence of functions that accept an array of nodes
	of shape either ``(n, nnodes)`` or ``(nnodes,)`` and return
	an array of shape ``(nnodes,)`` with the value of the
	basis function at the node.
	"""

(not exactly, since the initial basis is fixed, but morally)

[inducer]

https://arxiv.org/pdf/2202.04843 offers another competitor. (pointed out by @xywei)

cc @ShawnL00