Add Support for Polynomial-Polynomial Addition/Subtraction in the Newton Basis

In extending Minterpy capabilities to carry out polynomial algebra, we should add polynomial-polynomial addition/subtraction. Basically, we would like to define the following operation

>>> nwt_poly_1
<minterpy.polynomials.newton_polynomial.NewtonPolynomial at 0x17f6bc460>
>>> nwt_poly_2
<minterpy.polynomials.newton_polynomial.NewtonPolynomial at 0x17f6bc5b0>
>>> nwt_poly_1 + nwt_poly_2
<minterpy.polynomials.newton_polynomial.NewtonPolynomial at 0x15ef3be80>

Currently, polynomial in the canonical basis and Lagrange may be added or subtracted with a polynomial in the canonical basis and respectively, Lagrange basis (though the later is still carried out inconsistently; a separate issue will be opened). However, polynomials of a given basis are closed under addition/subtraction. This means that a polynomial in any supported basis may be added with another polynomial in the same basis and the result is still a polynomial (of the same basis). This rule should apply to all basis.

However, the Lagrange or Newton basis (interpolatory bases) are different than say the canonical or Chebyshev basis in that they depend on the generating points of the interpolating grid. Even if the generating function for these points are the same, the generating points may not be nested (in general they won't be). The default Chebyshev-Lobatto nodes, for instance are not nested. This means that the basis of Newton polynomial of a given degree differs depending the maximum degree in the multi-index set. For instance, consider a one-dimensional polynomial with the following multi-index set:

MultiIndexSet(
  array([[0],
         [1],
         [2],
         [3],
         [4]]),
  lp_degree=1.0
)

and another with the following multi-index set:

MultiIndexSet(
  array([[0],
         [1],
         [2],
         [3]]),
  lp_degree=1.0
)

The generating points for the first are:

1.0
-1.0
-6.123233995736766e-17
0.7071067811865475
-0.7071067811865476

while the second are:

1.0
-1.0
0.49999999999999983
-0.5000000000000001

This means the Newton polynomial basis of the 3nd-degree for the first polynomial is:

(x - 1) (x + 1) (x - 0)

but the same basis for the second polynomial is:

(x - 1) (x + 1) (x - 0.5)

These two bases are not the same even if both are Chebyshev-Lobatto nodes. Therefore, polynomial/addition subtraction method adopted for the canonical or (later on) Chebyshev basis where we add/subtract the coefficients whose the multi-index set elements match cannot be used for the Newton or Lagrange basis.

So what to do? Below are the steps required to add/subtract Newton polynomial by a another Newton polynomial:

• Create a union between the two multi-index sets (see Issue #124 (closed))
• Create a grid of the union multi-index set
• Compute the value at the unisolvent nodes
• Add the values together (they become the Lagrange coefficients)
• Do a DDS on the Lagrange coefficients to obtain the Newton coefficients
• Create a new instance of Newton polynomial with the new multi-index set and coefficients

Note that the subtraction is addition by the additive inverse (i.e., the negative) so we don't need a separate implementation for that.

By using DDS, we assume that the resulting multi-index set is downward-closed. We may relax this, but it would complicate the code and is not how Newton polynomials are used anyhow (they are interpolating polynomial that must be defined on a downward-closed multi-index set). Furthermore, we also assume that the generating function for both polynomials are the same; otherwise, an exception should be raised. Currently, however, we have no way to check for this.