BiRoots

In Matlab roots of a polynomial are computed as eigenvalues of the companion matrix. In this toolbox this approach is generalized to a systems of bivariate polynomials and the roots are computed as eigenvalues of an appropriate two-parameter eigenvalue problem.
Square matrices A,B,C form a determinantal representation (or linearization) of a bivariate polynomial p(x,y) if det(A+x*B+y*C) = p(x,y). To find solutions for the system p1(x,y)=0, p2(x,y)=0, where p1 and p2 are bivariate polynomials, we first find determinantal representations for p1 and p2 such that
det(A1 + x*B1 + y*C1) = p1(x,y),
det(A2 + x*B2 + y*C2) = p2(x,y).
This gives a two-parameter eigenvalue problem
A1*z1 + x*B1*z1 + y*C1*z1 = 0,
A2*z2 + x*B2*z2 + y*C2*z2 = 0,

where we are looking for an eigenvalue (x,y) and nonzero eigenvectors z1,z2 such that the above system is satisfied. The eigenvalues of the above two-parameter eigenvalue problem, which are computed by the MultiParEig toolbox, are the solutions of the system p1(x,y)=0, p2(x,y)=0.

Several determinantal representations for bivariate polynomials are implemented (see bipoly_detrep as the main method):
- 1: bipoly_detrep_123: Lin1 (matrices of asymptotic size n^2/4) - no computation [1]
- 2: bipoly_detrep_123: Lin2 (matrices of asymptotic size n^2/6) [1]
- 3: bipoly_detrep_123: Lin2 with special construction for degree 3 and 4 (matrices of asymptotic size n^2/6) [1]
- 4: bipoly_detrep_quantic - n x n matrices for degree 5 or less [2]
- 5: bipoly_detret_nxn - n x n matrices for a square-free polynomial [3]
- 6: bipoly_detrep_unif - (2n-1) x (2n-1) matrices with no computation [4]
Default method is 4 for degree 5 or less, 5 for degree 9 or less, and 6 otherwise

If the size of the determinantal representation is larger than the degree of the polynomial, then the obtained two-parameter eigenvalue problem is singular and is solved numerically by a staircase-type algorithm. Due to large matrices of the linearization and the difficulties of detecting the correct numerical rank in the staircase algorithm, this approach works efficiently and accurately only for polynomials of small order. It is recommended to use the toolbox only for polynomials of degree 13 or less. It can also fail for some polynomials of smaller degree, where double precision is not sufficient for this approach.

References:

[1] B. Plestenjak, M. E. Hochstenbach: Roots of bivariate polynomial systems via determinantal representations, SIAM J. Sci. Comput. 38 (2016) A765-A788
[2] A. Buckley, B. Plestenjak: Simple determinantal representations of up to quintic bivariate polynomials, arXiv.1609.00498
[3] B. Plestenjak: Minimal determinantal representations of bivariate polynomials, arXiv.1607.03969
[4] A. Boralevi, J. van Doornmalen, J. Draisma, M. E. Hochstenbach, B. Plestenjak: Uniform determinantal representations, arXiv.1607.04873

Example: To solve the system p1(x,y)=0, p2(x,y)=0, where

p1(x,y) = 1 + 2*x + 3*y + 4*x^2 + 5*x*y + 6*y^2 + 7*x^3 + 8*x^2*y + 9*x*y^2 + 10*y^3,
p2(x,y) = 10 + 9*x + 8*y + 7*x^2 + 6*x*y + 5*y^2 + 4*x^3 + 3*x^2*y + 2*x*y^2 + y^3,

use

P1 = [ 1 3 6 10; 2 5 9 0; 4 8 0 0; 7 0 0 0];
P2 = [10 8 5 1; 9 6 2 0; 7 3 0 0; 4 0 0 0];
[x,y] = biroots(P1,P2)

To obtain matrices [A,B,C] such that det(A+x*B+y*C) = p1(x,y), use

[A,B,C] = bipoly_detrep(P1).