# A planar quadratic clipping method for computing a root of a polynomial in an interval

## A planar quadratic clipping method for computing a root of a polynomial in an interval

Computers & Graphics 46 (2015) 89–98 Contents lists available at ScienceDirect Computers & Graphics journal homepage: www.elsevier.com/locate/cag S...
Computers & Graphics 46 (2015) 89–98

Contents lists available at ScienceDirect

Computers & Graphics journal homepage: www.elsevier.com/locate/cag

SMI 2014

A planar quadratic clipping method for computing a root of a polynomial in an interval Xiao-Diao Chen a,b,n, Weiyin Ma b a b

Hangzhou Dianzi University, Hangzhou 310037, PR China Department of MBE, City University of Hong Kong, Hong Kong, PR China

art ic l e i nf o

a b s t r a c t

Article history: Received 7 July 2014 Received in revised form 20 August 2014 Accepted 16 September 2014 Available online 2 October 2014

Keywords: Optimal approximation order Root ﬁnding Quadratic clipping Convergence rate

1. Introduction Computing the solutions of (systems of) non-linear equations is frequently needed for solving various geometric and graphics problems [1,2,8,9,12,13,15,20,22],. Typical examples include curve/surface intersection and surface rendering by ray-tracing [7,15,18,19], point projection, which is to compute the closest point on a curve or surface to a given point , collision detection [4,14], and bisectors/ medial axes computation . The method is particularly useful for the ﬂy computation in connection with the above applications in both geometry processing and graphics rendering. For solving a polynomial equation, the Descartes rule [5,21] and the Sturm theorem [10,11] are two well-known techniques for isolating real roots. A collection of related references can be found in . A Bézier curve is bounded by the convex hull of its control polygon, the corresponding roots are then bounded by the roots of the convex hull. Note that the number of the zeros of a Bézier function is less or equal to that of its control polygon, the zeros of its control polygon can be used for isolating the zeros of a Bézier function . There are also several clipping algorithms developed for ﬁnding the roots of polynomials [1,15,17,23]. In 1990, Sederberg and Nishita

n

90

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

approximate the curve ðt; f ðtÞÞ in R2 space. Three constructive interpolation methods QT1, QT2 and QT3 are proposed, in which the quadratic curve interpolates either two points and their two directional tangent vectors, or three points and one directional tangent vector. QT1 and QT2 are used for convex cases that the curve ðt; f ðtÞÞ has no inﬂexion point, while QT3 is applied for other cases. In principle, QT1 and QT2 can be used for bounding f(t) within a linear computation complexity, which is much better than that of previous clipping methods Oðn2 Þ. Numerical examples illustrate the effectiveness and the effect of the new method. The remainder of this paper is organized as follows. Section 2 introduces three constructive quadratic tangent methods with optimal approximation order that are used as the foundation of the proposed root ﬁnding method. Section 3 presents the new root ﬁnding method named planar quadratic clipping, while Section 4 discusses the computation complexity and convergence of the proposed method. Some further examples with comparisons are given in Section 5, and the conclusions are drawn at the end of this paper.

Given a smooth function f(t) and an interval ½a; b. For the sake of convenience, we transfer the parameter interval ½a; b into ½0; 1 0 0 by using f ðtÞ ¼ f ða þ htÞ and f ðtÞ ¼ f ða þhtÞh, where h ¼ b  a. Without loss of generality, suppose that the given interval is ½0; 1. Let pm ¼ ð0:5; f ð0:5ÞÞ ¼ ð0:5; pm Þ;

p1 ¼ ð1; f ð1ÞÞ ¼ ð1; p1 Þ;

v0 ¼ ð1; f ð0ÞÞ ¼ ð1; d0 Þ;

0

λ0 ¼ detðp1 p0 ; v0 Þ;

0

v1 ¼ ð1; f ð1ÞÞ ¼ ð1; d1 Þ;

λ1 ¼ detðp1  p0 ; v1 Þ;

where “det” denotes determinant of the two vectors. Suppose that the quadratic Bézier curves are Aj ðuÞ ¼ ðx j ðuÞ; y j ðuÞÞ ¼ p0 B20 ðuÞ þ qj B21 ðuÞ þ p1 B22 ðuÞ;

j ¼ 1; 2; 3;

v1 v0

p0

sponding computation complexity is linear, which is much better than those of previous clipping methods. 2.2. The ﬁrst planar quadratic clipping method (QT1) Suppose that λ0  λ1 o 0, e.g., when CðtÞ is convex within ½0; 1, see also Fig. 1. Let q1 be the intersection point between the two lines p0 þ v0 t and p1 þ v1 t, which can be computed by solving the following constraints:

0

vm ¼ ð1; f ð0:5ÞÞ ¼ ð1; dm Þ;

p1

Fig. 1. Illustration of QT1. The curve in solid black is the curve CðtÞ, while the curve in dashed blue is the quadratic Bézier curve with control points p0 , q1 and p1 . (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this article.)

2. Constructive planar quadratic clipping methods

p0 ¼ ð0; f ð0ÞÞ ¼ ð0; p0 Þ;

q1

ð1Þ

q1 ¼ p0 þ α0 v0 ¼ p1  α1 v1 ;

ð2Þ

where α0 and α1 are two unknown real numbers. From Eq. (2), we have that A1 ðiÞ ¼ pi ; A01 ðiÞ ¼ 2αi vi ;

i ¼ 0; 1;

ð3Þ

which means that A1 ðuÞ interpolates two end-points and two directional tangent vectors of the two end-points. It can be veriﬁed that the values of α0 and α1 are positive. In fact, if the length of interval ½a; b tends to be zero, then x 1 ðuÞ tends to be u, which means that the value of αi tends to be 0.5, i¼ 0, 1. In this case, let ϕ1 ðtÞ ¼ t þ ρ1;0 ðt 1Þt, y 1 ðtÞ ¼ p0 B20 ðuÞ þ ρ1;1 B21 ðuÞ þ p1 B21 ðuÞ and Y 1 ðtÞ ¼ y 1 ðϕ1 ðtÞÞ, where ρ1;i is a unknown. The equation system pi ¼ Y 1 ðiÞ;

0

di ¼ Y 1 ðiÞ;

i ¼ 0; 1;

ð4Þ

can be simpliﬁed into 8 2 > < 2ðp1  p0 Þρ1;0 þ ðd0  d1 Þρ1;0 þ 2p0 þ d1  2p1 þ d0 ¼ 0; d1 > : ρ1;1 ¼ p1  2ð1 þ ρ Þ : 1;0 There may be two values of ρ1;0 , and we select the one which is of the smaller absolute value. Let H 1 ðtÞ ¼ f ðtÞ Y 1 ðtÞ. Combing Eq. (4) with [6, Theorem 3.5.1, p. 67, Chapter 3.5], there exists ξ1 ðtÞ A ½0; 1 such that ð4Þ

f ðtÞ  Y 1 ðtÞ ¼

H 1 ðξ1 ðtÞÞ ðt 1Þ2 t 2 : 24

ð5Þ

2.3. The second planar quadratic clipping method (QT2) Let υ0 ¼ detðpm  p0 ; vm Þ, υ1 ¼ detðp1  pm ; vm Þ and υ2 ¼ detðp1  p0 ; vm Þ. Suppose that υ0  υ1 o 0, e.g., when CðtÞ is convex within

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

91

r0

vm

q2

p1

v0

pm

p1

pm q3 p0

v1 r1

p0 Fig. 2. Illustration of QT2. The curve in solid black is the curve CðtÞ, while the curve in dashed blue is the quadratic Bézier curve with control points p0 , q2 and p1 . (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this article.)

½0; 1, see also Fig. 2. We look for a quadratic curve A2 ðuÞ together with two unknowns β and um such that A2 ð0Þ ¼ p0 ;

A2 ð1Þ ¼ p1 ;

A2 ðum Þ ¼ pm ;

A02 ðum Þ ¼ β vm ;

ð6Þ

which means that A2 ðuÞ interpolates three points and one directional tangent vector at the middle point. The numerator of the divided differences of A2 ðuÞ on ½0; um ; um ; 1 becomes A2 ½0; um ; um ; 1 ¼ ðu2m  um ÞA02 ðum Þ þ ðA2 ð1Þ  A2 ð0ÞÞu2m  2ðA2 ðmÞ  A2 ð0ÞÞum þðA2 ðum Þ  A2 ð0ÞÞ ¼

ðu2m  um Þ

β

2.4. The third planar quadratic clipping method (QT3) The QT3 method is to ﬁnd A3 ðuÞ interpolating three points and one directional tangent vector at one end point. Suppose that λ0  λ1 4 0, e.g., when CðtÞ is concave within ½0; 1, see Fig. 3. Without loss of generality, suppose that λ0 4 0. Let λ2 ¼ detðpm  p0 ; v0 Þ and λ3 ¼ detðp1  pm ; v1 Þ. If pm locates inside the triangle △p0 p1 r0 such that 0 o λ2 o λ0 , we try to construct the quadratic Bézier curve A3 ðuÞ such that A3 ð0Þ ¼ p0 ;

A3 ð1Þ ¼ p1 ;

ð7Þ

ð9Þ

pm  p0 B20 ðum Þ  p1 B22 ðum Þ B21 ðum Þ

:

Similar to Section 2.2, we can ﬁnd ϕ2 ðtÞ ¼ t þ ρ2;0 ðt  1Þt, y 2 ðtÞ ¼ f ð0ÞB20 ðuÞ þ ρ2;1 B21 ðuÞ þ f ð1ÞB21 ðuÞ and Y 2 ðtÞ ¼ y 2 ðϕ2 ðtÞÞ such that p0 ¼ Y 2 ð0Þ;

pm ¼ Y 2 ð0:5Þ;

p1 ¼ Y 2 ð1Þ;

dm ¼

0 Y 2 ð0:5Þ;

¼ ðu 2m  u m Þγ 0 v0  u 2m ðp1 p0 Þ þ ðpm  p0 Þ ¼ 0:

ð10Þ

where the two unknowns ρ2;0 and ρ2;1 are determined by 8 ðp0  p1  dm Þρ22;0 þ 4ðp0 þ p1  2pm Þρ2;0 > > > > < þ4ðdm þp0  p1 Þ ¼ 0; 2ðp0  p1 þ dm Þ þ ðp0 þ p1 Þρ2;0 > > > : > : ρ2;1 ¼ 2ρ2;0

R2 ðu m Þ ¼  λ0 u 2m þ λ2 ¼ 0:

pm ¼ Y 3 ð0:5Þ

and

0

d0 ¼ Y 3 ð0Þ:

ð15Þ

0 3 ðtÞ Z0

is not satisﬁed for all t A ½0; 1, we set y 3 ðtÞ as the If ϕ y-coordinate of A3 ðuÞ, and ϕ3 ðtÞ ¼ t þ ρ3;0 ðt  1Þt=ðρ3;1 t þ 1Þ such that

ϕ3 ð1Þ ¼ 1;

ϕ3 ð0:5Þ ¼ u m and ϕ03 ð0Þ ¼ γ 0 1 ;

where u m and γ0 are determined by Eq. (12). Let Y 3 ðtÞ ¼ y 3 ðϕ3 ðtÞÞ and H 3 ðtÞ ¼ f ðtÞ  Y 3 ðtÞ. We have that

χ 3 ð0Þ ¼ 0;

ð11Þ

χ 3 ð0:5Þ ¼ 0;

χ 3 ð1Þ ¼ 0 and χ 03 ð0Þ ¼ 0;

ð16Þ

where χ 3 ¼ H 3 ; H^ 3 . Similarly, there exists ξ3 ðtÞ A ½0; 1 such that

ð4Þ

H 2 ðξ2 ðtÞÞ ðt  0:5Þ2 tðt  1Þ: 24

ð14Þ

From pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃthe assumption that 0o λ2 o λ0 , the value of u m is λ2 =λ0 A ð0; 1Þ. Similarly, let ϕ3 ðtÞ ¼ t þ ρ3;0 ðt 1Þt, y 3 ðtÞ ¼ p0 B20 ðuÞ þ ρ3;2 B21 ðuÞ þ p1 B21 ðuÞ and Y 3 ðtÞ ¼ y 3 ðϕ3 ðtÞÞ. It can be veriﬁed that ϕ3 ð0Þ ¼ 0, ϕ3 ð1Þ ¼ 1, p0 ¼ Y 3 ð0Þ and p1 ¼ Y 3 ð1Þ. In many cases, we can set ρ3;1 ¼ 0 and ﬁnd two unknowns ρ3;0 and ρ3;2 such that

ϕ3 ð0Þ ¼ 0;

Combing Eq. (10) with [6, Theorem 3.5.1, p. 67, Chapter 3.5], there exists ξ2 ðtÞ A ½0; 1 such that H 2 ðtÞ ¼ f ðtÞ  Y 2 ðtÞ ¼

ð13Þ

By taking the determinant with v0 , Eq. (13) can be simpliﬁed into

Finally, from the constraint A2 ðum Þ ¼ pm , we obtain q2 ¼

ð12Þ

þ ðA3 ðmÞ  A3 ð0ÞÞ

ð8Þ

From the assumption that υ0  υ1 o 0, we obtain 1 um ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ A ð0; 1Þ:  υ1 =υ0 þ 1

A03 ð0Þ ¼ γ 0 v0 :

A3 ½0; 0; u m ; 1 ¼ ðu 2m u m ÞA03 ð0Þ  u 2m ðA3 ð1Þ  A3 ð0ÞÞ

By taking the determinant with vm , Eq. (7) can be simpliﬁed as R2 ðmÞ ¼ υ0  2υ0 um þ υ2 u2m ¼ υ0 ð1  um Þ2 þ υ1 u2m ¼ 0:

A3 ðu m Þ ¼ pm ;

where γ0 is a positive real number determined by Eq. (12). In fact, the numerator of the divided differences of A3 ðuÞ on ½0; 0; u m ; 1 becomes

vm þ ðp1  p0 Þu2m  2ðpm  p0 Þum

þ ðpm  p0 Þ ¼ 0:

Fig. 3. Illustration of QT3. r0 (or r1 ) is the foot point of p1 (or p0 ) to the tangent line L0 ðp0 ; v0 Þ (or L1 ðp1 ; v1 Þ), the curve in solid black and the curve in dashed blue are the curve CðtÞ and the quadratic Bézier curve with control points p0 , q3 and p1 , respectively. (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this article.)

χ 3 ðtÞ ¼

χ ð4Þ 3 ðξ3 ðtÞÞ 24

t 2 ðt  0:5Þðt  1Þ;

χ 3 ¼ H 3 ; H^ 3 :

ð17Þ

92

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

Otherwise, if pm locates inside the triangle △p0 p1 r1 such that 0 o λ3 o λ1 , we try to construct A3 ð0Þ ¼ p0 ;

A3 ð1Þ ¼ p1 ;

A3 ðu m Þ ¼ pm ;

A03 ð1Þ ¼ γ 1 v1 ;

ð18Þ

where γ1 is a positive real number determined by Eq. (18). In fact, the numerator of the divided differences of A3 ðuÞ on ½0; u^ m ; 1; 1 becomes 2

ð4Þ

¼ 0:

ð19Þ

R2 ðu m Þ ¼ λ1 ðu m 1Þ2  λ3 ¼ 0:

ð20Þ

From the assumption that 0 o λ3 o λ1 , the value of u^ m is pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1  λ3 =λ1 A ð0; 1Þ. In this case, let ϕ4 ðtÞ ¼ t þ ρ4;0 ðt  1Þt, y 4 ðtÞ ¼ p0 B20 ðuÞ þ ρ4;2 B21 ðuÞ þ p1 B21 ðuÞ and Y 4 ðtÞ ¼ y 4 ðϕ4 ðtÞÞ. It can be veriﬁed that ϕ4 ð0Þ ¼ 0, ϕ4 ð1Þ ¼ 1, p0 ¼ Y 4 ð0Þ and p1 ¼ Y 4 ð1Þ. In many cases, we can ﬁnd two unknowns ρ4;0 and ρ4;2 such that 0 d1 ¼ Y 4 ð1Þ:

ð21Þ

0 4 ðtÞ Z 0

is not satisﬁed for all t A ½0; 1, we set y 4 ðtÞ as the yIf ϕ coordinate of A3 ðuÞ, and ϕ4 ðtÞ ¼ t þ ρ4;0 ðt 1Þt=ðρ4;1 t þ1Þ such that

ϕ4 ð0:5Þ ¼ u m and ϕ04 ð1Þ ¼ γ 1 1 ;

where A3 ðuÞ, u m and γ1 are determined by Eq. (18). Let Y 4 ðtÞ ¼ y 4 ðϕ4 ðtÞÞ and H 4 ðtÞ ¼ f ðtÞ  Y 4 ðtÞ. We have that

χ 4 ð0:5Þ ¼ 0;

χ 4 ð1Þ ¼ 0 and χ 04 ð1Þ ¼ 0;

ð22Þ

where χ 4 ¼ H 4 ; H^ 4 . Similarly, there exists ξ4 ðtÞ A ½0; 1 such that

χ 4 ðtÞ ¼

χ 4ð4Þ ðξ4 ðtÞÞ 24

t 2 ðt  0:5Þðt  1Þ;

χ 4 ¼ H 4 ; H^ 4 :

κ1 is a constant.

κ1 can be computed in

ð4Þ

By taking the determinant with v1 , Eq. (19) can be simpliﬁed into

χ 4 ð0Þ ¼ 0;

ðtÞ  κ 1 , where

If y 1 ðtÞ is rewritten as At þ 2Bt þ C in ð1; 3Þ,

2

ϕ4 ð1Þ ¼ 1;

ð4Þ

2

¼ ð  u^ m þ u^ m Þγ 1 v1 þ ðu^ m  1Þ2 ðp1  p0 Þ  ðp1  pm Þ

ϕ4 ð0Þ ¼ 0;

ð4Þ

Thus, Y 1 ðtÞ  κ 1 and H 1 ðtÞ ¼ f

ðA3 ð1Þ A3 ð0ÞÞ  ðA3 ð1Þ A3 ðmÞÞ

and

Without loss of generality, we assume that λ0  λ1 a 0. Suppose 0 that we obtain the values of f(i) and f ðiÞ, i ¼ a; b, which costs ð4n  2; 4n  2Þ.  Case 1: λ0  λ1 o 0. We use QT1 to obtain the reference polynomial Y 1 ðtÞ ¼ y 1 ðϕ1 ðtÞÞ,

where Y 1 ðtÞ, y 1 ðtÞ and ϕ1 ðtÞ are of degree 4, 2 and 2, respectively.

A3 ½0; u^ m ; 1; 1 ¼ ð  u^ m þ u^ m ÞA03 ð1Þ þ ðu^ m  1Þ2

pm ¼ Y 4 ð0:5Þ

3.1. The ﬁrst two steps for obtaining two bounding polynomials

ð23Þ

Fig. 3 illustrates the condition of 0 o λ2 o λ0 , and the illustration of 0 o λ3 o λ1 can be done in a similar way. In Fig. 3, the directional vector ðp1  r0 Þ is perpendicular to v0 . If the inner point pm locates into the triangle Δðp0 ; p1 ; r0 Þ, it can be veriﬁed that 0 o λ2 o λ0 . In principle, when the length h of the interval ½a; b is small enough, the middle point pm (or even the whole curve segment) locates in the quadrilateral p0 r0 p1 r1 , and thus, the quadratic curve A3 ðuÞ determined by either Eq. (12) or Eq. (18) exists. Remark 1. If p0 , p1 and pm are collinear, or λ0  λ1 ¼ 0, we simply set the quadratic curve as the line p0 ð1 tÞ þ p1 t. Usually, we can extend the idea of the planar clipping method by using cubic polynomial curves and achieve a much better approximation effect and a higher convergence rate, which is beyond the topic of this paper.

ð2; 0Þ. The differentiation f ðtÞ within ½a; b can be rewritten into Bernstein form within ð2ðn  4Þðn 3Þ; ðn  4Þðn  3ÞÞ. Thus, an upper bound

ð4Þ

τ1 and a lower bound τ2 of H 1 ðtÞ can be obtained from

ð2ðn  4Þðn  3Þ þ 3; ðn  4Þðn  3Þ þ 2Þ. Note that 0 r ðt  1Þ2 t 2 r 1 1 16 ; 8 t A ½0; 1, let K 1 ¼ 384 , τ r;1 ¼ maxf0; K 1 τ 1 g and τ l;1 ¼ min f0; K 1 τ2 g, we have the two bounding polynomials Y l;1 ðtÞ ¼ τ l;1 þY 1 ðtÞ r f ðtÞ r τ r;1 þ Y 1 ðtÞ ¼ Y r;1 ðtÞ:

ð24Þ

Suppose that u j;1 A ½0; 1 is a root of y 1 ðuÞ þ τ j;1 , then the corresponding root t ⋆ j;1 of Y j;1 ðtÞ can be computed by solving a quadratic equation

ϕ1 ðtÞ  u j;1 ¼ 0, j ¼ l; r.

Remark 2. For the case that ð4Þ

ð  1Þi H 1 ðtÞ r 0

and

ð4Þ

ð  1Þi H 2 ðtÞ r 0;

i ¼ 0 or i ¼ 1; 8 t A ½0; 1;

combining Eq. (5) with Eq. (11), we have that ð  1Þi Y 1 ðtÞ Z ð  1Þi f ðtÞ Z ð  1Þi Y 2 ðtÞ, which means that Y 1 ðtÞ and Y 2 ðtÞ can be used as the two bounding polynomials. The total computation time of each clipping process is O(n). In this case, the new method is called as the optimized one (see also Opt in Fig. 6).  Case 2: λ0  λ1 4 0. In this case, we ﬁrstly use QT3 to obtain the reference polynomial Y j ðtÞ ¼ y j ðϕj ðtÞÞ, where j¼3 or 4. Without loss of generality, suppose that j¼3. Let τ3 and τ4 be an upper bound and a lower ð4Þ bound of H 3 ðtÞ within ½0; 1, respectively. Note that pﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃ 51 17 þ 107 51 17  107 rðt  1Þðt  0:5Þt 2 r ¼ K2 K3 ¼  8192 8192 for all t A ½0; 1, we have the two bounding polynomials Y l;3 ðtÞ ¼ τ l;3 þY 3 ðtÞ r f ðtÞ r τ r;3 þ Y 3 ðtÞ ¼ Y r;3 ðtÞ;

ð25Þ

where τ r;3 ¼ maxf0; τ3 K 2 ; K 3 τ4 g and τ l;3 ¼ minf0; τ4 K 2 ; τ3 K 3 g. Similarly, suppose that u j;3 A ½0; 1 is a root of y 3 ðuÞ þ τ j;3 , then the corresponding root t ⋆ j;3 of Y j;3 ðtÞ can be computed by solving a quadratic equation ϕ3 ðtÞ  u j;3 ¼ 0, j ¼ l; r. Remark 3. Suppose the length ½a; b is small enough such that the value of γ0 tends to be 1, i.e., γ 0 ¼ 1 þ oð1Þ, where “oð1Þ” denotes inﬁnitively small. The reparameterization function ϕ3 ðtÞ can be a cubic or other polynomial satisfying

ϕ3 ð0Þ ¼ 0;

ϕ3 ð1Þ ¼ 1;

ϕ3 ð0:5Þ ¼ u m and ϕ03 ð0Þ ¼ γ 0 1 :

Let x^ 3 ðuÞ be a cubic polynomial determined by 3. The new quadratic clipping method for root ﬁnding Similar to the k-clipping method, each clipping process of the new quadratic clipping has three steps: (1) compute the reference polynomial; (2) estimate the bounding polynomials; and (3) bound the roots of f(t) by using the roots of the bounding polynomials. 0 Given a polynomial f(t) of degree n Z4, the interval ½a; b, f ðtÞ and ð4Þ f ðtÞ as well (by pre-computation once). We try to obtain a smaller subinterval by using one clipping process, which is as follows. For the sake of convenience, let the number pair ðT 1 ; T 2 Þ denotes T1 multiply/ division operations and T2 minus/add operations.

x^ 3 ð0Þ ¼ 0;

x^ 3 ð1Þ ¼ 1;

x^ 3 ðu m Þ ¼ 0:5

and

0 x^ 3 ð0Þ ¼ γ 0 :

It can be checked that ϕ3 ðx^ 3 ðuÞÞ ¼ u þ ðu  1Þðu  u m Þu2  oð1Þ and ϕ3 ðx^ 3 ðu j;3 ÞÞ  u j;3 ¼ oð1Þ  0. Thus, the root of ϕ3 ðtÞ  u j;3 ¼ 0 can be approximated by x^ 3 ðu j;3 Þ. 3.2. The third step for computing resulting subintervals The third step is to compute the subinterval ½a1 ; b1  D½0; 1 bounding the roots of f ðtÞ, then ½a þa1 ðb  aÞ; a þb1 ðb aÞ D ½a; b is the resulting subinterval bounding the root of f(t). Suppose that

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

the two bounding polynomials obtained in the ﬁrst two steps are Y l ðtÞ r f ðtÞ r Y r ðtÞ. We rearrange the roots of Y l ðtÞ and Y r ðtÞ and obtain that 0 ¼ t 0 ot 1 o⋯ o t d ¼ 1. Let Ωi be the subinterval ½t i ; t i þ 1 . If d ¼1, i.e., there is no root of Y l ðtÞ and Y r ðtÞ within ½0; 1, similar to that of the k-clipping method, we subdivide ½0; 1 and obtain two subintervals ½0; 0:5 and ½0:5; 1. Otherwise, if f ðtÞ Z Y l ðtÞ Z 0 or 0 ZY r ðtÞ Z f ðtÞ; 8 t A Ωi , note that the two inequalities are usually not satisﬁed at the same time, f ðtÞ has no root within Ωi, which means that Ωi can be removed. For the remaining case that both of Y l ðtÞ and Y r ðtÞ have roots within ½0; 1, note that f ðt l;i Þ Z Y l ðt l;i Þ ¼ 0 and f ðt r;i Þ r Y r ðt r;i Þ ¼ 0, from the middle ⋆ ⋆ value theorem, there exists t between t l;i and t r;i such that f ðt Þ ¼ 0, which means that the remaining subinterval has at least one root of f ðtÞ. Finally, the resulting subintervals are obtained in the above way. 3.3. Outline of the algorithm

Algorithm 1. Finding all of the roots of polynomial f(t) in ½a; b.

and

ε2 .

6. If the number of roots

κ is less than 2, put ½a; ða þ bÞ=2 and

½ða þ bÞ=2; b into Ω, and goto Step (2); Otherwise, remove the subinterval ½t i ; t i þ 1  which contains no roots of f(t), put the remaining subintervals into Ω, and goto Step (2). 7. End: Output Φ.

We show two examples to illustrate one clipping step of the new algorithm. Example 1 (A convex case). Given f ðtÞ ¼ ðt 1=3Þð2  tÞðt þ5Þ2 ; t A ½0; 1, which has a single root t¼1/3 in ½0; 1, see also Fig. 4(a). Note that λ0  λ1 ¼  682=3 o0, we ﬁrstly use QT1 to obtain   124u  29u2 950 þ3844u þ 1526u2 ; ; ðx 1 ðuÞ; y 1 ðuÞÞ ¼ 95 57

ϕ1 ðtÞ ¼ t þ 0:20124ðt  1ÞÞt:

In the new clipping method, the following clipping algorithm is iteratively executed until there are no roots of f(t) or the length of the subinterval is less than the given tolerance.

Input: The given polynomial f(t), interval ½a; b, tolerances

93

ε1

It can be veriﬁed that  0:0045 r H 1 ðtÞ ¼ f ðtÞ  Y 1 ðtÞ ¼ f ðtÞ  y 1 ðϕ1 ðtÞÞ r 0. The roots of the two bounding polynomials Y 1 ðtÞ and Y 1 ðtÞ  0:0045 form the resulting subinterval ½0:333304; 0:333400 of length 0.000096. In this case, the condition in Remark 2 is satisﬁed. From QT2, we construct

ϕ2 ðtÞ ¼ t þ 0:2012ðt  1Þt;

Output: The set Φ containing the roots of f(t) within ½a; b. 1. Compute the zeros of the control polygon of f(t), divide ½a; b into several subintervals by using these zeros, put all of the

y 2 ðtÞ ¼  50=3 þ 64:6788t  24:0122t 2 : Then Y i ðtÞ ¼ y i ðϕi ðtÞÞ bound f(t), 8 t A ½0; 1, i¼1, 2, see also Fig. 4(b). The corresponding subinterval is ½0:333304; 0:333336 of length 0.000032.

subintervals into Ω, and goto next step. 2. If Ω is empty, goto Step (7); otherwise, goto next step. 3. Pop up an element (i.e., an interval denoted as ½a1 ; b1 ) of Ω. If b1  a1 o ε1 , put ða1 þ b1 Þ=2 into Φ; otherwise, goto next step. 4. If jf ða1 Þj o ε2 , put a1 into Φ, and update the ½a1 ; b1  as ½a1 þ ε1 ; b1 . If jf ðb1 Þj o ε2 , put b1 into Φ, and update the ½a1 ; b1  as ½a1 ; b1  ε1 . Goto next step. 5. Compute the corresponding lower and upper bounds of f(t) by using Eq. (24), Eq. (25) or other similar ones, respectively. Solve the roots of the bounds and turn them into the values of parameter t to obtain a ¼ t 0 o t 1 o ⋯ o t κ ¼ b. Goto next step.

Example 2 (A concave case). Given f ðtÞ ¼ ððt  1=3Þ2 þ1Þðt  1=3Þ ðt þ 2Þ; t A ½0; 1, which has a single root t¼ 1/3 in ½0; 1, see Fig. 5 (a). In this case, λ0 ¼  4=3 and λ1 ¼ 8=3, so λ0  λ1 4 0. From QT3, note that point pm locates inside the triangle △p0 p1 r1 , we utilize Eq. (18) to obtain ðx 3 ðuÞ; y 3 ðuÞÞ ¼ ð  0:6254t 2 þ 1:6254t; 2:5682t 2  8:2719t þ 50=27Þ;

ϕ4 ðtÞ ¼ t 

0:2280ðt  1Þt t  1:3653

or ϕ4 ðtÞ ¼ t þ tðt  1Þð0:2695 þ 0:6086tÞ:

It can be veriﬁed that 0:413535 Z H 3 ðtÞ ¼ f ðtÞ  y 3 ðϕ4 ðtÞÞ Z 0:1878. The roots of the two bounding polynomials lead to the resulting subinterval ½0:3151; 0:4331 of length 0.1180. By the estimation method in Remark 3, the corresponding subinterval is ½0:3213; 0:4341 of length 0.1128.

H¯2 (t)

f (t)

t

t H¯1 (t)

Fig. 4. Example 1: (a) plots of the curve ðt; f ðtÞÞ in solid black, ðx 1 ðuÞ; y 1 ðuÞÞ in dashed blue and ðx 2 ðuÞ; y 2 ðuÞÞ in dashed red; (b) the error plots of H 1 ðtÞ in blue and H 2 ðtÞ in red. (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this article.)

94

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

Fig. 5. Example 2: (a) plots of the curve ðt; f ðtÞÞ in solid black and ðx 3 ðuÞ; y 3 ðuÞÞ in dashed red; (b) the error plots between f(t) and the reference polynomials in M2, M3 and Mn. (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this article.)

It shows that Mn has the least computation time among M2, M3 and Mn.

Table 1 Comparisons on length and approximation error. Method

Length Error

Example 1

Example 2

4.2. Analysis of convergence rate

M2

M3

Mn

M2

M3

Mn

0.0411 0.4976

0.0030 0.0142

0.000096 0.0017

0.1602 0.2142

0.0303 0.0142

0.1180 0.1945

For the sake of convenience, assume that t ⋆ is the unique root of f(t) in the interval ½a; b, m ¼ ða þ bÞ=2, and m r3 is the multiplicity of t ⋆ such that ðiÞ

f ðt ⋆ Þ ¼ 0; Let Mk and Mn be the k-clipping method and the new quadratic clipping method, respectively, k¼ 2, 3. The approximation error and the length of the resulting subinterval in the different methods are listed in Table 1. It shows that M3 achieves the best approximation error in Example 2 while Mn achieves the best approximation error in Example 1, both M3 and Mn have better approximation error than that of M2.

4. Computational complexity and convergence rate 4.1. Analysis of computational complexity Each clipping process in M2, M3 and Mn is similar to each other, which has three steps: (1) to ﬁnd a reference polynomial; (2) to compute two bounding polynomials and their roots within the given interval; and (3) to ﬁnd the resulting subintervals containing the roots of f(t). Suppose that the degree n Z4. We discuss the case λ0  λ1 o 0 by using QT1, and the other case λ0  λ1 4 0 can be done in a similar way. One needs to compute the values of p0, p1, d0 and d1 before determining y1 ðtÞ and ϕ1 ðtÞ. It can be veriﬁed that y1 ðtÞ and ϕ1 ðtÞ can be computed in ð4n þ 4; 4n þ 7Þ, and then the two bounding polynomials can be computed in ð2n2  14n þ 29; n2  7n þ 16Þ. A quadratic polynomial has at most two roots, and in the worst case, the roots of the two bounding polynomials can be computed in ð6; 18Þ and four “√” operations. When the condition in Remark 2 is satisﬁed, Y 1 ðtÞ and Y 2 ðtÞ can bound f(t) within the given subinterval, which leads to an optimized computation time ð6n þ 23; 6n þ 53Þ and ten “√” operations in the worst case. Fig. 6 plots the operation number curves of both multiply/division operation and add/minus operation in M2, M3, Mn and the optimized one mentioned in Remark 2 (denoted as “Opt”).

i ¼ 0; 1; …; m  1;

and

f

ðmÞ

ðt ⋆ Þ a 0:

ð26Þ

We have the following Lemma. Lemma 1. Assuming that h is small enough, we have   ðmÞ    1f ðt ⋆ Þ ⋆ m ðt  t Þ  r f ðtÞ; 8 t A ½a; b:   2 m!

ð27Þ

Proof. From Taylor expansion of f(t) at t ⋆ , there exists ξ5 ðtÞ A ½a; b such that f ðtÞ ¼

f

ðmÞ

ðm þ 1Þ

ðt ⋆ Þ f ðξ5 ðtÞÞ ðt  t ⋆ Þm þ ðt t ⋆ Þm þ 1 ; m! ðm þ 1Þ!

8 t A ½a; b:

From the assumption that h is small enough such that     f ðm þ 1Þ ðξ ðtÞÞ  f ðmÞ ðt ⋆ Þ    5 ⋆  ðt  t Þ r   ; 8 t A ½a; b:    2  mþ1

ð28Þ

ð29Þ

Combining Eq. (28) with Eq. (29), we obtain Eq. (27). Thus, we have completed the proof. □ Firstly, we discuss the case that λ0  λ1 o 0, e.g., the convex case. From QT1, we obtain y1 ðtÞ, the reparameterization function ϕ1 ðtÞ and then the reference polynomial Y 1 ðtÞ ¼ y1 ðϕ1 ðtÞÞ such that f ðiÞ ¼ Y 1 ðiÞ;

0

f ðiÞ ¼ Y 01 ðiÞ;

i ¼ a; b:

ð30Þ

From Eq. (30), there exists ξ6 ðtÞ A ½a; b such that H 1 ðtÞ ¼ f ðtÞ  Y 1 ðtÞ ¼

H 1ð4Þ ðξ6 ðtÞÞ ðt  aÞ2 ðt  bÞ2 : 4!

ð31Þ

Let η1 ¼ supa r t r b jðH 1ð4Þ ðtÞÞ=384j, combining with Eq. (31), we have that jH 1 ðtÞj r η1 ðb  aÞ4 ¼ η1 h : 4

ð32Þ

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

Nop

95

Nop

n

n

Fig. 6. Plots of operation numbers (Nop) in each clipping process of M2, M3 and Mn: (a) multiply/division operation; and (b) add/minus operation.

From Eq. (40), there exists ξ8 ðtÞ A ½a; b such that

Let t ⋆ j;l A ½a; b is the root of Y j;1 ðtÞ; j ¼ l; r. We have that jf ðt

⋆ ⋆ Þ  f ðt ⋆ j;l Þj ¼ jY j;1 ðt j;l Þ  f ðt j;l Þj ⋆ ⋆ ⋆ rjY 1 ðt ⋆ j;l Þ  f ðt j;l Þj þ jY j;1 ðt j;l Þ  Y 1 ðt j;l Þj r2η1 h : 4

ð33Þ

Combining Eq. (27) with Eq. (33), we have that   ðmÞ  1f ðt ⋆ Þ ⋆ 4 ⋆ m ðt j;l  t Þ  r 2η1 h ;   2 m!

4η1  m! jf

ðmÞ

ðt ⋆ Þj

Let η3 ¼ supa r t r b jðH ð4Þ 3 ðtÞÞ=384j, combining with Eq. (41), we have that ð42Þ

Let t ⋆ j;3 A ½a; b is the root of Y j;3 ðtÞ; j ¼ l; r. We have that ⋆ ⋆ jf ðt ⋆ Þ  f ðt ⋆ j;3 Þj ¼ jY j;3 ðt j;3 Þ  f ðt j;3 Þj

4

h :

ð34Þ

f ðmÞ ¼ Y 2 ðmÞ;

0

f ðmÞ ¼ Y 02 ðmÞ;

i ¼ a; b:

ð35Þ

From Eq. (35), there exists ξ7 ðtÞ A ½a; b such that H 2 ðtÞ ¼ f ðtÞ  Y 2 ðtÞ ¼

ð41Þ

4

Similarly, from QT2, we obtain y2 ðtÞ, the reparameterization function ϕ2 ðtÞ and then the reference polynomial Y 2 ðtÞ ¼ y2 ðϕ2 ðtÞÞ such that f ðiÞ ¼ Y 2 ðiÞ;

H ð4Þ 3 ðξ8 ðtÞÞ ðt aÞ2 ðt  bÞðt mÞ: 4!

jH 3 ðtÞj r η3 ðb  aÞ4 ¼ η3 h :

which leads to ⋆ m jðt ⋆ j;l  t Þ j r

H 3 ðtÞ ¼ f ðtÞ  Y 3 ðtÞ ¼

4!

ðt  aÞðt  bÞðt  mÞ2 :

ð36Þ

Let η2 ¼ supa r t r b jðH ð4Þ 2 ðtÞÞ=384j, combining with Eq. (36), we have that jH 2 ðtÞj r η2 ðb aÞ ¼ η2 h : 4

4

4

ð37Þ

ð43Þ

Combining Eq. (27) with Eq. (43), we have that   ðmÞ  1f ðt ⋆ Þ ⋆  4 ðt j;3  t ⋆ Þm  r2η3 h ;   2 m! which leads to ⋆ m jðt ⋆ j;3  t Þ j r

H ð4Þ 2 ð 7 ðtÞÞ

ξ

⋆ ⋆ ⋆ rjY 3 ðt ⋆ j;3 Þ  f ðt j;3 Þj þ jY j;3 ðt j;3 Þ  Y 3 ðt j;3 Þj r2η3 h :

4η3  m! jf

ðmÞ

ðt ⋆ Þj

4

h :

ð44Þ

From the above discussions, by combining Eqs. (34) and (39) with Eq. (44), we have the following theorem. Theorem 1. The new quadratic clipping method (by using QT1, QT2 or QT3) is capable of achieving a convergence rate of 4/m for solving the root t ⋆ of multiplicity m r 3.

Let t ⋆ j;2 A ½a; b is the root of Y j;2 ðtÞ; j ¼ l; r. We have that ⋆ ⋆ jf ðt ⋆ Þ  f ðt ⋆ j;2 Þj ¼ jY j;2 ðt j;2 Þ  f ðt j;2 Þj ⋆ ⋆ ⋆ r jY 2 ðt ⋆ j;2 Þ  f ðt j;2 Þj þ jY j;2 ðt j;2 Þ  Y 2 ðt j;2 Þj r2η2 h : 4

ð38Þ

Combining Eq. (27) with Eq. (38), we have that   ðmÞ  1f ðt ⋆ Þ ⋆  4 ðt j;2 t ⋆ Þm  r 2η2 h ;   2 m! which leads to ⋆ m jðt ⋆ j;2  t Þ j r

4η2  m! jf

ðmÞ

ðt ⋆ Þj

4

h :

ð39Þ

Secondly, we discuss the case λ0  λ1 40, e.g., a concave case. Without loss of generality, by using QT3, suppose that we obtain y3 ðtÞ, the reparameterization function ϕ3 ðtÞ and then the reference polynomial Y 3 ðtÞ ¼ y3 ðϕ3 ðtÞÞ such that f ðiÞ ¼ Y 3 ðiÞ;

0

f ðaÞ ¼ Y 03 ðaÞ;

i ¼ a; m; b:

ð40Þ

5. Further examples with comparison and discussions Let Mk and Mn be the k-clipping method in [1,15] and the new quadratic clipping method, k ¼2, 3. For the sake of convenience, when the condition in Remark 2 is satisﬁed, the optimized one is still denoted as Mn. We utilize the Maple software for testing more examples on PC with CPU 2.2 G and 4 G memory. In this section, the average computation time of the ﬁrst clipping process is tested by using 17 valid digits after the point, which is similar to the ﬂoat number, and the corresponding unit is millisecond. In the following tables, αi e  β i denotes the real number Li ¼ αi  10  βi , which is also the length of the resulting subinterval of the i-th clipping process. And then, the corresponding convergence rate can be simply approximated by β i þ 1 =β i . From Theorem 1, we can deduce the distribution of the roots from the information of the intervals N f½ai ; bi gi ¼ 0 . For example, suppose that Li ¼ bi  ai ¼ αi e β i , where

96

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

Table 2 (Example 3). Length Li of the resulting subinterval at different N. f(t)

f 1 ðk ¼ 1Þ

f 2 ðk ¼ 1Þ

f 3 ðk ¼ 1Þ

N

M2

M3

Mn

M2

M3

Mn

M2

M3

Mn

1 2 3 4 5 cr Time

4e  2 2e  6 6e  19 1e  56 8e  170 3 2.247

3e  3 2e  13 1e  53 4e  215 1e  860 4 3.994

3e  5 1e  23 2e  96 4e  388 6e  1555 4 0.2932

1e  2 2e  8 9e  26 5e  78 9e  235 3 2.652

9e  3 8e  11 5e  43 7e  172 3e  687 4 4.336

1e  2 2e  10 3e  41 4e  165 2e  660 4 0.327

1e  2 2e  7 5e  22 7e  66 2e  197 3 3.916

1e  2 1e  9 3e  38 2e  152 8e  596 4 5.772

3e  2 3e  8 7e  32 2e  127 6e  509 4 0.374

f(t)

f 4 ðk ¼ 2Þ

N

M2

M3

Mn

M2

M3

Mn

M2

M3

Mn

1 2 3 4 5 cr Time

2e  1 2e  2 7e  4 5e  6 2e  9 1.5 2.543

1e  1 1e  3 2e  7 4e  15 2e  30 2 4.118

3e  2 9e  6 3e  14 9e  35 3e  84 2.4 0.3541

2e  1 2e  2 1e  3 8e  6 6e  9 1.5 3.623

7e  2 4e  4 1e  8 1e  17 2e  35 2 4.196

5e  2 8e  5 1e  11 4e  28 7e  68 2.4 0.4025

3e  1 6e  2 5e  3 1e  4 6e  7 1.5 4.166

2e  2 2e  4 5e  9 2e  18 5e  37 2 5.679

7e  2 3e  4 4e  10 4e  24 5e  58 2.4 0.4337

f 5 ðk ¼ 2Þ

f 6 ðk ¼ 2Þ

Table 3 (Example 4) Length Li of the resulting subinterval at different N. f(t)

f 7 ðk ¼ 2Þ

f 8 ðk ¼ 2Þ

N

M2

M3

Mn

M2

M3

Mn

M2

M3

Mn

1 2 3 4 5 6 cr Time

2e  1 2e  2 8e  4 5e  6 3e  9 3e  14 1.5 3.104

1e  1 1e  3 2e  7 5e  15 3e  30 8e  61 2 5.070

3e  2 9e  6 3e  14 1e  34 3e  84 9e  204 2.4 0.336

1e  1 7e  3 1e  4 4e  7 7e  11 1e  16 1.5 3.542

8e  2 4e  4 1e  8 2e  17 2e  35 4e  71 2 5.211

6e  2 8e  5 1e  11 4e  28 8e  68 9e  164 2.4 0.363

1e  1 2e  2 7e  4 7e  6 7e  9 2e  13 1.5 5.023

6e  2 3e  4 6e  9 3e  18 5e  37 2e  74 2 6.474

7e  2 3e  4 3e  10 2e  24 8e  59 9e  142 2.4 0.429

the interval ½ai ; bi  is obtained from the i-th iteration step, 1 r αi o 10 and βi is a positive integer. If the closest integer number to βi þ 1 =βi is 4 for all i Z N, where N is a positive integer, then there is a single root in the interval ½ai ; bi ; if the closest integer number to β i þ 1 =β i is 2 for all i Z N, then there is a double root in the interval ½ai ; bi ; if the closest integer number to β i þ 1 =β i is 2 for N 1 Z iZ N, but becomes 1 for i Z N1 , there are possibly two (near double) roots in the interval ½ai ; bi . Example 3. Given six polynomials f i ðtÞ; i ¼ 1; …; 6, within ½0; 1, 8 > f ðtÞ ¼ ðt  1=3Þð2  tÞðt þ 5Þ2 ; > > 1 > > 3 4 > > > f 2 ðtÞ ¼ ðt  1=3Þð2  tÞ ðt þ5Þ ; > > 5 < f ðtÞ ¼ ðt  1=3Þð2  tÞ ðt þ5Þ10 ; 3 2 > > f 4 ðtÞ ¼ ðt  1=3Þ ð4  tÞðt þ7Þ; > > > > > f 5 ðtÞ ¼ ðt  1=3Þ2 ðt 4Þ3 ðt þ 5Þ2 ðt þ7Þ; > > > : f ðtÞ ¼ ðt  1=3Þ2 ðt 4Þ7 ðt þ 5Þ6 ðt þ7Þ; 6

where the corresponding curve Ci ðtÞ ¼ ðt; f i ðtÞÞ is convex nearby the root t ⋆ ¼ 1=3. Table 2 shows computation time and length of the resulting subinterval from M2, M3 and Mn. In Table 2, k means the multiplicity of the root, and cr denotes the convergence rate. It shows that both M3 and Mn achieve convergence rate 4/k while M2 is of 3/k. M3 obtains the smallest length in cases f 2 ðtÞ and f 3 ðtÞ, while Mn obtains the smallest length in other four cases. The optimized method can work for these six cases, and the average computation time of Mn is about 10% of those of M2 and M3.

f 9 ðk ¼ 2Þ

Example 4. We apply the algorithm to the three polynomials f 7 ðtÞ ¼ ðt  1=3Þ2 ð4  tÞðt þ7Þ, f 8 ðtÞ ¼ ðt  1=3Þ2 ðt 4Þ3 ðt þ 5Þ2 ðt þ7Þ and f 9 ðtÞ ¼ ðt  1=3Þ2 ðt  4Þ7 ðt þ5Þ6 ðt þ 7Þ, in order to compute the double root 1/3 in the interval ½0; 1. Table 3 shows the resulting lengths of the subintervals and the computation time in M 2 , M 3 and M n . These three cases show that M n has the smallest length and is the fastest among M 2 , M 3 and M n . The convergence rates of M 2 , M 3 and M n tend to be 3/2, 2 and 2.4, which shows that M n has the best performance on the convergence rate. As mentioned in , the k-clipping method is usually used for computing the cases where there are at most k real roots (counting the multiplicity) of the given polynomial f(t) in an interval. Example 6 shows a triple root case, in which f ″ðt ⋆ Þ ¼ 0 and the curve ðt; f ðtÞÞ is concave nearby the root t ⋆ . Example 5. We apply M2, M3 and Mn to the polynomial f ðtÞ ¼ ðt  1=3Þ3 ðt  3Þ, in order to compute the triple root t ⋆ ¼ 1=3 in the interval ½0; 1. As shown in Fig. 7(a), M2, M3 and Mn work for this case, and M3 has the best approximation effect, while Mn has a slightly better approximation effect than that of M2. We have also tested the numerical robustness of the planar quadratic clipping method (PQCM). Similar to the quadratic clipping method, PQCM needs to solve quadratic polynomial equations, which is trivial and stable to solve. Similar to the cubic clipping method in , we have also used PQCM to compute the

X.-D. Chen, W. Ma / Computers & Graphics 46 (2015) 89–98

97

Fig. 7. Example 5 (a triple root): (a) plots of the curve ðt; f ðtÞÞ in solid black and the reference curves in M2, M3 and Mn, respectively; (b) the error plots between f(t) and the corresponding reference polynomials.

roots of the Wilkinson polynomial 20

WðxÞ ¼ ∏ ðx  iÞ

no. CityU 118512), Defense Industrial Technology Development Program, the National Science Foundation of China (61003194, 61370166 and 61379072).

i¼1

with Bernstein–Bézier representation within the domain interval ½0; 25 with the same error tolerance ε1 ¼ 10  3 . In this case, PQCM can also generate 20 intervals containing the roots whose lengths are less than ε1. In other word, it can achieve a comparable result to that of the cubic clipping method in .