# Multivariate normal form solver

Ph. Trebuchet

We consider here the problem of computing the solutions of a system

in the case where such a system has a finite number of solutions in . The method developped in [PhT02] is implemented here. It proceeds in two steps

• Computation of the generalised normal form modulo the ideal .
• Computation of the roots from eigencomputation.
The stopping criterion for these generalized normal forms is the following:

 Theorem: Let be a vector space in connected to the constant polynomial {{Any monomial is of the form with and some in .}}. If is the vector subspace generated by , is a linear map such that is the identity on , we define for the maps The two following properties are equivalent: For all , . , where is the ideal generated by the kernel of If this holds, the -reduction along is canonical.

From this normal form , we deduce the roots of the system as follows. We use the properties of the operators of multiplication by elements of . For any , we define

We also consider its transpose operator

The matrix of in the dual basis of a basis of is the transpose of the matrix of in .

Our approach to solve polynomial systems is based on the following fundamental theorem:

 Theorem: [AuSt88] [BMr98]Assume that . We have Let . The eigenvalues of the operator (and its transpose ) are . The common eigenvectors of are (up to a scalar) .

Since is a basis of , the coordinates of in the dual basis of are . Thus if contains (which is often the case), we deduce the following algorithm:

simple roots)}

Let such that if (which is generically the case) and be the multiplication matrix by in the basis of .

• Compute the eigenvectors of .
• For each eigenvector with , compute and return the point .
The case of multiple roots is treated similarly, by applying a strategy similar the one described in [CGT97]. For more details on the algorithm, see {{mt:gnfpss-issac-05}}.

## Implementation

The function which implements all these steps is:
template<class C>
Seq<VectDse<std::complex<double> > > solve(L,Newmac<C>);

where

• C is the type of coefficient in which we want to compute the normal form,
• L is the sequence of polynomials to solve.
See synaps/msolve/Newmac.h.

## Example

#include <synaps/arithm/gmp.h>
#include <synaps/mpol.h>
#include <synaps/base/Seq.h>
#include <synaps/msolve/Newmac.h>

typedef MPol<RR> mpol_t;

int main (int argc, char **argv)
{
using namespace std;

mpol_t p("x0^2*x1-x0*x1-1"), q("x0^2-x1^3-1");

std::list<mpol_t> I; I.push_back(p); I.push_back(q);

cout << solve(I, Newmac<RR>())<<endl;

//| [ (0.421134,0.779225) (-1.15111,0.166259) ]
//| [ (0.421134,-0.779225) (-1.15111,-0.166259) ]
//| [ (1.24525,-0.588209) (-0.0526987,1.13817) ]
//| [ (-0.424222,0.478962) (0.405487,0.957865) ]
//| [ (1.24525,0.588209) (-0.0526987,-1.13817) ]
//| [ (-0.424222,-0.478962) (0.405487,-0.957865) ]
//| [ (1.5648,3.86247e-19) (1.13148,8.13152e-19) ]
//| [ (-1.04912,-1.30962e-18) (0.465166,-1.41031e-19) ]

// The result is a sequence of vectors with  coefficients of type
// complex<double>. The columns correspond to the coordinates
// of the roots (here x0, x1). The rows correspond to the roots
// of the system (here 8 roots, with 2 which are real).
}



SYNAPS DOCUMENTATION