SYNAPS (SYmbolic Numeric ApplicationS) » Tutorial

Polynomial solvers

The library provides several type of solvers of polynomial equations. The general signature for these functions is

 Seq<S> solve (E, Mth<S,...>(...))
 S      RootOf(E, Mth<S,...>(...))

where

E is the equation(s) type,

Mth<S,...>(...) is the class which specifies the solving method. The first template argument is specifying the type of solution we want (eg. a sequence of values of roots of type S. Some other parameters may be completing the solver definition. Other non-static parameters such as eg. the precision could also be given in (...).

A solve function computes a sequence of solutions, where as RootOf computes one specific root.

Univariate solvers

Univariate solvers take as input a polynomial of the type UPolDse. We illustrate their use on several exemples, now.

Suppose for instance that we are manipulating univariate polynomials with integer coefficients (ZZ) and that we want to compute the real roots by the Aberth method. Then after including the file:

#include <synaps/usolve/Aberth.h>

we can use the instruction:

 UpolDse<ZZ> p("x^7-2324325131254354*x+1");
 Seq<RR> sol= solve(p, Aberth<RR>());

We obtain

$4.30232 e^{- 16}, - 4.99422, 4.99422$

Instead, if we use Aberth<CC>(), we obtain the sequence of all the $23$ complex roots. Notice that these are approximate values of the actual roots. The way to decide wether we have a real or a complex roots depends on a threshold. The level of approximation is controled by the class Aberth.

#include <synaps/init.h>
#ifdef SYNAPS_WITH_GMP
#include <synaps/arithm/gmp.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/usolve/Aberth.h>

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  UPolDse<QQ> p("x^7-2324325131254354*x-1");
  cout<<p<<endl;
  Seq<RR> sol= solve(p, SlvAberth<RR>());
  cout<<sol<<endl;

  cout<<p(sol[1])<<endl;
}
#else
int main(int argc, char** argv) {}
#endif //SYNAPS_WITH_GMP

Exercise: Compute the real (resp. complexe) roots of the polynomial $W_n (x) = \prod_1^n (x - i) \in \mathbbm{Q}[x]$ for n=20. Same question for $W_{20} (x) + \frac{1}{10^{12}} x^{12}$ . $\boxdot$

If one is intereseted in the root of polynomial in an interval, subdivision solvers may be used. The following solver is using the Bernstein basis representation in order to approximate or isolate the real roots of a polynomial. It requires the following files:

#include <synaps/usolve/bezier/SlvBzStd.h>
#include <synaps/usolve/bezier/SlvBzBdg.h>
#include <synaps/usolve/bezier/IslBzStd.h>

and is illustrated now:

 UpolDse<ZZ> p("x^7-2324325131254354*x-1");
 cout<< solve(p, SlvBzStd<double>(1e-6)) << endl; 
 cout<< solve(p, SlvBzBdg<double>(1e-6), -400., 500.) << endl;

The parameter 1e-6 of the class SlvBzStd gives the level of approximation of the roots. In this computation we look at the roots of p in the interval $[- 400, 500]$

. The two solvers yields different approximation of the roots:

$\begin{eqnarray*} - 363.912 & - 4.76837 e^{- 07} & 363.912\\ - 363.956 & - 2.78912 e^{- 08} & 363.956 \end{eqnarray*}$

In the first case, the input coefficients are rounded to double and the subvision algorithm is applied on this approximate polynomial. Its is a fast solver over the double but should be used with caution, for certified computation when the input coefficients are big.

In the second case SlvBzBdg<double> a lower and upper envelloppe with double coefficients of the polynomial is computed. These enveloppes are used to approximate the roots or to detect if the precision is enough and to zoom if necessary. This approximation is certified.

We can simply be interested by intervals containing one root:

 cout<< solve(p, IslBzStd<QQ>(1e-6)) << endl;

This computes a sequence of intervals with end values in $\mathbbm{Q}$ such that $[a_i, b_i]$

contains only one root of p:

$- 2324325131254355, - 1, - 1, 0, 1, 2324325131254355$

The naming convention, which we use is the following: Slv... is solver which approximate the real roots, Isl... is a solver which isolate the real roots.

#include <synaps/init.h>
#ifdef SYNAPS_WITH_GMP
#include <synaps/arithm/gmp.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/usolve/bezier/SlvBzStd.h>
#include <synaps/usolve/bezier/SlvBzBdg.h>
#include <synaps/usolve/bezier/IslBzStd.h>

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  UPolDse<ZZ> p("x^7-2324325131254354*x-1");
  cout<<p<<endl;

  cout<<solve(p, SlvBzStd<double,double>())<<endl;
  cout<<solve(p, SlvBzBdg<double,QQ>(), -400., 500.)<<endl;
  cout<<solve(p, IslBzStd<QQ,QQ>())<<endl;
  
}
#else
int main(int argc, char** argv) {}
#endif //SYNAPS_WITH_GMP

Roots of polynomials over $\mathbbm{Z}$ can be manipulated exactly, as algebraic numbers. Their internal representation will be the defining polynomial and an interval of isolation. It requires the following file:

#include <synaps/usolve/Algebraic.h>

Here is an example, where we compute the real roots with rational isolating interval and evaluate the sign of a polynomial at the first root:

 UpolDse<ZZ> p("x^7-2314*x-1"),q("x^3-2*x+1");
 typedef algebraic::real<QQ> algebraic_t;
 Seq<algebraic_t> sol= solve(p, Algebraic<QQ>());
 cout<<SignAt(q,sol[0])<<endl;

It yields

$\begin{eqnarray*} \tmop{root}_\tmop{of} (x^7 - 2324325131254354) x - 1, 1) \tmop{in} & [- \frac{2324325131254355}{439804651104}, \frac{- 2324325131254355}{8796093022208}], & \tmop{approx} = - 363.956,,\\ \tmop{root}_\tmop{of} (x^7 - 2324325131254354 x - 1, 1) \tmop{in} & [0, 2324325131254355], & \tmop{approx} = 363.956,\\ \tmop{root}_\tmop{of} (x^7 - 2324325131254354) x - 1, 1) \tmop{in} & [ \frac{- 2324325131254355}{8796093022208}, 0], & \tmop{approx} = - 3.00412 e^{- 11} \end{eqnarray*}$

#include <synaps/init.h>
#ifdef SYNAPS_WITH_GMP
#include <synaps/arithm/gmp.h>
#include <synaps/upol/UPolDse.h>
#include <synaps/usolve/algebraic/sign_fct.h>
#include <synaps/usolve/Sturm.h>

typedef ALGEBRAIC::root_of<ZZ> algebraic_t;

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  UPolDse<ZZ> p("x^7-2324325131254354*x-1"), q("x^3-2*x-1");
  cout<<p<<endl;
  Seq<algebraic_t> sol = solve(p, Sturm<ZZ>());
  cout<<sol<<endl;
  cout<<sign_at(p,sol[1])<<","<<sign_at(q,sol[1])<<endl;
  
}
#else
int main(int argc, char** argv) {}
#endif //SYNAPS_WITH_GMP

Multivariate solvers

Two type of multivariate solvers are available in the library:

a normal form solver, generalising Grobner basis computation and

a subdivision solver using the representation of multivariate polynomials in the Bernstein basis.

The normal form solver can be used from the file:

#include <synaps/msolve/Newmac.h>

as follows:

 typedef MPol<RR> mpol_t;
 std::list<mpol_t> E;
 E.push_back(mpol_t("x^2+y^2+z^2-1")); 
 E.push_back(mpol_t("x^2+2*x*y-y^2+z^2-1")); 
 E.push_back(mpol_t("3*x^2-x*y+4*y^2+z^2-2*x-3*y+z-2"));

 Seq<VectDse<std:complex<long double> > > sol = solve(E, Newmac<RR>());
 cout<< sol <<endl;

It yields a sequence of vectors with complex numbers over the {{long double}}:

$\begin{eqnarray*} {}[ & (1.40281, 0.277651), (1.40281, 0.277651), (0.450884, - 1.72768) & ],\\ {}[ & (1.40281, - 0.277651), (1.40281, - 0.277651), (0.450884, 1.72768) & ],\\ {}[ & (- 0.305627, - 1.35784 e^{- 17}), (- 0.30 b 5627, - 2.26906 e^{- 17}), (- 0.901767, - 4.39707 e^{- 17}) & ],\\ {}[ & (- 0.576379, - 1.25406 e^{- 17}), (- 1.62598 e^{- 16}, 2.08405 e^{- 17}), (- 0.817183, - 2.11964 e^{- 17}) & ],\\ {}[ & (1.28819, - 0.274842), (9.78732 e^{- 17}, 3.27854 e^{- 16}), (0.408591, 0.866511) & ],\\ {}[ & (1.28819, 0.274842), (2.14442 e^{- 16}, - 6.61071 e^{- 17}), (0.408591, - 0.866511) & ],\\ {}[ & (1.40634 e^{- 09}, 1.43066 e^{- 08}), (- 9.37557 e^{- 10}, - 9.53772 e^{- 09}), (1, 2.382 e^{- 15}) & ],\\ {}[ & (- 1.40634^{- 09}, - 1.43066 e^{- 08}), (9.37557 e^{- 10}, 9.53772 e^{- 09}), (1, - 7.4936 e^{- 16}) & ] \end{eqnarray*}$

#include <synaps/init.h>
#ifdef SYNAPS_WITH_GMP
#include <synaps/arithm/gmp.h>
#include <synaps/mpol.h>
#include <synaps/msolve/Newmac.h>

typedef MPol<double> mpol_t;

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  std::list<mpol_t> E;
  E.push_back(mpol_t("x^2+y^2+z^2-1"));
  E.push_back(mpol_t("x^2+2*x*y-y^2+z^2-1"));
  E.push_back(mpol_t("3*x^2-x*y+4*y^2+z^2-2*x-3*y+z-2"));

  cout<<solve(E,Newmac<double>())<<endl;
  //  [(1.40281,0.277651),(1.40281,0.277651),(0.450884,-1.72768)], 
  //  [(1.40281,-0.277651),(1.40281,-0.277651),(0.450884,1.72768)], 
  //  [(-0.305627,-1.35784e-17),(-0.305627,-2.26906e-17),(-0.901767,-4.39707e-17)], 
  //  [(-0.576379,-1.25406e-17),(-1.62598e-16,2.08405e-17),(-0.817183,-2.11964e-17)], 
  //  [(1.28819,-0.274842),(9.78732e-17,3.27854e-16),(0.408591,0.866511)], 
  //  [(1.28819,0.274842),(2.14442e-16,-6.61071e-17),(0.408591,-0.866511)]
  //  [(1.40634e-09,1.43066e-08),(-9.37557e-10,-9.53772e-09),(1,2.382e-15)], 
  //  [(-1.40634e-09,-1.43066e-08),(9.37557e-10,9.53772e-09),(1,-7.4936e-16)], 
  
}
#else
int main(int argc, char** argv) {}
#endif //SYNAPS_WITH_GMP

The subdivision solver can be used from the file

#include <synaps/msolve/sbdslv.h>

We define the domain $[a_1, b_1] \times \cdots \times [a_n, b_n]$ where we search the roots by the vector:

  VectDse<double> domain(6,"-1 1.1 -1 1.1 -1 1.1");

and call the subdivision solver with the standard preconditionner:

 cout<<solve(E, SBDSLV<double,SBDSLV_RDRDL>(1e-6), domain)<<endl;

The first parameter double of SBDSLV is the type of coefficients during the computation,

The second parameter VectDse<double> is the type of domain,

The third parameter SBDSLV_RDRDL specifies the preconditioner,

The first argument domain specifies the domain where we search the solutions,

The second argument eps=1e-6 is precision to which we stop the subdivision.

It yields a sequence of 3 real vectors:

$\begin{eqnarray*} {}[ & 5.72114 e^{- 08} 4.53431 e^{- 08} 1 & ]\\ {}[ & - 0.576379 - 2.21964 e^{- 10} - 0.817183 & ]\\ {}[ & - 0.305627 - 0.305627 - 0.901767 & ] \end{eqnarray*}$

#include <synaps/init.h>
#ifdef SYNAPS_WITH_GMP
#include <synaps/arithm/gmp.h>
#include <synaps/mpol.h>
#include <synaps/msolve/Newmac.h>

typedef MPol<double> mpol_t;

int main(int argc, char** argv) 
{
  using std::cout;
  using std::endl;
  std::list<mpol_t> E;
  E.push_back(mpol_t("x^2+y^2+z^2-1"));
  E.push_back(mpol_t("x^2+2*x*y-y^2+z^2-1"));
  E.push_back(mpol_t("3*x^2-x*y+4*y^2+z^2-2*x-3*y+z-2"));

  cout<<solve(E,Newmac<double>())<<endl;
  //  [(1.40281,0.277651),(1.40281,0.277651),(0.450884,-1.72768)], 
  //  [(1.40281,-0.277651),(1.40281,-0.277651),(0.450884,1.72768)], 
  //  [(-0.305627,-1.35784e-17),(-0.305627,-2.26906e-17),(-0.901767,-4.39707e-17)], 
  //  [(-0.576379,-1.25406e-17),(-1.62598e-16,2.08405e-17),(-0.817183,-2.11964e-17)], 
  //  [(1.28819,-0.274842),(9.78732e-17,3.27854e-16),(0.408591,0.866511)], 
  //  [(1.28819,0.274842),(2.14442e-16,-6.61071e-17),(0.408591,-0.866511)]
  //  [(1.40634e-09,1.43066e-08),(-9.37557e-10,-9.53772e-09),(1,2.382e-15)], 
  //  [(-1.40634e-09,-1.43066e-08),(9.37557e-10,9.53772e-09),(1,-7.4936e-16)], 
  
}
#else
int main(int argc, char** argv) {}
#endif //SYNAPS_WITH_GMP

SYNAPS DOCUMENTATION