Introduce a Method to whether the Domain of Two Polynomials Matches

Several algebraic operations involving polynomial, namely polynomial-polynomial multiplication and polynomial-polynomial subtraction/addition requires that the operands have a matching domain, both internal and user-defined (though not yet used across codebase). Currently, this check is carried out by the function _matching_internal_domain() implemented in minterpy.polynomials.canonical_polynomial (which itself is inappropriate).

Because such functionality shall be used for different polynomial bases (for instance in the current use case, the addition, subtraction, and multiplication of polynomial in the Lagrange basis use it), this function should be implemented as a method attached to MultivariatePolynomialSingleABC instead.

So, for example:

>>> lag_poly_1.has_matching_domain(lag_poly_2)
True

If both the internal domain and user domain of the two polynomials match.

Right now because all polynomials by default are defined on [-1, 1]^m as their internal domain, calling the method on any given polynomial with another should return True. In fact, we haven't seen a use case inside the codebase where internal_domain is defined outside [-1, 1]^m which by default is used if none is specified. In the future, however, we may relax the requirement that internal domain is [-1, 1] as not all polynomials are defined there (e.g., Hermite polynomials).

PS: Right now, we're still not sure how should user_domain be used across the code. The idea is that it contains the original (bounded) domain of an interpolated function before it is transformed to [-1, 1]^m. The transformation is assumed to be linear. But this is not used across Minterpy; we are strictly imposing that the domain of polynomial must be in [-1, 1]^m. For now, I'll just transfer the implementation of _matching_internal_domain() to the method.