SYNAPS (SYmbolic Numeric ApplicationS) » Univariate Solvers

Univariate Bezier Subdivision Solver

We use the representation of a polynomial $f$

in the Bernstein basis: $f (x) = \sum_{i = 0}^d b_i B_i^d (x)$ , where $B_i^d (x) = \left(\begin{array}{c} d\\ i \end{array}\right) x^i (1 - x)^{d - i}$ .The coefficients $\mathbf{b}= [b_i]_{i = 0, \ldots, d}$ in this basis are called the control coefficients. We denote by $V (\mathbf{b})$ the number of sign variation of sequence $\mathbf{b}$ , not taking into account the $0$

of this sequence. Descartes' rule applies on this representation as follows:

Proposition: The number of sign changes $V (\mathbf{b})$ of $\mathbf{b}= [b_i]_{i = 0, \ldots, d}$ bounds the number of real roots of f on $[0, 1]$ and is equal to it modulo 2.

If $V (\mathbf{b}) = 0$ , the number of real root on [0, 1] is 0

If $V (\mathbf{b}) = 1$ , the number of real root on [0, 1] is 1.

We also use de Casteljau subdivision algorithm to subdivide the representation. This yields the following algorithm:

Algorithm: Isolation of the roots of on the interval .
INPUT: A representation $(\mathbf{b}, [a, b])$ associate with $f$ .

Compute the number of sign changes $V (\mathbf{b})$ of the coefficient sequence $\mathbf{b}$ .

If $V (\mathbf{b}) = 0$ , remove the interval.

If $V (\mathbf{b}) > 1$ and the size $|b - a|$ is greater than $\varepsilon$ , subdivide the representation into two subrepresentations, corresponding to the two halves of the input interval and apply recursively the algorithm to them.

If $V (\mathbf{b}) > 1$ and the size $|b - a| < \varepsilon$ , output the $\varepsilon$ -root $\frac{a + b}{2}$ with mutliplicity $V (\mathbf{b})$ .

OUTPUT: a list of subintervals of $[a, b]$ containing exactly one real root of $f$ or intervals of size $\varepsilon$ containing a $\varepsilon$ -multiple root.

Several implementations, based on this scheme are available. The methods which compute the roots, within a given precision are named Slv..., those which isolate the roots are named Isl.... We describe now the solvers of these families available in the library.

Bezier Standard Approximation Solver

Seq<O> solve(const UPOL& p, SlvBzStd<O,C> mth);

It computes all the real roots, within a given precision. The parameter O precises the type of the solutions, and C is the number type used to perform this computation. It is not necessarly the type of the coefficients of p, but it should be a field type. The default value of C is double so that it can be used as SlvBzStd<O>().

The precision of the computation can be given as an argument of the constructor: SlvBzStd<O,C>(1e-10). The default type for this precision $\varepsilon$ is double and the default value is 1e-6. If one want to use another type of precision (say X), one can use instead SlvBzStd<O,C,X>(xeps);

This solver perform a binary subdivision of an initial interval containing the roots. It stops the recursion when the size of the subintervals is smaller that the precision $\varepsilon$ specificied by the class SlvBzStd<O,C> and output in the result the middle of the interval. The subintervals which are removed have no sign variation. The \ intervals which are kept may have a sign variation bigger that . The sequence s contains the approximated roots in increasing order.

It can be called to isolate the roots in an interval [a,b] (where a and b are of type O): \ solve(p,SlvBzStd<O,C>(),a,b);

It can be used to compute quickly numerical approximation of roots taking C=double for instance, or with exact arithmetic (eg. rational numbers C=QQ) to get certified approximation.

See also:: synaps/usolve/bezier/SlvBzStd.h

Bezier Standard Isolation Solver

Seq<Interval<O> > solve(const UPOL& p, IslBzStd<O,C> mth);

It isolates the roots of the polynomial p, that is, it computes a sequence of intervals containg a unique root, or an interval of size $\varepsilon$ , where $\varepsilon$ is the precision defined in in the class IslBzStd<O,C>, if p has mutiple roots or a cluster of roots. This precision $\varepsilon =$ eps is specified by IslBzStd<O,C>(eps).

The result is a sequence of intervals with endpoints of type O. The type C is the number type used during the computation. It is not necessarly the type of the coefficients of p, but it should be a field type. The default value of C is double so that it can be used as IslBzStd<O>().

This solver perform the same subdivision as for the previous solver, except that it stops the recursion when the sign variation is . If there is a multiple or a cluster of roots, the recursion will be continue until the size of the subinterval is smaller that the precision $\varepsilon$ specificied by the class IslBzStd<O,C>. The subintervals which are removed have no sign variation. The intervals kept, which are os size $\varepsilon$ may have a sign variation bigger that . The sequence s contains the subintervals sorted in increasing order. These intervals may share some ending points.

It can be called to isolate the roots in an interval [a,b] (where a and b are of type O): \ solve(p,IslBzStd<O,C>(),a,b);

See also:: synaps/usolve/bezier/IslBzStd.h

Bezier Integer Approximation Solver

Seq<O> solve(const UPOL& p, SlvInteger<O> mth);

Same behavior as the solver SlvBzStd, except that the internal coeffcient type used to perform the subdivision is the integer type ZZ. The de Casteljau algorithm is adapted so that the representation of the internal polynomial remains with integer coefficients.

The same convention as before are used: O is the output type of the solutions; the precision $\varepsilon =$ eps can be specified by SlvBzInteger<O>(eps).

It can also be called to isolate the roots in an interval [a,b] (where a and b are of type O): solve(p,SlvBzInteger<O>(),a,b);

See also:: synaps/usolve/bezier/IslBzInteger.h

Bezier Bounding Approximation Solver

Seq<O> solve(const UPOL& p, SlvBzBgd<O> mth);

It computes all the real roots, within a given precision $\varepsilon$ . The parameter O specifies the output type, for the approximate solutions.

The precision $\varepsilon =$ eps of the approximation can be given as an argument of the constructor: SlvBzBdg<O>(eps). The default type for this precision $\varepsilon$ is double and the default value is 1e-6.

This solver assume that the coefficients are of the exact type. It performs a transformation in the Bernstein basis, using this exact arithmetic. Then the coefficients are scaled and rounded below and above to double. The subdivision is performed, exploiting the coefficient sign variations on this upper and lower enveloppes of the input polynomial. If these polynomial enveloppes do not allow us to approximate the roots within the precision $\varepsilon$ , we obtain a suspect subinterval (of size bigger than $\varepsilon$ ) on which the solver SlvBzInteger is applied.

See also:: synaps/usolve/bezier/SlvBzBdg.h

Rockwood Approximation Solver

O solve(const UPOL& p, SlvRkw<O,C> mth);

It computes an approximation of the first positive root of the polynomial p, using Rockwood's method. The internal type used during the computation is C, which default value is double.

The precision $\varepsilon =$ eps of the approximation is an argument of the constructor: SlvRkw<O,C>(eps). The default type for this precision $\varepsilon$ is double and the default value is 1e-6.

See also:: synaps/usolve/bezier/SlvRkw.h

Example

#include <synaps/upol.h>
#include <synaps/arithm/QQ.h>
#include <synaps/usolve/bezier/SlvBzStd.h>
#include <synaps/usolve/bezier/IslBzStd.h>
#include <synaps/usolve/bezier/SlvRkw.h>
int main()
{
  using std::cout; using std::endl;

  UPolDse<ZZ> P("2*x0^7+3*x0^6+2*x0^5+1*x0^4-1*x0^3-2*x0^2-3*x0+1");
  cout <<"Testing: "<<P<<endl;
  //| Testing: 2*x^7+3*x^6+2*x^5+1*x^4-1*x^3-2*x^2-3*x+1

  cout <<solve(P,SlvBzStd<double,double>(1e-6))<<endl;
  //| -1.35142, 0.278212, 0.880286

  cout <<solve(P,IslBzStd<double,double>(1e-6))<<endl;
  //| [ -1.35143,-1.35142 ], [ 0.278212,0.278213 ], 
  //| [ 0.880286,0.880287 ]

  cout << solve(P,SlvRkw<double>(1e-6)) <<endl;
  //| 0.278212
}

SYNAPS DOCUMENTATION