Remove all usage of monomial basis in bezier._algebraic_intersection

[dhermes]

By re-phrasing all computations as eigenvalue problems / rank computations, the code will be

• more stable (converting from Bernstein basis is unstable)
• more "Fortran"-izable (the numpy polynomial procedures are a not-so-thin wrapper around e.g. dgeev)

---

This is actually related to #38, as the code in locate_point() does not generalize nicely off of [0, 1] due to the use of normalize_polynomial().

[dhermes]

See: https://gist.github.com/dhermes/8c177036e426ed6fb936943ebb01b5fb

[dhermes]

As a "simple"-ish example of what can be done, consider:

f(s) = 60 b06 + 45 b16 + 32 b26 + 21 b36 + 12 b46 + 5 b56
     = 30 s^2 - 90 s + 60
     = 30 (s - 2) (s - 1)

i.e.

>>> import numpy as np
>>> import bezier
>>>
>>> nodes = np.asfortranarray([
...     [60.0],
...     [45.0],
...     [32.0],
...     [21.0],
...     [12.0],
...     [ 5.0],
...     [ 0.0],
... ])
>>> polynomial = bezier.Curve(nodes, degree=1, _copy=False)
>>> s_vals = np.arange(7.0)
>>> s_vals
array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.])
>>> def f(s):
...     return 30 * (s - 2) * (s - 1)
... 
>>> f(s_vals)
array([  60.,   -0.,    0.,   60.,  180.,  360.,  600.])
>>> polynomial.evaluate_multi(s_vals).flatten()
array([  60.,    0.,    0.,   60.,  180.,  360.,  600.])

By defining t = s / (1 - s) <==> s = t / (1 + t), we can rewrite:

f(s) = (1 - s) (1 - s)^5 [60 + 270 t + 480 t^2 + 420 t^3 + 180 t^4 + 30 t^5]

(we "donate" the binomial coefficients to the Bernstein coefficients to get "generalized Bernstein coefficients"). NOTE that one of the (1 - s) factors is distinguished from the (1 - s)^5. In the "general" case we'd just peel off (1 - s)^6, but since the coefficient of b66 is 0, we'd have a degenerate polynomial in t (i.e. degree 6 with lead term 0). Hence we already know s = 1 is a root.

Now we have a companion matrix for the polynomial in t:

M = [[ -6, -14, -16,  -9,  -2],
     [  1,   0,   0,   0,   0],
     [  0,   1,   0,   0,   0],
     [  0,   0,   1,   0,   0],
     [  0,   0,   0,   1,   0]]

• M has eigenvalues -2, -1, -1, -1, -1.
• -1 corresponds to a point at infinity of multiplicity four: s = (-1) / (1 + (-1)) = infinity. This is a symptom of the degree elevation (i.e. the true degree of the polynomial is 2, but we have a degree 6 Bézier curve / polynomial)
• t = -2 ==> s = (-2) / (1 - 2) = 2

All together, we find s = -2, -1 are the roots (without ever having formed the monomial coefficients or attempting to degree reduce).

---

$ git log -1 --pretty=%H
26b029c5cccc3dbada8c9f99f83243d7d1621950

[dhermes]

Continuing the example from before, a few notes:

---

Actually computing the eigenvalues of M from above is easier said than done:

>>> import numpy as np
>>> import scipy.linalg
>>> M = scipy.linalg.companion([30., 180., 420., 480., 270., 60.])
>>> M
array([[ -6., -14., -16.,  -9.,  -2.],
       [  1.,   0.,   0.,   0.,   0.],
       [  0.,   1.,   0.,   0.,   0.],
       [  0.,   0.,   1.,   0.,   0.],
       [  0.,   0.,   0.,   1.,   0.]])
>>>
>>> eigvals = np.linalg.eigvals(M)
>>> eigvals.real
array([-2.        , -1.00019175, -1.00019175, -0.99980825, -0.99980825])
>>> eigvals.imag
array([ 0.        ,  0.00019181, -0.00019181,  0.00019169, -0.00019169])
>>> np.abs(eigvals[1:] + 1)
array([ 0.00027121,  0.00027121,  0.00027113,  0.00027113])
>>> eigvals[0]
(-2.0000000000000311+0j)
>>> eigvals[0].real.hex()
'-0x1.0000000000046p+1'

---

In the non-elevated case we have:

f(s) = 60 b02 + 30 b12

which becomes

f(s) = (1 - s) (1 - s) [60 + 60 t]

and has a simple companion matrix: [[-1]].