Multivariate normal form solver

Ph. Trebuchet

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

\[ \left\{ \begin{array}{l} f_1 (x_1, \ldots, x_n) = 0\\ \vdots\\ f_s (x_1, \ldots, x_n) = 0 \end{array} \right. \]

in the case where such a system has a finite number of solutions in $\mathbbm{C}^n$. The method developped in [PhT02] is implemented here. It proceeds in two steps

The stopping criterion for these generalized normal forms $N$ is the following:

Theorem: Let $B$ be a vector space in $R = \mathbbm{R} [x_1, \ldots, x_n]$ connected to the constant polynomial $1$ {{Any monomial $\mathbf{x}^{\alpha} \neq 1 \in B$ is of the form $x_i \mathbf{x}^{\beta}$ with $\mathbf{x}^{\beta} \in B$ and some $i$ in $\{1, \ldots, n\}$.}}. If $B^+$ is the vector subspace generated by $B \cup x_1 \hspace{0.25em} B \cup \ldots \cup x_n \hspace{0.25em} B$, $N : B^+ \rightarrow B$ is a linear map such that $N$ is the identity on $B$, we define for $i = 1, \ldots, n,$ the maps

\begin{eqnarray*} M_i : B & \rightarrow & B\\ b & \mapsto & M_i (b) : = N (x_i b) . \end{eqnarray*}

The two following properties are equivalent:

  • For all $1 \leq i, j \leq n$, $M_i \circ M_j = M_j \circ M_i$.
  • $R = B \oplus I$, where $I$ is the ideal generated by the kernel of $N$
If this holds, the $B$-reduction along $\ker (N)$ is canonical.

From this normal form $N$, we deduce the roots of the system as follows. We use the properties of the operators of multiplication by elements of $\mathcal{A} = R / (f_1, .., f_s))$. For any $a \in \mathcal{A}$, we define

\begin{eqnarray*} M_a : \mathcal{A} & \rightarrow & \mathcal{A}\\ b & \mapsto & M_a (b) \equiv a b \end{eqnarray*}

We also consider its transpose operator

\begin{eqnarray*} M_a^t : \widehat{\mathcal{A}} & \rightarrow & \widehat{\mathcal{A}}\\ \Lambda & \mapsto & M_a (\Lambda) \equiv \Lambda \circ M_a \end{eqnarray*}

The matrix of $M_a^{\mathtt{t}}$ in the dual basis of a basis $B$ of $\mathcal{A}$ is the transpose of the matrix of $M_a$ in $B$.

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

Theorem: [AuSt88] [BMr98]Assume that $\mathcal{Z} (I) =\{\zeta_1, \ldots, \zeta_d \}$. We have

  • Let $a \in \mathcal{A}$. The eigenvalues of the operator $M_a$ (and its transpose $M^{\mathtt{t}}_a$) are $a (\zeta_1), \ldots, a (\zeta_d)$.
  • The common eigenvectors of $(M^{\mathtt{t}}_a)_{a \in \mathcal{A}}$ are (up to a scalar) $\mathbf{1}_{\zeta_1}, \ldots, \mathbf{1}_{\zeta_d}$.

Since $\mathbf{x}^E = ( \mathbf{x}^{\alpha})_{\alpha \in E}$ is a basis of $\mathcal{A}$, the coordinates of $\mathbf{1}_{\zeta_i}$ in the dual basis of $\mathbf{x}^E$ are $(\zeta_i^{\alpha})_{\alpha \in E}$. Thus if $\mathbf{x}^E$ contains $1, x_1, \ldots, x_n$ (which is often the case), we deduce the following algorithm:

simple roots)}

Let $a \in \mathcal{A}$ such that $a (\zeta_i) \neq a (\zeta_j)$ if $i \neq j$ (which is generically the case) and $M_a$ be the multiplication matrix by $a$ in the basis $\mathbf{x}^E = (1, x_1, \ldots, x_n, \ldots)$ of $\mathcal{A}$.

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}}.


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

See synaps/msolve/Newmac.h.



#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).