- Email: [email protected]
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
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 efficiency of the new method. The method is particularly useful for the fly computation in many geometry processing and graphics rendering applications. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Optimal approximation order Root finding 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 fly 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 finding 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 firstly 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 inflexion 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 finding method. Section 3 presents the new root finding 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 first 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 figure 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 find 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 inflexion 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 verified 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 simplified 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 figure 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 find 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 find ϕ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 satisfied 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 pffiffiffiffiffiffiffiffiffiffiffiffithe 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 verified 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 find 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 simplified 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 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 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 simplified 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 figure 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 pffiffiffiffiffiffiffiffiffiffiffiffi 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 verified that ϕ4 ð0Þ ¼ 0, ϕ4 ð1Þ ¼ 1, p0 ¼ Y 4 ð0Þ and p1 ¼ Y 4 ð1Þ. In many cases, we can find two unknowns ρ4;0 and ρ4;2 such that 0 d1 ¼ Y 4 ð1Þ:
ð21Þ
0 4 ðtÞ Z 0
is not satisfied 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 simplified 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 first 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 verified 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 firstly 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 pffiffiffiffiffiffi pffiffiffiffiffiffi 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 infinitively 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 finding 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 first 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 satisfied 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 firstly 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 verified 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 satisfied. 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 verified 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 figure 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 figure 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 find a reference polynomial; (2) to compute two bounding polynomials and their roots within the given interval; and (3) to find 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 verified 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 satisfied, 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 satisfied, 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 first clipping process is tested by using 17 valid digits after the point, which is similar to the float 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 finding 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 efficiency. 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. Efficient 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 fitting – 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.
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 efficiency of the new method. The method is particularly useful for the fly computation in many geometry processing and graphics rendering applications. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Optimal approximation order Root finding 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 fly 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 finding 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 firstly 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 inflexion 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 finding method. Section 3 presents the new root finding 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 first 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 figure 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 find 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 inflexion 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 verified 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 simplified 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 figure 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 find 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 find ϕ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 satisfied 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 pffiffiffiffiffiffiffiffiffiffiffiffithe 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 verified 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 find 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 simplified 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 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 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 simplified 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 figure 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 pffiffiffiffiffiffiffiffiffiffiffiffi 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 verified that ϕ4 ð0Þ ¼ 0, ϕ4 ð1Þ ¼ 1, p0 ¼ Y 4 ð0Þ and p1 ¼ Y 4 ð1Þ. In many cases, we can find two unknowns ρ4;0 and ρ4;2 such that 0 d1 ¼ Y 4 ð1Þ:
ð21Þ
0 4 ðtÞ Z 0
is not satisfied 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 simplified 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 first 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 verified 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 firstly 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 pffiffiffiffiffiffi pffiffiffiffiffiffi 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 infinitively 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 finding 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 first 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 satisfied 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 firstly 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 verified 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 satisfied. 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 verified 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 figure 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 figure 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 find a reference polynomial; (2) to compute two bounding polynomials and their roots within the given interval; and (3) to find 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 verified 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 satisfied, 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 satisfied, 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 first clipping process is tested by using 17 valid digits after the point, which is similar to the float 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 finding 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 efficiency. 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. Efficient 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 fitting – 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.