- Email: [email protected]

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Rational cubic clipping with linear complexity for computing roots of polynomials Xiao-diao Chen a,b,∗, Weiyin Ma b a b

School of Computer, Hangzhou Dianzi University, Hangzhou, PR China Department of MBE, City University of Hong Kong, Hongkong

a r t i c l e

i n f o

Keywords: Approximation order Fast cubic clipping method Root-ﬁnding Convergence rate Linear computational complexity

a b s t r a c t Many problems in computer aided geometric design and computer graphics can be turned into a root-ﬁnding problem of polynomial equations. Among various clipping methods, the ones based on the Bernstein–Bézier form have good numerical stability. One of such clipping methods is the k-clipping method, where k = 2, 3 and often called a cubic clipping method when k = 3. It utilizes O(n2 ) time to ﬁnd two polynomials of degree k bounding the given polynomial f(t) of degree n, and achieves a convergence rate of k + 1 for a single root. The roots of the bounding polynomials of degree k are then used for bounding the roots of f(t). This paper presents a rational cubic clipping method for ﬁnding two bounding cubics within O(n) time, which can achieve a higher convergence rate 5 than that of 4 of the previous cubic clipping method. When the bounding cubics are obtained, the remaining operations are the same as those of previous cubic clipping method. Numerical examples show the eﬃciency and the convergence rate of the new method. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Many problems in computer aided geometric design and computer graphics can be turned into a root-ﬁnding problem of polynomial equations, such as robotics [15], computer aided design and manufacturing [8,19], including curve/surface intersection [7,13,18,19], point projection [2], collision detection [5,12], and bisectors/medial axes computation [8]. In principle, a system of polynomial equations of multiple variables can also be turned into a univariate polynomial equation by using the resultant theory. This paper discusses the root-ﬁnding problem of a univaraite polynomial equation within an interval. In many references, the given polynomial f(t) of degree n is written into its power series [10,14,17,20]. However, the Bernstein– Bézier form of f(t) has good stability with respect to perturbations of the coeﬃcients [9,11]. Several methods based on the Bernstein–Bézier form are developed [1,13,16,21], which are proved to achieve good numerical stability. Note that the number of zeros of a Bézier function is less or equal to that of its control polygon, the method in [16] utilizes the corresponding control polygon to approximation f(t), in which the zeros of the control polygon is used to approximate the zeros of f(t) from one side. The method in [16] achieves the convergence rate 2 for a simple root. In principle, a B-spline (or 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, i.e., the roots are bounded from two sides. The corresponding approximation order of the convex hull is 2. Comparing with the method in [16], the k-clipping method in [1,13] bounds the zeros of f(t) by using the zeros of two bounding polynomials of degree k, which ∗

Corresponding author at: School of Computer, Hangzhou Dianzi University, Hangzhou, PR China. Tel.: +86 571 86915140. E-mail addresses: [email protected], [email protected] (X.-d. Chen), [email protected] (W. Ma).

http://dx.doi.org/10.1016/j.amc.2015.10.054 0096-3003/© 2015 Elsevier Inc. All rights reserved.

1052

X.-d. Chen, W. Ma / Applied Mathematics and Computation 273 (2016) 1051–1058

achieves a higher approximation order k + 1, where k is usually set as 2 or 3. Recently, Chen et. al presented a planar quadratic clipping method by using rational quadratic polynomials to approximate (t, f(t)) in R2 space, which achieves the convergence rate 4 for a simple root [3]. In principle, for a case that (t, f(t)) is convex within the given interval, the rational quadratic clipping method in [4] can achieve the convergence rate 5 for a single root. However, in case when the curve (t, f(t)) is not convex within [a, b], the denominators of the rational quadratic polynomials for bounding f(t) may have one or more zeros within [a, b], which leads to a bad approximation effect between f(t) and its bounding polynomials. In principle, the above clipping methods have mainly two steps for each clipping process, i.e., one is to compute the two bounding polynomials whose computation time is dominant, and the other is to obtain the resulting subintervals based on the roots of the bounding polynomials within O(n) time. In the worst case, the above methods need O(n2 ) time to compute the two bounding polynomials during each clipping process, and the total computational complexity is O(n2 ). This paper presents a fast cubic clipping method of linear computational complexity. The main contribution is to directly construct the two bounding cubics within O(n) time. It is proved that the bounding cubics achieve the approximation order 4 to f(t), which is the same as that of previous cubic clipping method in [13]. Once the bounding cubics are obtained, the remaining operations are the same as that of the cubic clipping method in [13]. Numerical examples conﬁrms higher computational eﬃciency of the new method. The remainder of this paper is organized as follows. Section 2 provides an outline of the proposed clipping method and discusses two key issues of a general clipping method. Section 3 presents the method for constructing two bounding cubics within a linear time. Section 4 illustrates the algorithm with several numerical examples and other related discussions. Conclusions are drawn at the end of this paper. 2. Clipping methods 2.1. Outline of clipping methods The basic idea of clipping methods is to iteratively clip the parts of an interval containing no roots of the given polynomial f(t). In each clipping process, there are mainly two steps: (1) The ﬁrst step is to compute two bounding polynomials, which is a key issue of a clipping method. (2) The second step is to divide the given interval into several subintervals by using the roots of the bounding polynomials, and clip the subintervals containing no roots of f(t). In principle, different clipping methods may obtain different bounding polynomials. Once two bounding polynomials are obtained such that g1 (t) ≤ f(t) ≤ g2 (t). The remaining work is to compute the roots of gi (t), to rearrange the roots in order such that a = t0 < t1 < t2 < · · · < tl = b, and to remove the subintervals containing no roots of f(t). Solving a polynomial equation of degree 3 or less seems to be trivial. The details of judging whether or not a subinterval i contains one or more roots of f(t) are as follows. It can be divided into three cases: • Case 1: If g1 (t) > 0 for all t ∈ i , i can be removed. From the assumption, we have that f(t) ≥ g1 (t) > 0, for all t ∈ i , which means that i contains no root of f(t) and can thus be removed. • Case 2: If g2 (t) < 0 for all t ∈ i , i can be removed. From the assumption, we have that f(t) ≤ g2 (t) < 0, for all t ∈ i , which means that i contains no root of f(t) and can also be removed. • Case 3: g1 (t) ≤ 0 and g2 (t) ≥ 0, i is kept for further dealing. 2.2. The analysis of convergence rate The analysis of convergence rate is another key issue of a clipping method, which depends on the approximation order of the bounding polynomials. We have the following theorem. Theorem 1. Suppose that the two bounding polynomials gi (t), i = 1, 2, achieve the approximation order m to f(t) within interval [a, . b] whose length is small enough for satisfying Eq. (2), then the corresponding convergence rate for a root t of multiplicity k is m k Proof. From the assumption that t has multiplicity k, we have that

f ( j) (t ) = 0,

for

j = 0, 1, . . . , k − 1,

and

f (k) (t ) = 0.

(1)

Combining Taylor’s expansion with Eq. (1), there exists a small enough η > 0 and ξ 1 (t) such that

f (k) (t ) f (k+1) (ξ1 (t )) k k+1 | f (t ) − f (t )| = (t − t ) + (t − t ) k! (k + 1)! 1 f (k) (t ) k > (t − t ) , ∀ t ∈ [t − η, t + η].

2

k!

(2)

X.-d. Chen, W. Ma / Applied Mathematics and Computation 273 (2016) 1051–1058

1053

Let ti be a root of gi (t). From Eq. (2), we have that

1 f (k) (t ) k | f ( ) − f (t )| > (ti − t ) , i = 1, 2. 2 k! ti

(3)

From the assumption that gi (t) achieves the approximation order m to f(t) within [a, b], there exists a constant ci such that

| f (t ) − gi (t )| < ci (b − a)m , for all t ∈ [a, b].

(4)

Combining Eq. (3) and Eq. (4), we have that

1 f (k) (t ) k < | f (ti ) − f (t )| = | f (ti )| ( t − t ) i 2 k! = | f (ti ) − gi (ti )| < ci (b − a)m .

From Eq. (5), we obtain that

| −t | < ti

2ci k! f (k) (t )

1k

(5)

(b − a) k , m

(6)

which means that the corresponding convergence rate is

m . k

For the sake of convenience, we introduce Theorem 3.5.1 in Page 67, Chapter 3.5 of [6] as follows. Theorem 2. Let w0 , w1 , . . . , wr be r + 1 distinct points in [a, b], and n0 , . . . , nr be r + 1 integers ≥ 0. Let N = n0 + · · · + nr + r. Suppose that g(t) is a polynomial of degree N such that

g(i) (wκ ) = f (i) (wκ ),

i = 0, . . . , nκ ,

κ = 0, . . . , r.

Then there exists ξ 2 (t) ∈ [a, b] such that

f (t ) − g(t ) =

f (N+1) (ξ2 (t )) (t − wi )ni . (N + 1)! r

i=0

In the following two sections, we provide two methods of linear complexity for estimating the two bounding polynomials, which leads to convergence rates 4/k and 6/k for a root of multiplicity k. 3. Constructing two bounding rational cubics within a linear time The idea of this section is to ﬁrstly construct two quartics such that Yl (t) ≤ f(t) ≤ Yr (t), and then to construct two rational cubics such that Rl (t) ≤ Yl (t) ≤ f(t) ≤ Yr (t) ≤ Rr (t). Since Yl (t) and Yr (t) achieve the approximation order 5 to approximate f(t), Rl (t) and Rr (t) achieve the approximation order 5 to approximate Yl (t) and Yr (t), respectively, so both Rl (t) and Rr (t) achieve the approximation order 5 to approximate f(t). From Theorem 1, the corresponding convergence rate is 5/k, wher k is the multiplicity of the root t of f(t). Suppose that the given polynomial f(t) is written as ni=0 ci t i , t ∈ [a, b]. If a < 0, let g(u) = f (u + a), and we compute the roots of g(u) within [0, b − a] instead, and ui + a is a root of f(t) within [a, b], where ui is a root of g(u) within [0, b − a]. Without loss of generality, suppose that a ≥ 0. For the sake of convenience, let Bm (t ) = Cim (b−t()b−a)(mt−a) be a Bernstein basis function of degree m i m m mapping to the interval [a, b], Yi,m (t ) = j=0 di,m B j (t ) be a polynomial of m, i = 1, 2, . . . , 4, t0 = a, t1 = (a + b)/2 and t2 = b. m−i

i

3.1. Computing two bounding quartics for the case that f(5) (t) ≥ 0, ∀t ∈ [a, b] This section assumes that f(5) (t) ≥ 0, ∀t ∈ [a, b], and introduces two quartic polynomials Y1, 4 (t) and Y2, 4 (t) which are determined by

Y1,4 (ti ) = f (ti ), Y2,4 (ti ) = f (ti ),

(t ) = f (t ), Y1,4 κ κ (t ) = f (t ), Y2,4 κ κ

i = 0, 1, 2, i = 0, 1, 2,

and and

κ = 0, 1; κ = 1, 2.

(7)

Combining Eq. (7) with Theorem 2, we obtain that

⎧ f (5)(t ) ⎪ ⎨ f (t ) − Y1,4 (t ) = (t − t0 )2 (t − t1 )2 (t − t2 ) ≤ 0; 5!

⎪ ⎩ f (t ) − Y (t ) = f (5)(t ) (t − t )(t − t )2 (t − t )2 ≥ 0. 2,4 0 1 2 5!

Note that

f(5) (t)

≥ 0, from Eq. (8), we have that

Y2,4 (t ) ≤ f (t ) ≤ Y1,4 (t ),

for all t ∈ [a, b],

which means that Y1,4 (t) and Y2,4 (t) bound f(t) within [a, b].

(8)

1054

X.-d. Chen, W. Ma / Applied Mathematics and Computation 273 (2016) 1051–1058

Remark 1. Suppose that f(4) (t) ≥ 0, ∀t ∈ [a, b], Y1,3 (t) and Y2,3 (t) are two cubics satisfying that

(t ) = f (t ), Y1,3 1 1 (t ) = f (t ), Y2,3 i i

Y1,3 (ti ) = f (ti ), Y2,3 (ti ) = f (ti ),

i = 0, 1, 2; i = 0, 2;

(9)

then we have that Y2, 3 (t) ≤ f(t) ≤ Y1, 3 (t). The above method is a cubic clipping method with linear computational complexity. 3.2. Computing two bounding quartics for a general case Suppose that + = {0 ≤ i ≤ n | ci ≥ 0} and − = {0 ≤ i ≤ n | ci < 0}. Let

f+ (t ) =

ci t i

f− (t ) =

and

i∈+

ci t i ,

i∈−

(κ)

(κ)

It can be veriﬁed that f+ (t ) ≥ 0 and f− (t ) ≤ 0, for all t ≥ 0, where κ = 0, 1, 2, . . .. From the method in Section 3.1, we can directly construct four quartics Yi,4 , i = 1, 2, 3, 4 such that

Y2,4 (t ) ≤ f+ (t ) ≤ Y1,4 (t )

and Y4,4 (t ) ≤ f− (t ) ≤ Y3,4 (t ).

(10)

From Eq. (10), we have that

Yl (t ) = Y2,4 (t ) + Y4,4 (t ) ≤ f+ (t ) + f− (t ) = f (t ) ≤ Y1,4 (t ) + Y3,4 (t ) = Yr (t ),

(11)

which means that two quartics Yl (t) and Yr (t) bound f(t) within [a, b] for a general case. Let M1 = sup

| f+(5) (t )| 5!

a≤t≤b

and M2 = sup a≤t≤b

| f−(5) (t )| 5!

| f (t ) − Y j (t )| ≤ (M1 + M2 )(b − a)5 ,

. From Eq. (11), we have that

j = l, r,

(12)

which means that both Yl (t) and Yr (t) achieve the approximation order 5 to f(t) within [a, b]. (5)

(5)

Remark 2. There may be other better methods for dividing f(t) into two parts fneg (t) and fpos (t) such that fneq (t ) ≤ 0 and f pos (t ) ≥

0, for all t ∈ [a, b], which can be used to further improve the approximation effect. Speciﬁcally, if ( − 1) j f (5) (t ) ≤ 0 for all t ∈ [a, b], where j = 1 or j = 2, there is no need for dividing f(t), which will lead to an improved result. 3.3. Constructing two rational cubics for bounding a quartic Given a quartic Y j (t ) =

Rκ , j (t ) =

4

3

i=0

3 i=0 ri,κ , j Bi

t + r4,κ , j

di, j t i within [a, b], j = l, r, we construct two rational cubics

(t )

=

Xκ , j (t ) , Wκ , j (t )

κ = 1, 2, 3, 4, and j = l, r,

(13)

to bound Y(t). Let Hκ , j (t ) = Wκ , j (t )Y j (t ) − Xκ , j (t ). The values of ri, κ , j can be determined by

H1, j (ti ) = 0, H2, j (ti ) = 0,

(t ) = 0, H1, j κ (t ) = 0, H2, j κ

i = 0, 1, 2, i = 0, 1, 2,

κ = 0, 1; κ = 1, 2;

and and

j = l, r.

(14)

Combining Eq. (14) with Theorem 2, we obtain that

H1, j (t ) = d4, j (t − t0 )2 (t − t1 )2 (t − t2 ); H2, j (t ) = d4, j (t − t0 )(t − t1 )2 (t − t2 )2 ;

j = l, r.

(15)

From Eq. (15), we have that H1, j (t) · H2, j (t) ≤ 0. If W1, j (t) or W2, j (t) has a root within [a, b], we use the cubic clipping method of linear computational complexity mentioned in Remark 1 to ﬁnd the corresponding bounding polynomials of f(t) instead. Otherwise, both W1, j (t) and W2, j (t) have no roots within [a, b], i.e., W1, j (t) · W2, j (t) > 0, ∀t ∈ [a, b], we obtain that

H1, j (t ) H2, j (t ) · = (Y j (t ) − R1, j (t )) · (Y j (t ) − R2, j (t )) ≤ 0, W1, j (t ) W2, j (t )

(16)

which means that R1, j (t) and R2, j (t) bound Yj (t). In this case, without loss of generality, suppose that Wi, j (t) > 0. Let wi, j = inf Wi, j (t ) > 0. We have that

a≤t≤b

|Y j (t ) − Ri, j (t )| ≤

d4, j (b − a)5 , wi, j

i = 1, 2,

and

j = l, r.

(17)

By using the above method, we can construct two rational cubics Rl (t ) = min{R1,l (t ), R2,l (t )} and Rr (t ) = max{R1,r (t ), R2,r (t )} such that Rl (t) ≤ Yl (t) ≤ f(t) ≤ Yr (t) ≤ Rr (t). Combining Eq. (12) with Eq. (17), we have that

X.-d. Chen, W. Ma / Applied Mathematics and Computation 273 (2016) 1051–1058

(a)

1055

(b)

Fig. 1. Example 1 - ﬁnding roots of f1 (t ) = (t − 1/3) · (4 − t )6 , t ∈ [0, 1]: (a) the plots of f1 (t) and its bounding curves from different methods; and (b) the error plots between f1 (t) and its bounding curves in different methods.

| f (t ) − R j (t )| ≤ | f (t ) − Y j (t )| + |Y j (t ) − R j (t )| ≤ (M1 + M2 + j )(b − a)5 , j = l, r, d

where j = max{ w4, j , 1, j

d4, j w2, j

(18)

}. From Eq. (18), the bounding rational cubics Rl (t) and Rr (t) achieve the approximation order 5 to f(t).

Combining with Theorem 1, we have the following theorem. Theorem 3. The new method by using Rl (t) and Rr (t) can achieve the convergence rate 5/k for a root t , where k is the multiplicity of the root t . 3.4. Illustrations This section shows an example to illustrate the rational clipping method of linear computation complexity. For the sake of convenience, let M3 , M4 , Mr and MI be the cubic clipping method in [13], the quartic clipping method by using Eq. (11), the rational cubic clipping method by using Eq. (14) and the improved method mentioned in Remark 2, respectively. Example 1. Given f1 (t ) = (t − 1/3) · (4 − t )6 , t ∈ [0, 1]. Firstly, we compute the bounding quartics by using M4 , and then compute the bounding rational cubics by using Mr . As shown in Fig. 1(a), the given curve f(t) is in solid black, while the curves in dotted green, in solid blue and in dashed red are the bounding curves from M3 , M4 and Mr , respectively. The error plots between f(t) and the bounding curves are shown in Fig. 1(b). The maximum error between f(t) and the bounding curves of M3 , M4 and Mr are 22.044, 3.605 and 12.054, respectively, it shows that the errors of both M4 and Mr are much less than that of M3 . The length of the resulting subintervals from M3 and Mr are 1.66e-2 and 2.15e-3, and it shows that Mr achieves much tighter bounds than (5) that of M3 . In this case, f1 (t ) ≥ 0, ∀t ∈ [0, 1], by using MI mentioned in Remark 2, we can obtain an improved result by directly bounding f1 (t), the corresponding maximum error and the length of the resulting subinterval are 2.613 and 4.08e-4, which are much better than those of M3 , M4 and Mr . More comparison results between M3 and Mr are listed in Table 1, it shows that M3 and Mr achieve the convergence rates 4 and 5, respectively, and Mr has much better approximation effect and much higher computation eﬃciency than those of M3 . As shown in Table 1, when the length b − a of the given interval [a, b] tends to be zero, the results of both Mr and MI tends to be the same. 4. Further numerical examples and related discussions 4.1. Comparison on computational complexity In a clipping step of previous cubic clipping method M3 , it ﬁrstly ﬁnds a cubic polynomial g(t) to approximate f(t), and then turns the error polynomial f (t ) − g(t ) into a Bézier form within [a, b] for estimating the error bounds ε1 ≤ f (t ) − g(t ) ≤ ε2 within

1056

X.-d. Chen, W. Ma / Applied Mathematics and Computation 273 (2016) 1051–1058

(a)

(b)

(c)

Fig. 2. Plots of Examples 2–4 (a) f2 (t); (b) f3 (t); and (c) f4 (t). Table 1 Comparison of errors and convergence rates for Example 1. Method

l1

l2

l3

l4

l5

Time(ms)

CR

M3 Mr MI

1.6e-2 2.1e-3 4.0e-4

1.3e-9 1.0e-16 2.5e-17

6.8e-38 2.3e-84 5.4e-85

3.9e-151 5.8e-424 1.3e-424

8.7e-600 / /

30.74 14.54 6.87

4 5 5

M3 : The cubic-clipping method in [13]. Mr : The rational cubic clipping method in this paper. MI : The improved method mentioned in Remark 2. Time: The average computation time of a clipping step (in millisecond). CR: Convergence rate. Table 2 Comparisons on computational complexity. Method

M3 2

Time

O(n )

Mr

MI

O(n)

O(n)

Table 3 Average computation time of a clipping step. Example

Fig. 1

Fig. 2(a)

Fig. 2(b)

Fig. 2(c)

n T3 Tr TI

7 30.74 14.54 6.87

16 45.40 16.16 14.76

22 83.15 28.24 26.84

13 42.43 20.12 16.22

O(n2 ) time. Thus, g(t ) + ε1 and g(t ) + ε2 bound f(t). Finally, one solves the roots of the bounding cubics in O(1) time, and utilizes these roots for clipping the subintervals which contain no roots of f(t). The total computational complexity of M3 is O(n2 ). By verifying the signs of the coeﬃcients of the items in power series, one can easily obtain both f+ (t ) and f− (t ) in O(n) time. (5) (5) In each clipping step, the rational cubic method has the following steps. Firstly, note that f+ (t ) ≥ 0 and − f− (t ) ≥ 0, the values of f(ti ) and f (tκ ) in Eq. (7) can be computed in O(n) time, by using Formula (7) whose solution can be explicitly represented in O(1) time, one can compute the bounding quartics Yi, 4 (t) in O(n), i = 1, 2, 3, 4. Secondly, by using Formula (11), one obtains quartics Yl (t) and Yr (t) for bounding f(t) in O(1) time. Thirdly, by using the method in Section 3.3, one can compute the rational cubics R1, j (t) and R2, j (t) for bounding the quartic Yj in O(1) time, j = l, r. Thus, we have that min {R1, l (t), R2, l (t)} ≤ Yl (t) ≤ f(t) ≤ Yr (t) ≤ max {R1, r (t), R2, r (t)}. Finally, similarly as that of M3 , one ﬁnds the roots of the bounding rational cubics in O(1) time, and utilizes these roots for clipping the subintervals which contain no roots of f(t). The total computational complexity is O(n). By using the improved method MI , which directly bounds f(t) without dividing it into f+ (t ) and f− (t ), one can also do a clipping step in O(n) time. The comparison results are listed in Table 2. We have tested the four examples in Figs. 1 and 2, the average computation time of a clipping step of the methods M3 , Mr and MI , i.e., T3 , Tr and TI , respectively, are listed in Table 3. It shows that Mr and MI achieve a much better computation eﬃciency than that of M3 . For a large n, the computational time of Mr and

X.-d. Chen, W. Ma / Applied Mathematics and Computation 273 (2016) 1051–1058

1057

Table 4 Comparison on errors and convergence rates for Examples 2–4. Exam

Method

l1

l2

l3

l4

l5

Time(ms)

CR

f2 (t)

M3 Mr MI M3 Mr MI M3 Mr MI

1.6e-2 7.0e-2 7.0e-2 1.0e-1 1.0e-1 1.0e-1 2.1e-1 2.1e-1 2.1e-1

1.2e-9 3.7e-8 2.8e-8 1.4e-3 1.3e-3 1.3e-3 4.9e-2 4.2e-2 4.2e-2

3.4e-38 1.5e-39 2.1e-40 6.5e-13 1.9e-13 8.0e-16 2.5e-3 2.0e-4 2.0e-4

2.2e-152 1.9e-196 5.1e-201 5.1e-51 1.4e-60 6.8e-77 5.4e-6 2.8e-11 2.7e-11

5.5e-607 3.5e-981 / 1.8e-203 1.3e-301 3.0e-382 2.4e-11 9.9e-32 8.6e-32

45.40 16.16 14.76 83.15 28.24 26.84 42.43 20.12 16.22

4 5 5 4 5 5 2 2.5 2.5

f3 (t)

f4 (t)

that of MI are very close to each other, but are much better than that of M3 . For cases that 7 ≤ n ≤ 22, Tr and TI are about r ∈ [1/3, 1/2] times of T3 . 4.2. Further numerical examples and discussions In this section, we show several examples for comparing the rational clipping method Mr with the cubic clipping method M3 in [13]. All of the examples are tested by using the Maple software on PC with 2.5G CPU and 4G Memory, the number of the digits after point is 16, and the unit is millisecond. Note that there may be three or more roots within the given interval [a, b], and the corresponding convergence rate is near 1 or less, which converges very slowly. To avoid this, in this section, at the beginning, if there are three or more zeros of the control polygon, we use these zeros of the control polygon to divide the given interval into several subintervals, which can isolate most of the real roots within [a, b]. Examples 2–4. We have tested the following three examples within [0, 1], where f2 (t) and f3 (t) have one and six single roots respectively, while f4 (t) has a double root and two single roots.

f2 (t ) = (t − 1/3) · (2 − t )5 · (t + 5)10 ; f3 (t ) = (t − 0.10001) · (t − 0.1000002) · (t − 0.234) · (t − 0.45) · (t − 0.5) · (t − 0.8) · (t + 0.1)2 · (t − 2)14 ; f4 (t ) = (t − 1/3)2 · (t − 1/4) · (t − 3/4) · (t + 5)7 · (t − 6)2 . We have tested the single root t = 1/3 of f2 (t), the single root t = 0.234 of f3 (t) and the double root t = 1/3 of f4 (t), and the corresponding results of M3 , Mr and MI are listed in Table 4. It shows that M3 achieves the convergence rate 4/k, while both Mr and MI achieve the convergence rate 5/k. The computational eﬃciencies of Mr and MI are much better than that of M3 . 4.3. Numerical robustness The numerical robustness of Mr and MI is similar to that of the cubic clipping method M3 in [13]. The Mr method also uses the Cardano formula to solve cubic polynomial equations, which is the same to M3 . The stability of BB representation is also applicable to Mr , which has the same robustness to that of M3 . Similarly as done in [13], we have tested both M3 and Mr to compute the roots of the Wilkinson polynomial

W (x) =

20

(x − i),

i=0

and both M3 and Mr achieve a comparable or almost the same robustness. 5. Conclusions This paper presents a rational cubic clipping method for ﬁnding the roots of a polynomial f(t). It can achieve a convergence rate of 5 for a single root by using rational cubics as bounding polynomials. Different from previous clipping methods in [1,13] for computing two bounding polynomials in O(n2 ) time, the new rational cubic clipping method Mr directly constructs two rational cubic polynomials, which can be used to bound f(t) for cases that the length of [a, b] is small enough and leads to an O(n) computational complexity. Numerical examples also show that the new rational cubic clipping method can achieve a higher convergence rate, much better approximation effect and much higher computation eﬃciency than that of previous clipping methods in [1,13]. Acknowledgment This research was partially supported by City University of Hong Kong (SRG 7004245) and the National Science Foundation of China (61370218, 61370166).

1058

X.-d. Chen, W. Ma / Applied Mathematics and Computation 273 (2016) 1051–1058

References ˇ B. Jüttler, Computing roots of polynomials by quadratic clipping, Comput. Aided Geom. Des. 24 (2007) 125–141. [1] M. Barton, [2] X.D. Chen, J.H. Yong, G.Z. Wang, J.C. Paul, G. Xu, Computing the minimum distance between a point and a NURBS curve, Comput.-Aided Des. 40 (10-11) (2008) 1051–1054. [3] X.D. Chen, W.Y. Ma, A planar quadratic clipping method for computing a root of a polynomial in an interval, Comput. Gra. 46 (2) (2015) 89–98. [4] X.D. Chen, W.Y. Ma, Y.T. Ye, A rational quadratic clipping method for computing roots of a polynomial, J. Comput. Appl. Math. (2015). (submiting for publication). [5] Y.K. Choi, W.P. Wang, Y. Liu, M.S. Kim, Continuous collision detection for two moving elliptic disks, IEEE Trans. Robot. 22 (2) (2006) 213–224. [6] P.J. Davis, Interpolation and approximation, Dover Publications, New York, 1975. [7] A. Efremov, V. Havran, H.-P. Seidel, Robust and numerically stable Bzier clipping method for ray tracing NURBS surfaces, in: Proceedings of the 21st Spring Conference on Computer Graphics, New York, ACM Press, pp. 127–135. [8] G. Elber, M.-S. Kim, Geometric constraint solver using multivariate rational spline functions, In: Proceedings of the Sixth ACM/IEEE Symposium on Solid Modeling and Applications, Ann Arbor, Michigan, pp. 1–10. [9] R.T. Farouki, T.N.T. Goodman, On the optimal stability of the Bernstein basis, Math. Comput. 126 (65) (1996) 1553–1566. [10] E. Isaacson, H.B. Keller, Analysis of Numerical Methods, Wiley, New York, 1966. [11] B. Jüttler, The dual basis functions of the Bernstein polynomials, Adv. Comput. Math. 8 (1998) 345–352. [12] M. Lin, S. Gottschalk, Collision detection between geometric models: a survey, in: Proceedings of IMA Conference on Mathematics of Surfaces, Birmingham, UK, 1998, pp. 37–56. [13] L.G. Liu, L. Zhang, B.B. Lin, G.J. Wang, Fast approach for computing roots of polynomials using cubic clipping, Comput. Aided Geom. Des. 26 (2009) 547–559. [14] J.M. McNamee, 1993–2002. Bibliographies on roots of polynomials. J. Comp. Appl. Math. 47, 391–394; 110, 305–306; 142, 433–434, available at http:// www1.elsevier.com/homepage/sac/cam/mcnamee/. [15] J.-P. Merlet, Parallel robots. Solid Mechanics and Its Applications, Kluwer, 2005. [16] K. Morken, M. Reimers, An unconditionally convergent method for computing zeros of splines and polynomials, Math. Comput. 76 (258) (2007) 845–865. [17] B. Mourrain, J.-P. Pavone, Subdivision methods for solving polynomial equations, Technical report, no. 5658, INRIA Sophia Antipolis, 2005. http://www.inria. fr/rrrt/rr-5658.html. [18] T. Nishita, T.W. Sederberg, M. Kakimoto, Ray tracing trimmed rational surface patches, In: Proceedings of Siggraph. ACM, pp. 337–345. [19] N.M. Patrikalakis, T. Maekawa, Intersection problems, in: FarinG. HoschekJ. KimM.-S. (Eds.), Handbook of Computer Aided Geometric Design. Elsevier, pp. 623–649. [20] M. Reuter, T. Mikkelsen, E. Sherbrooke, T. Maekawa, N. Patrikalakis, Solving nonlinear polynomial systems in the barycentric Bernstein basis, Vis. Comput. 24 (3) (2007) 187–200. [21] T.W. Sederberg, T. Nishita, Curve intersection using Bzier clipping, Comput.-Aided Des. 22 (9) (1990) 538–549.