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 $f$ on the interval $[a, b]$.

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);

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

Bezier Standard Isolation Solver

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

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

Bezier Integer Approximation Solver

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

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

Bezier Bounding Approximation Solver

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

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

Rockwood Approximation Solver

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

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

Example

o+

#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
logo