- Email: [email protected]

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

This paper presents a new quadratic clipping method for computing a root of a polynomial f(t) of degree n within an interval. Different from the traditional one in R1 space, it derives three quadratic curves in R2 space for approximating ðt; f ðtÞÞ instead, which leads to a higher approximation order. Two bounding polynomials are then computed in Oðn2 Þ for bounding the roots of f(t) within the interval. The new clipping method achieves a convergence rate of 4 for a single root, compared with that of 3 from traditional method using quadratic polynomial approximation in R1 space. When f(t) is convex within the interval, the two bounding polynomials are able to be directly constructed without error estimation, which leads to computational complexity O(n). Numerical examples show the approximation effect and efﬁciency of the new method. The method is particularly useful for the ﬂy computation in many geometry processing and graphics rendering applications. & 2014 Elsevier Ltd. All rights reserved.

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 [3], collision detection [4,14], and bisectors/ medial axes computation [8]. 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 [16]. 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 [17]. There are also several clipping algorithms developed for ﬁnding the roots of polynomials [1,15,17,23]. In 1990, Sederberg and Nishita

n

Corresponding author. E-mail addresses: [email protected] (X.-D. Chen), [email protected] (W. Ma). http://dx.doi.org/10.1016/j.cag.2014.09.014 0097-8493/& 2014 Elsevier Ltd. All rights reserved.

proposed a Bézier clipping technique, which is to identify regions of the polynomial which do not include the part of the roots by using the convex hull property of Bézier curves [23]. The Bézier clipping method was proved to be quadratically convergent in [24]. Mørken and Reimers presented a linear approximation method [17], in which the given B-Spline (or Bézier) curve is approximated by its control polygon, which achieves a convergence rate of 2. Bartoň and Jüttler presented a quadratic clipping method for computing all the roots of a univariate polynomial equation within an interval [1]. In the quadratic clipping method in [1], the original univariate polynomial is ﬁrstly approximated by a quadratic polynomial based on degree reduction. Then two quadratic polynomials bounding the original polynomial within an interval are estimated, and the roots of the original polynomial are also bounded by the roots of the two bounding polynomials. The quadratic clipping method achieves a convergence rate of 3 for a single root and a super-linear rate of 3/2 for double roots, which is faster than the Bézier clipping method. Later, Liu et al. extended the quadratic clipping technique and proposed a cubic clipping method, which uses a cubic polynomial to approximate the original polynomial based on degree reduction and achieves a convergence rate of 4 for single roots [15]. This paper presents a planar quadratic clipping method for computing a root of a polynomial within an interval, which can achieve an optimal convergence rate of 4 for a single root by using quadratic polynomial curves. Suppose that the given polynomial f (t) is of degree n. Instead of approximating f(t) in R1 space by using degree reduction, we use a planar quadratic polynomial curve to

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Þ

which are used to approximate the curve CðtÞ ¼ ðt; f ðtÞÞ within the interval ½0; 1, where B20 ðuÞ ¼ ð1 uÞ2 , B21 ðuÞ ¼ 2ð1 uÞu and B22 ðuÞ ¼ u2 are Bernstein polynomials. 2.1. Outline of the idea In the previous quadratic clipping method, it utilizes a quadratic polynomial y(t) to approximate f(t), while in the planar quadratic clipping method, it utilizes a quadratic polynomial curve ðx2 ðuÞ; y2 ðuÞÞ to approximate ðt; f ðtÞÞ in R2 space instead, which needs to estimate the mapping function between parameters u and t. In principle, the planar quadratic clipping method is to ﬁnd a quadratic reparameterization function ϕðtÞ ¼ t þ ρ1 ðt 1Þt and a quadratic polynomial yðtÞ ¼ ρ2 þ ρ3 t þ ρ4 t 2 such that yðϕðtÞÞ obtains the approximation 0 order of 4 to f(t), where ϕ ðtÞ Z 0, for all t A ½0; 1. Note that there are four unknowns ρi ; i ¼ 1; 2; 3; 4 in yðϕðtÞÞ, we can introduce four constraints such that yðϕðtÞÞ interpolates either two points and their two derivatives (see also QT1 in Section 2.2), or three points and one derivative of f(t) (see also QT2 in Section 2.3 and QT3 in Section 2.4). From [6, Theorem 3.5.1, p. 67, Chapter 3.5], yðϕðtÞÞ can obtain both the approximation order of 4 and the convergence rate of 4, which is better than that of previous quadratic clipping method 3. Generally speaking, QT1 and QT2 are used for convex cases that ðt; f ðtÞÞ has no inﬂexion point, while QT3 is applied for other cases that QT1 and QT2 fail. In principle, both QT1 and QT2 can construct two polynomials which can be used to bound f(t), and the corre-

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Þ 1f ð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Þ 1f ð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Þ 1f ð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Þ 1f ð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 [15], 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 [15], 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 [15].

6. Conclusions This paper presents a quadratic tangent method for ﬁnding the roots of a polynomial within an interval, which utilizes both the tangent idea and the reparameterization technique in R2 space instead of R1 space in the traditional k-clipping method. It is able to achieve a convergence rate of 4 for a single root which is the same as that of a higher order cubic clipping method in [15]. On the other hand, the new method tries to avoid the timeconsuming computation of the bounding polynomials, which can lead to a much higher computation efﬁciency. Numerical examples show that the convergence rate of the quadratic tangent method can be equal to that of the cubic clipping method which is higher than that of the traditional quadratic clipping method, while the new tangent method is able to perform much faster than those of previous methods in [1,15]. s for future work, several areas can be further addressed. When the parameter interval is too large, though the bounding polynomials can bound f(t), the bounding polynomials may fail to have roots in the interval. In this case, both the methods in [1,15] and the new quadratic tangent method may fail to obtain a much smaller interval. Note that a cubic curve can achieve a better approximation effect than a quadratic curve in most cases. In our future work, we will discuss the cubic tangent method which is promising to achieve higher convergence rate and better approximation effect than those of the quadratic tangent one.

Acknowledgments The research was partially supported by the Research Grants Council of Hong Kong Special Administrative Region, China (Grant

References [1] Bartǒ n M, Jüttler B. Computing roots of polynomials by quadratic clipping. Comput Aided Geom Des 2007;24:125–41. [2] Borwein P, Erdelyi T. Polynomials and polynomial inequalities. New York: Springer-Verlag; 1995. [3] Chen XD, Yong JH, Wang GZ, Paul JC, Xu G. Computing the minimum distance between a point and a NURBS curve. Comput-Aided Des 2008;40(10– 11):1051–4. [4] Choi YK, Wang WP, Liu Y, Kim MS. Continuous collision detection for two moving elliptic disks. IEEE Trans Robot 2006;22(2):213–24. [5] Collins GE, Akritas AG. Polynomial real root isolation using Descartes's rule of signs. In: Proceedings of the third ACM symposium on symbolic and algebraic computation. New York, US; 1976. p. 272–5. [6] Davis PJ. Interpolation and approximation. New York: Dover Publications; 1975. [7] Efremov A, Havran V, Seidel H-P. Robust and numerically stable Bézier clipping method for ray tracing NURBS surfaces. In: The 21st Spring conference on computer graphics. New York: ACM Press; 2005. p. 127–35. [8] Elber G, Kim M-S. Geometric constraint solver using multivariate rational spline functions. In: The sixth ACM/IEEE symposium on solid modeling and applications. Ann Arbor, Michigan; 2001. p. 1–10. [9] Floater M. High order approximation of rational curves by polynomial curves. Comput Aided Geom Des 2006;23(8):621–8. [10] Hook DG, McAree PR. Using Sturm sequences to bracket real roots of polynomial equations. In: Glassner AS, editor. Graphics gems I; 1990. p. 416–22. [11] Isaacson E, Keller HB. Analysis of numerical methods. New York: Wiley; 1966. [12] Jaklič G, Kozak J, Krajnc M, Žagar E. On geometric interpolation by planar parametric polynomial curves. Math Comput 2007;76(260):1981–93. [13] Jüttler B. The dual basis functions of the Bernstein polynomials. Adv Comput Math 1998;8:345–52. [14] Lin M, Gottschalk S. Collision detection between geometric models: a survey. In: Proceedings of the IMA conference on mathematics of surfaces. Birmingham, UK; 1998. p. 37–56. [15] Liu LG, Zhang L, Lin BB, Wang GJ. Fast approach for computing roots of polynomials using cubic clipping. Comput Aided Geom Des 2009;26:547–59. [16] McNamee JM. Bibliographies on roots of polynomials. J Comp Appl Math 1993–2002;47:391–4 p. 110, 305–306; 142, 433–434, available at: 〈http:// www1.elsevier.com/homepage/sac/cam/mcnamee/〉. [17] Mørken K, Reimers M. An unconditionally convergent method for computing zeros of splines and polynomials[J]. Math Comput 2007;76(258):845–65. [18] Nishita T, Sederberg TW, Kakimoto M. Ray tracing trimmed rational surface patches. In: Proceedings of the Siggraph. ACM; 1990. p. 337–45. [19] Patrikalakis N.M., Maekawa T ,Ko K.H, Mukundan H. Surface to surface intersections. Comput.-Aided Des. Appl. 2004, 1 (1-4):449-458. [20] Reuter M, Mikkelsen T, Sherbrooke E, Maekawa T, Patrikalakis N. Solving nonlinear polynomial systems in the barycentric Bernstein basis. Vis Comput 2007;24(3):187–200.

98

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

[21] Rouillier F, Zimmermann P. Efﬁcient isolation of polynomial's real roots. J Comput Appl Math 2004;162(1):33–50. [22] Scherer K. Parametric polynomial curves of local approximation of order 8. In: Cohen A, Rabut C, Schumaker LL, editors. Curve and surface ﬁtting – Saint Malo 1999. Nashville: Vanderbilt University Press; 2000. p. 375–84.

[23] Sederberg TW, Nishita T. Curve intersection using Bézier clipping. ComputAided Des 1990;22(9):538–49. [24] Schulz C. Bézier clipping is quadratically convergent. Comput Aided Geom Des 2009;26(1):61–74.