nearly optimal
678
IEEE TRANSACTIONS ON COMPUTERS,
VOL. 53,
NO. 6,
JUNE 2004
Near Optimality of Chebyshev Interpolation for Elementary Function Computations Ren-Cang Li Abstract—A common practice for computing an elementary transcendental function in an libm implementation nowadays has two phases: reductions of input arguments to fall into a tiny interval and polynomial approximations for the function within the interval. Typically, the interval is made tiny enough so that polynomials of very high degree aren’t required for accurate approximations. Often, approximating polynomials as such are taken to be the best polynomials or any others such as the Chebyshev interpolating polynomials. The best polynomial of degree n has the property that the biggest difference between it and the function is smallest among all possible polynomials of degrees no higher than n. Thus, it is natural to choose the best polynomials over others. In this paper, it is proven that the best polynomial can only be more accurate by at most a fractional bit than the Chebyshev interpolating polynomial of the same degree in computing elementary functions or, in other words, the Chebyshev interpolating polynomials will do just as well as the best polynomials. Similar results were obtained in 1967 by Powell who, however, did not target elementary function computations in particular and placed no assumption on the function and, remarkably, whose results imply accuracy differences of no more than 2 to 3 bits in the context of this paper. Index Terms—Elementary function computation, libm, Chebyshev Interpolation, Remez, best polynomial, accuracy.
æ 1
INTRODUCTION
A
common practice for computing elementary transcendental functions fðxÞ has two phases (see [1], [2], [3], [4], [5], [7], [8], [9], [10], and references therein): An argument reduction that typically reduces x to fall into a very narrow interval ½; , e.g., 25 or smaller. 2. An efficient polynomial approximation1 to fðxÞ for x 2 ½; . Often, approximating polynomials are forced to match fðxÞ at x ¼ 0 for computational advantages. This is especially true when fð0Þ ¼ 0. Small ensures the use of approximating polynomials of low degrees and, consequently, little time to evaluate them. Often, is very tiny, especially for table-driven algorithms [3], [4], [5], [7], [8], [9], [10]. There are many ways to find approximating polynomials. Given the degree n of the polynomials, common ones are 1) finding the best polynomial pn ðxÞ that minimizes maxx2½; jfðxÞ pn ðxÞj (or a variation of this if fð0Þ is to be preserved; see below for details) and 2) finding the Chebyshev interpolating polynomial. Rational approximation is another alternative, but it is seldom used in elementary function computations because division is much slower than addition/multiplication and, moreover, the division can only be initiated after the availability of the denominator and/or the numerator. 1.
1. There may be exceptions to this, e.g., computing sin x may require polynomial approximations to both sin x and cos x, depending on algorithms.
. The author is with the Department of Mathematics, University of Kentucky, Lexington, KY 40506. E-mail: [email protected]. Manuscript received 3 Apr. 2003; revised 6 Nov. 2003; accepted 17 Nov. 2003. For information on obtaining reprints of this article, please send e-mail to: [email protected], and reference IEEECS Log Number 118532. 0018-9340/04/$20.00 ß 2004 IEEE
This paper compares the approximation accuracy achieved by the two polynomial approximation techniques. Earlier results along this line are due to Powell [6]. However, Powell’s paper went largely unnoticed in the literature of elementary function computations. Even so, through practice it has become folklore [3], [4] that the Chebyshev interpolating polynomial will be just as effective as the best polynomial with the same degree. In fact, most error analysis for the polynomial solution part in Markstein [3] is based on the Chebyshev polynomial interpolation. In this paper, we will show that the Chebyshev interpolating polynomial may be less accurate, if indeed, than the best polynomial of the same degree by only a fractional bit. The best polynomials are typically obtained through a slowly convergent iterative procedure—the Remez procedure [11], while we do have a close formula for the Chebyshev interpolating polynomials. This does not, however, necessarily place the best polynomial approach at a disadvantage since all approximating polynomials are precomputed in any libm implementations [1], [12] and, thus, the cost in computing approximating polynomials becomes somewhat irrelevant. Even so, we feel that our contribution here has its significance, that is, that simple Chebyshev interpolation is sufficient and that the common practice of feeding Chebyshev interpolation polynomials into the Remez procedure offers negligible improvement, if any, and thus may be considered unnecessary. Our results compare favorably to Powell’s results which, applied to the context here, say the difference in accuracy is no more than 2 to 3 bits. It is worth mentioning that Powell’s results were obtained, remarkably, without any assumption on fðxÞ. The rest of this paper is organized as follows: Section 2 briefly introduces the general interpolation theory, enough to serve our purpose later on. Chebyshev interpolation and Published by the IEEE Computer Society
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
best polynomial approximations are discussed in Sections 3 and 4 with emphasis on error expansions on these approximations. Our main theorems that compare the accuracies by the two approaches are given in Section 5, where a theorem of Powell’s will also be discussed. An asymptotic analysis in Section 6 on the main theorems confirms our claim made at the beginning of this paper and concludes that the approximation accuracies between the best polynomials and the Chebyshev interpolating polynomials differ by Oð2 Þ bits for odd and even functions and by OðÞ for others. Finally, in Section 8, the main theorems are applied to the most common elementary transcendental functions, such as expðxÞ, ln x, sin x, and so on.
1.1 Notation Throughout the remainder of this paper, the interval of interest is ½; on which fðxÞ is defined. It is assumed that fðxÞ is continuous and has continuous derivatives up to ðn þ 1Þth (or more as deemed necessary), where n is the degree of polynomials sought to approximate fðxÞ. If needed, 0 2 ð; Þ is also implied. Define f0 ðxÞ :¼
fðxÞ fð0Þ ; xm
ð1Þ
where m is the multiplicity of x ¼ 0 as a zero of fðxÞ fð0Þ, and “:=” is the assignment operator. Usually, m ¼ 1, but may not be, e.g., for cos x it is 2. Unknown points ; 2 ½; in error equations later may be different from one occurrence to another.
2
GENERAL POLYNOMIAL INTERPOLATION
At times, polynomials that preserve certain properties of fðxÞ are called for so that simple but important mathematical identities, such as log 1 ¼ 0, expð0Þ ¼ 1, sin 0 ¼ 0, and cosð0Þ ¼ 1, are numerically conserved. In elementary function computations, this typically happens when 0 2 ½; (and ¼ as well, but we shall not assume this for now). To make sure that identities as such remain valid in the face of roundoffs, we need to place2 x ¼ 0 to the list of interpolating points to ensure fð0Þ ¼ pn ð0Þ exactly! One way to implement this is to interpolate f0 ðxÞ at n m þ 1 other points, x~0 ; x~1 ; . . . ; x~nm ;
ð2Þ
be interpolating points in the interval. The interpolating polynomial pn ðxÞ is uniquely determined by pn ðxÞ and fðxÞ coinciding at those points in the sense that, if xi is repeated k times,
p~n ðxÞ :¼ fð0Þ þ xm hnm ðxÞ; a polynomial of degree n that can be seen to coincide with fðxÞ at x ¼ 0, as well as at x~i s. An application of (3) gives fðxÞ p~n ðxÞ ¼
ð5Þ
for some minfx; 0; x~0 ; . . . ; x~nm g maxfx; 0; x~0 ; . . . ; x~nm g:
3
CHEBYSHEV INTERPOLATION
The standard Chebyshev polynomial interpolation for a function uðtÞ is on the closed interval ½1; 1. The idea is to seek a polynomial qn ðtÞ of degree n that matches exactly uðtÞ at the zeros ti :¼ cos
2i þ 1 ; nþ1 2
i ¼ 0; 1; . . . ; n;
ð6Þ
of Tnþ1 ðtÞ cosðn þ 1Þ arccos t, the ðn þ 1Þth Chebyshev polynomial, i.e., uðti Þ ¼ qn ðti Þ for i ¼ 0; 1; . . . ; n. It can be proven that (see, e.g., [3, p. 99]) qn ðtÞ ¼
n X
k Tk ðtÞ;
ð7Þ
where
(So, if k ¼ 1, only the function value matches at xi .) If none of the points xi collapse, this is called Lagrange Interpolation. In the extreme case when all xi s collapse to one point, pn ðxÞ is just the truncated Taylor series of fðxÞ at the collapsed point. In any case, the difference of fðxÞ and pn ðxÞ is given as follows: If f ðnþ1Þ ðxÞ exists and is continuous, then n f ðnþ1Þ ðÞ Y ðx xi Þ; ðn þ 1Þ! i¼0
Y f ðnþ1Þ ðÞ m nm x ðx x~i Þ ðn þ 1Þ! i¼0
k¼0
ðjÞ pðjÞ n ðxi Þ ¼ f ðxi Þ for j ¼ 0; . . . k 1:
fðxÞ pn ðxÞ ¼
ð4Þ
by a polynomial hnm ðxÞ of degree n m. Then, set
Polynomial interpolation is an ancient idea that is still widely in use today to approximate a function. Let x0 ; x1 ; . . . ; xn
679
ð3Þ
for some minfx; x0 ; . . . ; xn g maxfx; x0 ; . . . ; xn g [13, p. 60]. This equation offers one way to find a good polynomial to make max jfðxÞ pn ðxÞj small. Since f ðnþ1Þ ðÞ is function-dependent and, often, Q is unknown, it is natural to look for xi s so that max j ni¼0 ðx xi Þj is minimized. It turns out this minimization problem can be solved elegantly: Pick xi s as the translated zeros of the ðn þ 1Þth Chebyshev polynomial (see Theorem 3.2 below).
k ¼
Pn 1 nþ1 Pi¼0 n 2 i¼0 nþ1
uðti Þ; if k ¼ 0; Tk ðti Þuðti Þ; otherwise:
ð8Þ
Theorem 3.1. If u is odd (even), then qn is also odd (even). Proof. This is a well-known property of Chebyshev interpolation. But, we shall give a proof here for completeness. Notice that Tnþ1 ðtÞ ¼ cosðn þ 1Þ arccosðtÞ ¼ cos½ðn þ 1Þð arccos tÞ ¼ ð1Þnþ1 Tnþ1 ðtÞ: 2. Mathematically, it does not have to be 0 and, in fact, we can preserve the function value at any given floating-point number so long as the function value is also a floating-point number. Such a generalization, however, is no more general than 0 itself because shifting of origin can always make it 0!
680
IEEE TRANSACTIONS ON COMPUTERS,
Tnþ1 is an even function if n þ 1 is even and an odd function otherwise. Therefore, ti and ti appear in pairs in the list (6). Suppose u is odd, i.e., uðtÞ ¼ uðtÞ. Let sn ðtÞ :¼ qn ðtÞ, also a polynomial of degree n. We have sn ðti Þ ¼ qn ðti Þ ¼ uðti Þ ¼ uðti Þ; which says that sn also interpolates u at the same points (6) as qn does. Since the interpolating polynomial is unique, we conclude that qn ðtÞ sn ðtÞ ¼ qn ðtÞ, i.e., qn is odd. Now, suppose u is even, i.e., uðtÞ ¼ uðtÞ. Let sn ðtÞ :¼ qn ðtÞ, a polynomial of degree n. We have sn ðti Þ ¼ qn ðti Þ ¼ uðti Þ ¼ uðti Þ; which says that sn also interpolates u at the same points (6) as qn does. The uniqueness of the interpolating polynomial leads to qn ðtÞ sn ðtÞ ¼ qn ðtÞ, i.e., qn is u t
even.
Now, we turn to the continuous function fðxÞ defined on an arbitrary closed interval ½; . The idea is first to devise a transformation x ¼ xðtÞ such that gðtÞ :¼ fðxðtÞÞ is defined on ½1; 1. One such commonly used transformation is x
: ½1; 1 ! t !
½; þ xðtÞ :¼ 2 tþ 2 :
The inverse transformation of xðtÞ is 2 þ x : tðxÞ :¼ 2
NO. 6,
JUNE 2004
n n Y nþ1 Y 2 þ x ti ðx xi Þ ¼ 2 2 i¼0 i¼0 nþ1 1 Tnþ1 ðtÞ; ¼ n 2 2 ð13Þ where the last equality holds because ti s are all the zeros of Tnþ1 whose leading coefficient is 2n , t is as defined by (10). Similarly, n Y nþ1 ðx yi Þ ¼ qnþ1 ðtÞ; 2 i¼0 where qnþ1 ðtÞ is polynomial in t of degree n with leading coefficient 1. It is known [13, p. 61] that, on ½1; 1, Tnþ1 ðtÞ=2n takes the least maximum absolute value among all polynomials of degree n þ 1 with leading coefficient 1. Therefore, (11) holds. u t It follows from (7) and (8) that the Chebyshev interpolating polynomial cn ðxÞ of degree n for fðxÞ on ½; is given by n X 2 þ x ; ð14Þ k Tk cn ðxÞ ¼ 2 k¼0
ð9Þ
VOL. 53,
k ¼
Pn 1 nþ1 Pi¼0 n 2 i¼0 nþ1
fðxi Þ; Tk ðti Þfðxi Þ;
if k ¼ 0; otherwise:
ð15Þ
The proof of Theorem 3.2 also yields ð10Þ
Write, for i ¼ 0; 1; . . . ; n,
fðxÞ cn ðxÞ ¼
n f ðnþ1Þ ðÞ Y ðx xi Þ ðn þ 1Þ! i¼0
1 nþ1 f ðnþ1Þ ðÞ Tnþ1 ðtÞ; ¼ n 2 2 ðn þ 1Þ!
þ ti þ ; xi :¼ xðti Þ :¼ 2 2
ð16Þ
ð17Þ
called the translated zeros in ½; of Tnþ1 . The polynomial
for some 2 ½; and t is defined as in (10).
cn ðxÞ of degree n such that fðxi Þ ¼ cn ðxi Þ for all i is called
Theorem 3.3. Suppose ¼ . Then, if f is odd (even), cn is odd (even).
the Chebyshev interpolating polynomial of degree n for fðxÞ. The general theory described in Section 2 still applies here, so the error equation (7) holds with xi assigned here and for some 2 ½; . Theorem 3.2. The translated zeros xi s of Tnþ1 satisfy Y Y n n max ðx xi Þ max ðx yi Þ x x i¼0 i¼0
ð11Þ
for any y0 ; y1 ; . . . ; yn 2 ½; . Proof. This is also well-known. Again, a proof is given for completeness; also, a key equation (13) is useful later. We have þ ti x xi ¼ x 2 2 2 þ x ti ; ¼ 2 2
ð12Þ
Proof. ¼ implies xðtÞ ¼ t. If fðxÞ is odd (even) in x, fðtÞ is odd (even) in t. Hence, by Theorem 3.1, cn ðtÞ is u t odd (even) in t and, thus, cn ðxÞ is odd (even) in x. Tnþ1 ðtÞ has a zero at t ¼ 0 for even n. Thus, if ¼ and n is even, the Chebyshev interpolating polynomial cn ðxÞ automatically satisfies cn ð0Þ ¼ fð0Þ. But, it is a different story when n is odd. As we pointed out in Section 2, we can ensure that fð0Þ is always recovered exactly by looking for a Chebyshev interpolating polynomial hnm ðxÞ of degree n m for f0 ðxÞ and then setting c~n ðxÞ :¼ fð0Þ þ xm hnm ðxÞ: Notice that, now, hnm ðxÞ coincides with f0 ðxÞ at the translated zeros x~i s of Tnmþ1 . It is interesting to notice that Tnmþ1 ð0Þ ¼ 0 also for odd n m þ 1; this implies x ¼ 0 appears m þ 1 times in the interpolating points for fðxÞ when ¼ ; so, for odd n m þ 1 and ¼ , we have
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
fð0Þ ¼ c~n ð0Þ; f 0 ð0Þ ¼ c~0n ð0Þ; . . . ; f ðmÞ ð0Þ ¼ c~ðmÞ n ð0Þ: With the arguments similar to those which lead to the error equation (17), we can prove fðxÞ c~n ðxÞ fðxÞ ½fð0Þ þ xm hnm ðxÞ Y f ðnþ1Þ ðÞ m nm x ¼ ðx x~i Þ ðn þ 1Þ! i¼0
¼
1 2nm
nmþ1 2
f ðnþ1Þ ðÞ m x Tnmþ1 ðtÞ; ðn þ 1Þ!
ð18Þ
fðxÞ cn ðxÞ ¼
2nþ1
f ðnþ2Þ ðÞ Tnþ2 ðtÞ; ðn þ 2Þ!
1 nmþ2 2nmþ1 f ðnþ2Þ ðÞ m x Tnmþ2 ðtÞ; ðn þ 2Þ!
ð20Þ
fðxÞ c~n ðxÞ ¼
ð21Þ
where ; 2 ½; and t is defined as in (6).
4
BEST APPROXIMATION
Given a continuous function fðxÞ defined on a closed interval ½; , the best polynomial approximation bn ðxÞ of degree n for fðxÞ is the one such that n :¼ max jfðxÞ bn ðxÞj :¼ min max jfðxÞ pn ðxÞj; x
such that fðzi Þ ¼ bn ðzi Þ for i ¼ 0; 1; . . . ; n. So, the best polynomial may well be thought of as a Lagrange interpolation polynomial with interpolating points not known a priori and, therefore,
ð19Þ
Theorem 3.4. Assume ¼ and f is odd (even). If n is odd (even), then nþ2
By the continuity of the functions, we conclude that, between every two consecutive yi s, fðxÞ bn ðxÞ has a zero, i.e., there are y0 < z0 < . . . < yi < zi < yiþ1 < . . . < zn < ynþ1
for some 2 ½; and t is defined as in (10). It can be seen that Theorem 3.3 holds with cn replaced by c~n . There is an important implication to this fact and Theorem 3.3. Assume ¼ and n is odd (or even), depending on whether f is odd (or even). For better error bounds, we shall always think of interpolating f (or f0 ) at the translated zeros of Tnþ2 (or of Tnmþ2 ), in contrast to Tnþ1 (or Tnmþ1 ) for the general f. In general, this will produce a Chebyshev interpolating polynomial of degree n þ 1, but, since it has to be odd (or even), its leading coefficient must be zero! and thus degenerates to a polynomial of degree n.
1
681
fðxÞ bn ðxÞ ¼
n f ðnþ1Þ ðÞ Y ðx zi Þ; ðn þ 1Þ! i¼0
ð22Þ
for some 2 ½; . Such reinterpretation is usually not useful from a numerical point of view, but it is a key step in deriving our main theorems in the next section. For the same reason as explained in Sections 2 and 3, it may also be advantageous to work with f0 ðxÞ for the sake of simple but important mathematical identities. Let hnm ðxÞ be the best polynomial for f0 ðxÞ in the sense that max jf0 ðxÞ hnm ðxÞj ¼ min max jf0 ðxÞ pnm ðxÞj pnm x
x
taken over all possible polynomials pnm ðxÞ of degree n m, and set b~n ðxÞ :¼ fð0Þ þ xm hnm ðxÞ; which approximates fðxÞ over the interval. Of course, then b~n ðxÞ may do worse than3 bn ðxÞ, but the trade off is that now fð0Þ ¼ b~n ð0Þ always! Similarly, we have the following error equation Y f ðnþ1Þ ð~Þ m nm fðxÞ b~n ðxÞ ¼ ðx z~i Þ; x ðn þ 1Þ! i¼0
ð23Þ
for some ~ 2 ½; and z~0 < z~1 < . . . < z~nm . Even though, now, b~n ðxÞ is no longer the best polynomial for fðxÞ, we still call it the best polynomial at places where no confusion may rise. It can be seen that the uniqueness of the best polynomial implies that bn and b~n are odd (even) if f is odd (even) and ¼ . This implies: Theorem 4.2. Assume ¼ and f is odd (even). If n is odd (even), then
pn x
where the min is taken over all possible polynomials pn ðxÞ of degree n. It is known [13], [14, p. 58], [15] that bn ðxÞ exists and is unique and is characterized by the following theorem, traditionally called the Chebyshev Theorem, although he did not prove it (see [14, p. 58] for a detailed history of it). Theorem 4.1 (Chebyshev). Assume f is not a polynomial of degree less than or equal to n. Then, there are at least n þ 2 points y0 < y1 < . . . < ynþ1 at which fðyi Þ bn ðyi Þ ¼ sð1Þi n for 0 i n þ 1, where s ¼ 1, independent of i. In other words, fðxÞ pn ðxÞ travels from one extreme point to another of the opposite sign for at least n þ 1 times.
bn ¼ bnþ1 ;
b~n ¼ b~nþ1 :
Consequently, fðxÞ bn ðxÞ ¼
þ1 f ðnþ2Þ ðÞ nY ðx zi Þ; ðn þ 2Þ! i¼0
ð24Þ
fðxÞ b~n ðxÞ ¼
Y f ðnþ2Þ ð~Þ m nmþ1 ðx z~i Þ; x ðn þ 2Þ! i¼0
ð25Þ
where zi ; ; z~i ; ~ 2 ½; . Proof. Consider the case where n is odd and f is odd only. The other case can be proven in the same way. All the 3. It may be true the other way. See Remark 5.1.
682
IEEE TRANSACTIONS ON COMPUTERS,
best polynomials of f are odd, so bnþ1 is odd. But, since n þ 1 is even, a necessary condition for bnþ1 to be odd is its leading coefficient is zero, i.e., bnþ1 ’s degree is no higher than n. Now, the uniqueness of the best polynomials implies bn ¼ bnþ1 . b~n ¼ b~nþ1 can be proven in the same way. u t
5
VOL. 53,
NO. 6,
JUNE 2004
5.1 Approximations by bn ðxÞ and cn ðxÞ Theorem 5.1. Let bn ðxÞ and cn ðxÞ be defined by (25) and (26). We have 1
cheb ðf ðnþ1Þ Þ; best
ð29Þ
1 cheb ðfÞ ðf ðnþ1Þ Þ: ðfÞ best
MAIN THEOREMS
We shall use BestApproxn ðfÞ :best approximation polynomial; ChebApproxn ðfÞ :Chebyshev interpolating polynomial; of degree n, as building blocks to construct polynomials. Specifically, we will study two sets of approximation polynomials of degree n defined either as bn ðxÞ :¼ BestApproxn ðfÞðxÞ;
ð25Þ
cn ðxÞ :¼ ChebApproxn ðfÞðxÞ;
ð26Þ
ð30Þ
Proof. best cheb holds since bn ðxÞ is the best among all polynomials of degree n. This proves the leftmost inequality in (29). We now turn to the rightmost inequality in (29). Our proof relies on (3) with translated zeros xi s of Tnþ1 , (22), and Theorem 3.2. We have Y n 1 cheb max jf ðnþ1Þ ðxÞj max ðx xi Þ; x ðn þ 1Þ! x i¼0 Y n 1 max ðx zi Þ; best min jf ðnþ1Þ ðxÞj x ðn þ 1Þ! x i¼0
or (when 0 2 ð; Þ of course) b~n ðxÞ :¼ fð0Þ þ xm BestApproxnm ðf0 ÞðxÞ;
ð27Þ
c~n ðxÞ :¼ fð0Þ þ xm ChebApproxnm ðf0 ÞðxÞ;
ð28Þ
which preserve fð0Þ exactly. Define absolute approximation errors best :¼ max jfðxÞ bn ðxÞj; x
cheb :¼ max jfðxÞ cn ðxÞj; x
and relative approximation errors fðxÞ bn ðxÞ ; best :¼ max x fðxÞ fðxÞ cn ðxÞ ; cheb :¼ max x fðxÞ and similarly define errors ~best , ~cheb , ~best , and ~cheb when bn ðxÞ and cn ðxÞ are replaced by b~n ðxÞ and c~n ðxÞ, respectively. For any function hðxÞ on ½; , define min jhðxÞj: ðhÞ :¼ max jhðxÞj x
x
Our theorems below are only useful when ðhÞ is around 1 or at least not too big for the h of interest. The worst case is when h has a zero in ½; and then ðhÞ ¼ þ1. Fortunately, in elementary function computations, the interval ½; is kept relatively narrow by choice and h varies little and stays away from zero and, in fact, often ðhÞ 1:2, as will be detailed in Section 7. Our theorems show the ratio of cheb over best lies in a tiny neighborhood of 1 unless fð0Þ ¼ 0, for which case the ratio of ~cheb over ~best enjoys the same property. This is significant because, as we mentioned in the introduction, it shows that simple Chebyshev interpolation is sufficient. It is also argued that, in the case of b~n or c~n , good approximations to f0 ðxÞ will generally yield good approximations to fðxÞ.
where the second inequality holds because of jwðxÞjjhðxÞj ðmin jwðxÞjÞjhðxÞj x
) max jwðxÞjjhðxÞj min jwðxÞj max jhðxÞj x
x
x
for two arbitrary functions w and h. The rightmost inequality in (29) is a consequence of the above two inequalities and Theorem 3.2. The inequalities in (30) are consequences of those in (29) and cheb cheb cheb ; maxx jfðxÞj minx jfðxÞj best best best : minx jfðxÞj maxx jfðxÞj t u Remark 5.1. For (29) to be of any use at all, ðf ðnþ1Þ Þ must not be too big and, for (30) to be useful, both ðf ðnþ1Þ Þ and ðfÞ must not be too big, which, if < 0 < , excludes the case fð0Þ ¼ 0. Although cheb best always by (29), there is no guarantee that cheb best by (30), except the two will be close under certain circumstances. For example, expðxÞ on ½; with ¼ 25 and n ¼ 5, cheb 244:453 < best 244:447 ; ~cheb 243:453 < ~best 243:447 : Odd (or even) functions are special. By Theorems 3.4 and 4.2, along the same lines of the above proof, we have the following theorem: Theorem 5.2. Let bn ðxÞ and cn ðxÞ be defined by (25) and (26). Suppose ¼ , f is even (odd), and n is even (odd). Then, 1
cheb ðf ðnþ2Þ Þ: best
ð31Þ
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
best0 :¼ max jf0 ðxÞ Bestnm ðf0 ÞðxÞj;
For even f, only and fð0Þ 6¼ 0, 1 cheb ðfÞ ðf ðnþ2Þ Þ: ðfÞ best
683
x
cheb0 :¼ max jf0 ðxÞ Chebnm ðf0 ÞðxÞj;
ð32Þ
5.2 Approximations by b~n ðxÞ and c~n ðxÞ It is natural to develop some version of Theorem 5.1 for b~n ðxÞ and c~n ðxÞ defined by (27) and (28). It is actually very hard to bound one of ~cheb and ~best by the other or one of ~cheb and ~best by the other. The reason for this is that, now, approximations to fðxÞ are obtained indirectly through f0 ðxÞ. First, we will still have max jf0 ðxÞ Bestnm ðf0 ÞðxÞj
x
x
and fðxÞ b~n ðxÞ ; best0 :¼ max x xm f0 ðxÞ fðxÞ c~n ðxÞ : cheb0 :¼ max x xm f ðxÞ 0
Theorem 5.3. Let b~n ðxÞ and c~n ðxÞ be defined by (27) and (28). We have 1
max jf0 ðxÞ Chebnm ðf0 ÞðxÞj;
cheb0 ðf ðnþ1Þ Þ; best0
ð33Þ
x
but this does not imply that one of ~cheb and ~best is always smaller than the other. Second, we also still have nm nm Y Y ðx x~i Þ max ðx z~i Þ; max x x i¼0 i¼0 but this does not imply that one of nm nm mY mY max x ðx x~i Þ and max x ðx z~i Þ x x i¼0 i¼0 is always smaller than the other, either, where the x~i s are the translated zeros of Tnmþ1 . Nevertheless, the following analysis indicates that a good approximation to f0 ðxÞ would yield a good approximation to fðxÞ. Let’s say that we have an approximation hðxÞ to f0 ðxÞ and hðxÞ ¼ f0 ðxÞ½1 þ ðxÞ; where ðxÞ is the relative error (and is bounded by a tiny number). Then, fð0Þ þ xm hðxÞ ¼ fð0Þ þ xm f0 ðxÞð1 þ ðxÞÞ ¼ fð0Þ þ xm f0 ðxÞ þ xm f0 ðxÞ ðxÞ xm f0 ðxÞ ðxÞ : ¼ fðxÞ 1 þ fðxÞ In the case fð0Þ ¼ 0, the induced relative error to fðxÞ is also ðxÞ and, thus, an approximation to f0 ðxÞ with a smaller relative error maxx jðxÞj leads to one with a smaller relative error for fðxÞ as well. However, when fð0Þ 6¼ 0, a smaller maxx jðxÞj does not necessarily imply a smaller m x f0 ðxÞ ðxÞ: max x fðxÞ But, since m m x f0 ðxÞ x f0 ðxÞ max jðxÞj; ðxÞ max max x x fðxÞ fðxÞ x unless the two sides of the above inequality differ dramatically, we still expect that a small relative error maxx jðxÞj would ultimately indicate a good approximation to fðxÞ. Set
1 cheb0 ðf0 Þ ðf ðnþ1Þ Þ: ðf0 Þ best0
ð34Þ
If fð0Þ ¼ 0, then (34) becomes 1 ~cheb ðf0 Þ ðf ðnþ1Þ Þ: ðf0 Þ ~best
ð35Þ
Proof. Equations (33) and (34) can be proven along the same lines as for Theorem 5.1, using (5) and (23). For (35), we notice that ~best ¼ best0 and ~cheb ¼ cheb0 for both the best and the Chebyshev polynomials when fð0Þ ¼ 0. u t As the counterpart to Theorem 5.2, we have: Theorem 5.4. Let b~n ðxÞ and c~n ðxÞ be defined by (27) and (28). Suppose ¼ , f is odd (even), and n is odd (even). Then, cheb0 ðf ðnþ2Þ Þ; best0
ð36Þ
1 cheb0 ðf0 Þ ðf ðnþ2Þ Þ: ðf0 Þ best0
ð37Þ
1
For the case fð0Þ ¼ 0, (37) becomes 1 ~cheb ðf0 Þ ðf ðnþ2Þ Þ: ðf0 Þ ~best
ð38Þ
5.3 Powell’s Theorem Using an entirely different technique, Powell [6] showed that n cheb 1 X ði þ 1=2Þ : ð39Þ n :¼ 1 þ tan 1 n þ 1 i¼0 2ðn þ 1Þ best Powell’s technique is quite elegant and needs no assumption on f being differentiable as ours does. Powell also showed that n grows rather slowly as n grows. In fact, 2:414 n 4:037
for 1 n 25:
Analogously to (30), (39) yields 1 cheb n ðfÞ: ðfÞ best
ð40Þ
684
IEEE TRANSACTIONS ON COMPUTERS,
Since 1:27 log2 ð n Þ 2:02 for 1 n 25, (40) is translated into differences of no more than 2 to 3 bits in accuracy if ðfÞ is nearly 1; see Section 6 for more analysis. Similar conclusions on approximations to f0 ðxÞ can be derived, but we omit the detail.
6
AN ASYMPTOTIC ANALYSIS
We shall now perform an asymptotic analysis to explain the implications of our main theorems in the context of elementary function computations, which is the original motivation for our work. In this application, usually
VOL. 53,
NO. 6,
JUNE 2004
2 f 0 ð0Þ < < 2 f 0 ð0Þ 2 f ðnþ2Þ ð0Þ ; þ ln 2 fð0Þ ln 2 fð0Þ ln 2 f ðnþ1Þ ð0Þ
ðmþ1Þ ðmþ1Þ f f 2 ð0Þ < 2 ð0Þ < 0 ðmÞ ðmÞ ðm þ 1Þ ln 2 f ð0Þ ðm þ 1Þ ln 2 f ð0Þ 2 f ðnþ2Þ ð0Þ þ ; ln 2 f ðnþ1Þ ð0Þ ð49Þ
and, if fð0Þ ¼ 0,
¼ ¼ > 0; < 1 is tiny; e:g:; 25 or smaller: ð41Þ This is especially true for table-based methods. Notice that
ðmþ1Þ ðmþ1Þ f f 2 ð0Þ < < 2 ð0Þ ~ ðmÞ ðmÞ ðm þ 1Þ ln 2 f ð0Þ ðm þ 1Þ ln 2 f ð0Þ 2 f ðnþ2Þ ð0Þ þ : ln 2 f ðnþ1Þ ð0Þ ð50Þ
log2 jðrel: err:Þj ¼ ðNo: of accurate bits in an approx:Þ: We are interested in knowing how big the differences, :¼ log2 cheb log2 best ;
ð42Þ
0 :¼ log2 cheb0 log2 best0 ;
ð43Þ
~ :¼ log2 ~cheb log2 ~best ;
ð44Þ
can get. As we shall see later, these differences are on the order of OðÞ, in general, and Oð2 Þ for even or odd functions. It follows from Theorem 5.1 and Theorem 5.3 that log2 ðfÞ log2 ðfÞ þ log2 ðf ðnþ1Þ Þ;
ð45Þ
log2 ðf0 Þ 0 log2 ðf0 Þ þ log2 ðf ðnþ1Þ Þ;
ð46Þ
In elementary function computations, usually these upper bounds are less than 1 and lower bounds bigger than 1 and, thus, the differences in the numbers of correct bits are under 1 bit! The bounds in (48)-(50) are of order OðÞ. This can be improved to order Oð2 Þ for odd (or even) function f. To see this, it suffices to show, under appropriate conditions, the s in both Theorems 5.2 and 5.4 are in a neighborhood of width Oð2 Þ of 1. For the purpose of showing this idea, we take Theorem 5.4 as an example. Under the conditions of Theorem 5.4, we have f0 ðxÞ ¼
f ðmÞ ð0Þ f ðmþ2Þ ð0Þ þ x2 þ Oðx4 Þ; m! ðm þ 2Þ!
f ðnþ2Þ ðxÞ ¼ f ðnþ2Þ ð0Þ þ x2
and, if fð0Þ ¼ 0, log2 ðf0 Þ ~ log2 ðf0 Þ þ log2 ðf ðnþ1Þ Þ:
ð47Þ
What do these inequalities tell us? Suppose fðxÞ is wellbehaved in ½; , as are most elementary mathematical functions. Since is very small, we have4 fðxÞ ¼ fð0Þ þ xf 0 ð0Þ þ Oðx2 Þ; f0 ðxÞ ¼
f ðmÞ ð0Þ f ðmþ1Þ ð0Þ þx þ Oðx2 Þ; m! ðm þ 1Þ!
f ðnþ1Þ ðxÞ ¼ f ðnþ1Þ ð0Þ þ xf ðnþ2Þ ð0Þ þ Oðx2 Þ; which imply 0 f ð0Þ þ Oð2 Þ; ðfÞ 1 þ 2 fð0Þ 2 f ðmþ1Þ ð0Þ þ Oð2 Þ; ðf0 Þ 1 þ m þ 1 f ðmÞ ð0Þ ðnþ2Þ f ð0Þ ðnþ1Þ ðf Þ 1 þ 2 ðnþ1Þ þ Oð2 Þ: ð0Þ f Therefore, ignoring the terms of order Oð2 Þ or higher, we have that 4. For f0 ðxÞ, we notice that m is the algebraic multiplicity of x ¼ 0 as a zero of fðxÞ fð0Þ.
ð48Þ
f ðnþ4Þ ð0Þ þ Oðx4 Þ; 2!
since, for odd (even) functions, even (odd) derivatives at x ¼ 0 are zero. Therefore, ðmþ2Þ f 22 ð0Þ þ Oð4 Þ; ðf0 Þ 1 þ ðmÞ ðm þ 1Þðm þ 2Þ f ð0Þ ðnþ4Þ f ð0Þ ðf ðnþ2Þ Þ 1 þ 2 ðnþ2Þ þ Oð4 Þ: ð0Þ f
7
APPLICATIONS TO ELEMENTARY FUNCTION COMPUTATIONS
In this section, we shall see the implications of the bounds we developed in the previous sections on the most common elementary functions—exponentials, logarithms, and trigonometric functions. Assume (41).
7.1 Functions ax , ax 1 for a > 1 Here, fðxÞ ¼ ax or ax 1, where a > 1, and, then, f0 ðxÞ ¼ ðax 1Þ=x. This includes all natural exponential functions as special cases. We have f 0 ðxÞ ¼ ax ln a, f ðnþ1Þ ðxÞ ¼ ax ðln aÞnþ1 , and f00 ðxÞ ¼
xax ln a ðax 1Þ : x2
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
685
TABLE 1 Bounds on Differences in Numbers of Accurate Bits to ax
It can be proven that f0 ðxÞ > 0 and f00 ðxÞ > 0, using their Maclaurin series. Therefore, for fðxÞ ¼ ax ðfÞ ¼ a2 ¼ ðf ðnþ1Þ Þ;
ðf0 Þ ¼ f0 ðÞ=f0 ðÞ ¼ a :
We have, by (45) and (46) 2 log2 a 4 log2 a;
ð51Þ
log2 a 0 3 log2 a;
ð52Þ
and, if fðxÞ ¼ ax 1 for which f0 ðxÞ is the same as above, polynomials b~n ðxÞ or c~n ðxÞ should be used, by (47) log2 a ~ 3 log2 a:
ð53Þ
Table 1 displays 4 log2 a as an upper bound on the differences in accurate bits by the two approximations as a and ¼ 2N vary.
7.2 Function loga ð1 þ xÞ for a > 1 Here, fðxÞ ¼ loga ð1 þ xÞ for a > 1. This includes the natural logarithm as a special case. Since fð0Þ ¼ 0, f0 ðxÞ ¼ loga ð1 þ xÞ=x and approximations b~n ðxÞ or c~n ðxÞ should be used. We have f 0 ðxÞ ¼ ð1 þ xÞ1 = ln a, and n! ð1 þ xÞn1 ; f ðnþ1Þ ðxÞ ¼ ð1Þ ln a x=ð1 þ xÞ lnð1 þ xÞ : f00 ðxÞ ¼ x2 ln a n
It can be seen that5 for jxj < 1, f0 ðxÞ > 0 and f00 ðxÞ < 0, using their Maclaurin series. Therefore, f0 ðÞ log ð1 Þ lnð1 Þ ; ¼ a f0 ðÞ loga ð1 þ Þ lnð1 þ Þ 1 þ nþ1 ðf ðnþ1Þ Þ ¼ : 1 ðf0 Þ ¼
We have bounds lnð1 Þ lnð1 Þ ~ log2 log2 lnð1 þ Þ lnð1 þ Þ 1þ : þ ðn þ 1Þ log2 1
7.3 Function sin x Here, fðxÞ ¼ sin x. Assume jxj < =4 or smaller and n is odd. Since fð0Þ ¼ 0, f0 ðxÞ ¼ sin x=x and approximations b~n ðxÞ or c~n ðxÞ should be used. We have f 0 ðxÞ ¼ cos x, f ð2kÞ ðxÞ ¼ ð1Þk sin x, f ð2kþ1Þ ðxÞ ¼ ð1Þk cos x, and f00 ðxÞ ¼
x cos x sin x : x2
It can be seen that for f0 ðxÞ > 0, and that f00 ðxÞ > 0 for x < 0 and f00 ðxÞ < 0 for 0 < x, using their Maclaurin series. Therefore, ðf0 Þ ¼
f0 ð0Þ ; ¼ f0 ðÞ sin
ðf ðnþ2Þ Þ ¼
1 : cos
By (47), we have bounds ~ log2 log2 ðcos Þ: log2 sin sin
ð55Þ
The bounds in (55), as well as those in (59) and (60) for cos x, are independent of n. Table 3 displays largest bounds log2 ðcos Þ; 1 :¼ log2 ð56Þ sin 2 :¼ 2 log2 cos ; 3 :¼ log2
2 log2 cos 2ð1 cos Þ
ð57Þ ð58Þ
in (55), and (59) and (60) below as upper bounds on the difference in accurate bits by the two approximations as ¼ 2N varies for both sin x and cos x.
ð54Þ
Table 2 displays the last quantity in (54) as an upper bound on the differences in accurate bits by the two approximations as n and ¼ 2N vary. More than 1 bit differences may occur only for ¼ 25 with n ¼ 10 or higher. Usually, is made much smaller than this. Current HP-UX implementation [12] and 5. They are actually true for any x > 1.
Intel’s implementation [1] for loga for a ¼ 2; e; 10 on Itanium processors [16] uses 28 .
7.4 Function cos x Now, fðxÞ ¼ cos x and then f0 ðxÞ ¼ ðcos x 1Þ=x2 . Assume jxj < =4 or smaller, and n is even. We have f 0 ðxÞ ¼ sin x, f ð2k1Þ ðxÞ ¼ ð1Þk sin x, f ð2kÞ ðxÞ ¼ ð1Þk cos x, and f00 ðxÞ ¼
x sin x 2 cos x þ 2 : x3
It can be proven that f0 ðxÞ < 0 and xf00 ðxÞ 0; so, f0 ðxÞ is an increasing function for x > 0 and a decreasing one for x < 0. Therefore,
686
IEEE TRANSACTIONS ON COMPUTERS,
VOL. 53,
NO. 6,
JUNE 2004
TABLE 2 Bounds on Differences in Numbers of Accurate Bits to loga x
TABLE 3 Bounds on Differences in Numbers of Accurate Bits to sin x and cos x
ðfÞ ¼
1 ¼ ðf ðnþ2Þ Þ; cos
ðf0 Þ ¼
f0 ð0Þ 2 ¼ : f0 ðÞ 2ð1 cos Þ
min jf ðnþ2Þ ðxÞj ¼ f ðnþ2Þ ð0Þ jxj
¼ ðn þ 2Þ!
We have, by (45) and (46), log2 cos 2 log2 cos ; log2
ð59Þ
2 2 0 log2 log2 cos : ð60Þ 2ð1 cos Þ 2ð1 cos Þ
Table 3 displays the upper bounds in (59) and (60) as 2 and 3 .
7.5 Function arcsin x Here, fðxÞ ¼ arcsin x for jxj 1=2 and, then, f0 ðxÞ ¼ arcsin x=x. Assume n is odd. The Maclaurin series of fðxÞ is [17, p. 81]
n!! ¼ ðn!!Þ2 ; ðn þ 1Þ!! ðn þ 2Þ
minjxj f0 ðxÞ ¼ f0 ð0Þ ¼ 1, and maxjxj f0 ðxÞ ¼ f0 ðÞ. Therefore, ðf ðnþ2Þ Þ ¼
1 2
ðn!!Þ
f ðnþ2Þ ðÞ;
ðf0 Þ ¼
f0 ðÞ arcsin ; ¼ f0 ð0Þ
which imply log2
arcsin arcsin ~ log2 1 f ðnþ2Þ ðÞ: þ log2 ðn!!Þ2
ð61Þ
Table 4 displays the last quantity in (61) as an upper bound 3
arcsin x ¼ x þ
5
2kþ1
1x 1 3x ð2k 1Þ!! x þ þ þ þ : 2 3 24 5 ð2kÞ!! 2k þ 1
Thus, maxjxj jf ðnþ2Þ ðxÞj ¼ f ðnþ2Þ ðÞ and
on the differences in accurate bits by the two approximations as n and ¼ 2N vary.
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
687
TABLE 4 Bounds on Differences in Numbers of Accurate Bits to arcsin x
8
CONCLUSIONS
[7]
This paper shows that the Chebyshev interpolating polynomials are just as good as the best polynomials in the context of table-driven methods for elementary function computations and, thus, the common practice of feeding Chebyshev interpolating polynomials into a Remez procedure offers negligible improvement, if any, and thus may be considered unnecessary. Our results improve the earlier results due to Powell [6] which, however, have more general applications.
[11]
ACKNOWLEDGMENTS
[12]
The author is grateful to Prof. W. Kahan of the University of California at Berkeley who brought Powell [6] to his attention and to Professor M.J.D. Powell of the University of Cambridge for sending the paper to him. He thanks the referees for their careful reading of the manuscript and constructive suggestions that improved the presentation. Part of this work was done while the author was on leave at Hewlett-Packard Company. He is grateful for help received from Jim Thomas, Jon Okada, and Peter Markstein of the HP Itanium floating-point and elementary math library team at Cupertino, California. This work was supported in part by the US National Science Foundation CAREER award under Grant No. CCR-9875201.
REFERENCES [1] [2] [3] [4] [5] [6]
J. Harrison, T. Kubaska, S. Story, and P.T.P. Tang, “The Computation of Transcendental Functions on the IA-64 Architecture,” Intel Technology J., no. Q4, pp. 1-7, Nov. 1999. R.-C. Li, S. Boldo, and M. Daumas, “Theorems on Efficient Argument Reductions,” Proc. 16th IEEE Symp. Computer Arithmetic, June 2003. P. Markstein, IA-64 and Elementary Functions: Speed and Precision, 2000. J.-M. Muller, Elementary Functions: Algorithms and Implementation. Boston, Basel, Berlin: Birkha˚user, 1997. K.C. Ng, “Argument Reduction for Huge Arguments: Good to the Last Bit,” SunPro, technical report, 1992, http://www.validlab. com/. M.J.D. Powell, “On the Maximum Errors of Polynomial Approximations Defined by Interpolation and by Least Squares Criteria,” The Computer J., vol. 9, no. 4, pp. 404-407, Feb. 1967.
[8] [9]
[10]
[13] [14] [15] [16]
[17]
P.T.P. Tang, “Table-Driven Implementation of the Exponential Function in IEEE Floating-Point Arithmetic,” ACM Trans. Math. Software, vol. 15, no. 2, pp. 144-157, June 1989. P.T.P. Tang, “Table-Driven Implementation of the Logarithm Function in IEEE Floating-Point Arithmetic,” ACM Trans. Math. Software, vol. 16, no. 4, pp. 378-400, Dec. 1990. P.T.P. Tang, “Table Lookup Algorithms for Elementary Functions and Their Error Analysis,” Proc. 10th IEEE Symp. Computer Arithmetic, P. Kornerup and D.W. Matula, eds., pp. 232-236, June 1991. P.T.P. Tang, “Table-Driven Implementation of the expm1 Function in IEEE Floating-Point Arithmetic,” ACM Trans. Math. Software, vol. 18, no. 2, pp. 211-222, June 1992. E. Remez, “Sur un Proce´de´ Convergent d’Approximations Successives pour De´terminer les Polynoˆmes d’Approximation,” C.R. Acade´mie des Sciences, Paris, vol. 198, 1934. R.-C. Li, P. Markstein, J. Okada, and J. Thomas, “The libm Library and Floating-Point Arithmetic in HP-UX for Itanium II,” June 2000, http://h21007.www2.hp.com/dspp/files/unprotected/ Itanium/FP_White_Paper_v2.pdf. E.W. Cheney, Approximation Theory, second ed. Providence, R.I.: AMS, 1981. R.A. DeVore and G.G. Lorentz, Constructive Approximation. Berlin, Heidelberg, New York: Springer-Verlag, 1993. G.G. Lorentz, Approximation of Functions, second ed. New York: Chelsea Publishing, 1986. INTEL, Intel IA-64 Architecture Software Developer’s Manual. Intel Corp., 2002, vol. 3: Instruction Set Reference, document No. 245319-002, http://developer.intel.com/design/itanium/ manuals/index.htm. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, M. Abramowitz and I.A Stegun, eds., ninth ed. New York: Dover Publications, 1970.
Ren-Cang Li received the BS degree in computational mathematics from Xiamen University, People’s Republic of China, in 1985, the MS degree, also in computational mathematics, from the Chinese Academy of Science in 1988, and the PhD degree in applied mathematics from the University of California at Berkeley in 1995. He is currently an associate professor in the Department of Mathematics, University of Kentucky, Lexington. He was awarded the 1995 Householder Fellowship in Scientific Computing by Oak Ridge National Laboratory, a Friedman memorial prize in Applied Mathematics from the University of California at Berkeley in 1996, and a CAREER award from the US National Science Foundation in 1999. His research interest includes elementary function computations, floating-point support for scientific computing, large and sparse linear systems, eigenvalue problems, and model reduction, and unconventional schemes for differential equations.
IEEE TRANSACTIONS ON COMPUTERS,
VOL. 53,
NO. 6,
JUNE 2004
Near Optimality of Chebyshev Interpolation for Elementary Function Computations Ren-Cang Li Abstract—A common practice for computing an elementary transcendental function in an libm implementation nowadays has two phases: reductions of input arguments to fall into a tiny interval and polynomial approximations for the function within the interval. Typically, the interval is made tiny enough so that polynomials of very high degree aren’t required for accurate approximations. Often, approximating polynomials as such are taken to be the best polynomials or any others such as the Chebyshev interpolating polynomials. The best polynomial of degree n has the property that the biggest difference between it and the function is smallest among all possible polynomials of degrees no higher than n. Thus, it is natural to choose the best polynomials over others. In this paper, it is proven that the best polynomial can only be more accurate by at most a fractional bit than the Chebyshev interpolating polynomial of the same degree in computing elementary functions or, in other words, the Chebyshev interpolating polynomials will do just as well as the best polynomials. Similar results were obtained in 1967 by Powell who, however, did not target elementary function computations in particular and placed no assumption on the function and, remarkably, whose results imply accuracy differences of no more than 2 to 3 bits in the context of this paper. Index Terms—Elementary function computation, libm, Chebyshev Interpolation, Remez, best polynomial, accuracy.
æ 1
INTRODUCTION
A
common practice for computing elementary transcendental functions fðxÞ has two phases (see [1], [2], [3], [4], [5], [7], [8], [9], [10], and references therein): An argument reduction that typically reduces x to fall into a very narrow interval ½; , e.g., 25 or smaller. 2. An efficient polynomial approximation1 to fðxÞ for x 2 ½; . Often, approximating polynomials are forced to match fðxÞ at x ¼ 0 for computational advantages. This is especially true when fð0Þ ¼ 0. Small ensures the use of approximating polynomials of low degrees and, consequently, little time to evaluate them. Often, is very tiny, especially for table-driven algorithms [3], [4], [5], [7], [8], [9], [10]. There are many ways to find approximating polynomials. Given the degree n of the polynomials, common ones are 1) finding the best polynomial pn ðxÞ that minimizes maxx2½; jfðxÞ pn ðxÞj (or a variation of this if fð0Þ is to be preserved; see below for details) and 2) finding the Chebyshev interpolating polynomial. Rational approximation is another alternative, but it is seldom used in elementary function computations because division is much slower than addition/multiplication and, moreover, the division can only be initiated after the availability of the denominator and/or the numerator. 1.
1. There may be exceptions to this, e.g., computing sin x may require polynomial approximations to both sin x and cos x, depending on algorithms.
. The author is with the Department of Mathematics, University of Kentucky, Lexington, KY 40506. E-mail: [email protected]. Manuscript received 3 Apr. 2003; revised 6 Nov. 2003; accepted 17 Nov. 2003. For information on obtaining reprints of this article, please send e-mail to: [email protected], and reference IEEECS Log Number 118532. 0018-9340/04/$20.00 ß 2004 IEEE
This paper compares the approximation accuracy achieved by the two polynomial approximation techniques. Earlier results along this line are due to Powell [6]. However, Powell’s paper went largely unnoticed in the literature of elementary function computations. Even so, through practice it has become folklore [3], [4] that the Chebyshev interpolating polynomial will be just as effective as the best polynomial with the same degree. In fact, most error analysis for the polynomial solution part in Markstein [3] is based on the Chebyshev polynomial interpolation. In this paper, we will show that the Chebyshev interpolating polynomial may be less accurate, if indeed, than the best polynomial of the same degree by only a fractional bit. The best polynomials are typically obtained through a slowly convergent iterative procedure—the Remez procedure [11], while we do have a close formula for the Chebyshev interpolating polynomials. This does not, however, necessarily place the best polynomial approach at a disadvantage since all approximating polynomials are precomputed in any libm implementations [1], [12] and, thus, the cost in computing approximating polynomials becomes somewhat irrelevant. Even so, we feel that our contribution here has its significance, that is, that simple Chebyshev interpolation is sufficient and that the common practice of feeding Chebyshev interpolation polynomials into the Remez procedure offers negligible improvement, if any, and thus may be considered unnecessary. Our results compare favorably to Powell’s results which, applied to the context here, say the difference in accuracy is no more than 2 to 3 bits. It is worth mentioning that Powell’s results were obtained, remarkably, without any assumption on fðxÞ. The rest of this paper is organized as follows: Section 2 briefly introduces the general interpolation theory, enough to serve our purpose later on. Chebyshev interpolation and Published by the IEEE Computer Society
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
best polynomial approximations are discussed in Sections 3 and 4 with emphasis on error expansions on these approximations. Our main theorems that compare the accuracies by the two approaches are given in Section 5, where a theorem of Powell’s will also be discussed. An asymptotic analysis in Section 6 on the main theorems confirms our claim made at the beginning of this paper and concludes that the approximation accuracies between the best polynomials and the Chebyshev interpolating polynomials differ by Oð2 Þ bits for odd and even functions and by OðÞ for others. Finally, in Section 8, the main theorems are applied to the most common elementary transcendental functions, such as expðxÞ, ln x, sin x, and so on.
1.1 Notation Throughout the remainder of this paper, the interval of interest is ½; on which fðxÞ is defined. It is assumed that fðxÞ is continuous and has continuous derivatives up to ðn þ 1Þth (or more as deemed necessary), where n is the degree of polynomials sought to approximate fðxÞ. If needed, 0 2 ð; Þ is also implied. Define f0 ðxÞ :¼
fðxÞ fð0Þ ; xm
ð1Þ
where m is the multiplicity of x ¼ 0 as a zero of fðxÞ fð0Þ, and “:=” is the assignment operator. Usually, m ¼ 1, but may not be, e.g., for cos x it is 2. Unknown points ; 2 ½; in error equations later may be different from one occurrence to another.
2
GENERAL POLYNOMIAL INTERPOLATION
At times, polynomials that preserve certain properties of fðxÞ are called for so that simple but important mathematical identities, such as log 1 ¼ 0, expð0Þ ¼ 1, sin 0 ¼ 0, and cosð0Þ ¼ 1, are numerically conserved. In elementary function computations, this typically happens when 0 2 ½; (and ¼ as well, but we shall not assume this for now). To make sure that identities as such remain valid in the face of roundoffs, we need to place2 x ¼ 0 to the list of interpolating points to ensure fð0Þ ¼ pn ð0Þ exactly! One way to implement this is to interpolate f0 ðxÞ at n m þ 1 other points, x~0 ; x~1 ; . . . ; x~nm ;
ð2Þ
be interpolating points in the interval. The interpolating polynomial pn ðxÞ is uniquely determined by pn ðxÞ and fðxÞ coinciding at those points in the sense that, if xi is repeated k times,
p~n ðxÞ :¼ fð0Þ þ xm hnm ðxÞ; a polynomial of degree n that can be seen to coincide with fðxÞ at x ¼ 0, as well as at x~i s. An application of (3) gives fðxÞ p~n ðxÞ ¼
ð5Þ
for some minfx; 0; x~0 ; . . . ; x~nm g maxfx; 0; x~0 ; . . . ; x~nm g:
3
CHEBYSHEV INTERPOLATION
The standard Chebyshev polynomial interpolation for a function uðtÞ is on the closed interval ½1; 1. The idea is to seek a polynomial qn ðtÞ of degree n that matches exactly uðtÞ at the zeros ti :¼ cos
2i þ 1 ; nþ1 2
i ¼ 0; 1; . . . ; n;
ð6Þ
of Tnþ1 ðtÞ cosðn þ 1Þ arccos t, the ðn þ 1Þth Chebyshev polynomial, i.e., uðti Þ ¼ qn ðti Þ for i ¼ 0; 1; . . . ; n. It can be proven that (see, e.g., [3, p. 99]) qn ðtÞ ¼
n X
k Tk ðtÞ;
ð7Þ
where
(So, if k ¼ 1, only the function value matches at xi .) If none of the points xi collapse, this is called Lagrange Interpolation. In the extreme case when all xi s collapse to one point, pn ðxÞ is just the truncated Taylor series of fðxÞ at the collapsed point. In any case, the difference of fðxÞ and pn ðxÞ is given as follows: If f ðnþ1Þ ðxÞ exists and is continuous, then n f ðnþ1Þ ðÞ Y ðx xi Þ; ðn þ 1Þ! i¼0
Y f ðnþ1Þ ðÞ m nm x ðx x~i Þ ðn þ 1Þ! i¼0
k¼0
ðjÞ pðjÞ n ðxi Þ ¼ f ðxi Þ for j ¼ 0; . . . k 1:
fðxÞ pn ðxÞ ¼
ð4Þ
by a polynomial hnm ðxÞ of degree n m. Then, set
Polynomial interpolation is an ancient idea that is still widely in use today to approximate a function. Let x0 ; x1 ; . . . ; xn
679
ð3Þ
for some minfx; x0 ; . . . ; xn g maxfx; x0 ; . . . ; xn g [13, p. 60]. This equation offers one way to find a good polynomial to make max jfðxÞ pn ðxÞj small. Since f ðnþ1Þ ðÞ is function-dependent and, often, Q is unknown, it is natural to look for xi s so that max j ni¼0 ðx xi Þj is minimized. It turns out this minimization problem can be solved elegantly: Pick xi s as the translated zeros of the ðn þ 1Þth Chebyshev polynomial (see Theorem 3.2 below).
k ¼
Pn 1 nþ1 Pi¼0 n 2 i¼0 nþ1
uðti Þ; if k ¼ 0; Tk ðti Þuðti Þ; otherwise:
ð8Þ
Theorem 3.1. If u is odd (even), then qn is also odd (even). Proof. This is a well-known property of Chebyshev interpolation. But, we shall give a proof here for completeness. Notice that Tnþ1 ðtÞ ¼ cosðn þ 1Þ arccosðtÞ ¼ cos½ðn þ 1Þð arccos tÞ ¼ ð1Þnþ1 Tnþ1 ðtÞ: 2. Mathematically, it does not have to be 0 and, in fact, we can preserve the function value at any given floating-point number so long as the function value is also a floating-point number. Such a generalization, however, is no more general than 0 itself because shifting of origin can always make it 0!
680
IEEE TRANSACTIONS ON COMPUTERS,
Tnþ1 is an even function if n þ 1 is even and an odd function otherwise. Therefore, ti and ti appear in pairs in the list (6). Suppose u is odd, i.e., uðtÞ ¼ uðtÞ. Let sn ðtÞ :¼ qn ðtÞ, also a polynomial of degree n. We have sn ðti Þ ¼ qn ðti Þ ¼ uðti Þ ¼ uðti Þ; which says that sn also interpolates u at the same points (6) as qn does. Since the interpolating polynomial is unique, we conclude that qn ðtÞ sn ðtÞ ¼ qn ðtÞ, i.e., qn is odd. Now, suppose u is even, i.e., uðtÞ ¼ uðtÞ. Let sn ðtÞ :¼ qn ðtÞ, a polynomial of degree n. We have sn ðti Þ ¼ qn ðti Þ ¼ uðti Þ ¼ uðti Þ; which says that sn also interpolates u at the same points (6) as qn does. The uniqueness of the interpolating polynomial leads to qn ðtÞ sn ðtÞ ¼ qn ðtÞ, i.e., qn is u t
even.
Now, we turn to the continuous function fðxÞ defined on an arbitrary closed interval ½; . The idea is first to devise a transformation x ¼ xðtÞ such that gðtÞ :¼ fðxðtÞÞ is defined on ½1; 1. One such commonly used transformation is x
: ½1; 1 ! t !
½; þ xðtÞ :¼ 2 tþ 2 :
The inverse transformation of xðtÞ is 2 þ x : tðxÞ :¼ 2
NO. 6,
JUNE 2004
n n Y nþ1 Y 2 þ x ti ðx xi Þ ¼ 2 2 i¼0 i¼0 nþ1 1 Tnþ1 ðtÞ; ¼ n 2 2 ð13Þ where the last equality holds because ti s are all the zeros of Tnþ1 whose leading coefficient is 2n , t is as defined by (10). Similarly, n Y nþ1 ðx yi Þ ¼ qnþ1 ðtÞ; 2 i¼0 where qnþ1 ðtÞ is polynomial in t of degree n with leading coefficient 1. It is known [13, p. 61] that, on ½1; 1, Tnþ1 ðtÞ=2n takes the least maximum absolute value among all polynomials of degree n þ 1 with leading coefficient 1. Therefore, (11) holds. u t It follows from (7) and (8) that the Chebyshev interpolating polynomial cn ðxÞ of degree n for fðxÞ on ½; is given by n X 2 þ x ; ð14Þ k Tk cn ðxÞ ¼ 2 k¼0
ð9Þ
VOL. 53,
k ¼
Pn 1 nþ1 Pi¼0 n 2 i¼0 nþ1
fðxi Þ; Tk ðti Þfðxi Þ;
if k ¼ 0; otherwise:
ð15Þ
The proof of Theorem 3.2 also yields ð10Þ
Write, for i ¼ 0; 1; . . . ; n,
fðxÞ cn ðxÞ ¼
n f ðnþ1Þ ðÞ Y ðx xi Þ ðn þ 1Þ! i¼0
1 nþ1 f ðnþ1Þ ðÞ Tnþ1 ðtÞ; ¼ n 2 2 ðn þ 1Þ!
þ ti þ ; xi :¼ xðti Þ :¼ 2 2
ð16Þ
ð17Þ
called the translated zeros in ½; of Tnþ1 . The polynomial
for some 2 ½; and t is defined as in (10).
cn ðxÞ of degree n such that fðxi Þ ¼ cn ðxi Þ for all i is called
Theorem 3.3. Suppose ¼ . Then, if f is odd (even), cn is odd (even).
the Chebyshev interpolating polynomial of degree n for fðxÞ. The general theory described in Section 2 still applies here, so the error equation (7) holds with xi assigned here and for some 2 ½; . Theorem 3.2. The translated zeros xi s of Tnþ1 satisfy Y Y n n max ðx xi Þ max ðx yi Þ x x i¼0 i¼0
ð11Þ
for any y0 ; y1 ; . . . ; yn 2 ½; . Proof. This is also well-known. Again, a proof is given for completeness; also, a key equation (13) is useful later. We have þ ti x xi ¼ x 2 2 2 þ x ti ; ¼ 2 2
ð12Þ
Proof. ¼ implies xðtÞ ¼ t. If fðxÞ is odd (even) in x, fðtÞ is odd (even) in t. Hence, by Theorem 3.1, cn ðtÞ is u t odd (even) in t and, thus, cn ðxÞ is odd (even) in x. Tnþ1 ðtÞ has a zero at t ¼ 0 for even n. Thus, if ¼ and n is even, the Chebyshev interpolating polynomial cn ðxÞ automatically satisfies cn ð0Þ ¼ fð0Þ. But, it is a different story when n is odd. As we pointed out in Section 2, we can ensure that fð0Þ is always recovered exactly by looking for a Chebyshev interpolating polynomial hnm ðxÞ of degree n m for f0 ðxÞ and then setting c~n ðxÞ :¼ fð0Þ þ xm hnm ðxÞ: Notice that, now, hnm ðxÞ coincides with f0 ðxÞ at the translated zeros x~i s of Tnmþ1 . It is interesting to notice that Tnmþ1 ð0Þ ¼ 0 also for odd n m þ 1; this implies x ¼ 0 appears m þ 1 times in the interpolating points for fðxÞ when ¼ ; so, for odd n m þ 1 and ¼ , we have
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
fð0Þ ¼ c~n ð0Þ; f 0 ð0Þ ¼ c~0n ð0Þ; . . . ; f ðmÞ ð0Þ ¼ c~ðmÞ n ð0Þ: With the arguments similar to those which lead to the error equation (17), we can prove fðxÞ c~n ðxÞ fðxÞ ½fð0Þ þ xm hnm ðxÞ Y f ðnþ1Þ ðÞ m nm x ¼ ðx x~i Þ ðn þ 1Þ! i¼0
¼
1 2nm
nmþ1 2
f ðnþ1Þ ðÞ m x Tnmþ1 ðtÞ; ðn þ 1Þ!
ð18Þ
fðxÞ cn ðxÞ ¼
2nþ1
f ðnþ2Þ ðÞ Tnþ2 ðtÞ; ðn þ 2Þ!
1 nmþ2 2nmþ1 f ðnþ2Þ ðÞ m x Tnmþ2 ðtÞ; ðn þ 2Þ!
ð20Þ
fðxÞ c~n ðxÞ ¼
ð21Þ
where ; 2 ½; and t is defined as in (6).
4
BEST APPROXIMATION
Given a continuous function fðxÞ defined on a closed interval ½; , the best polynomial approximation bn ðxÞ of degree n for fðxÞ is the one such that n :¼ max jfðxÞ bn ðxÞj :¼ min max jfðxÞ pn ðxÞj; x
such that fðzi Þ ¼ bn ðzi Þ for i ¼ 0; 1; . . . ; n. So, the best polynomial may well be thought of as a Lagrange interpolation polynomial with interpolating points not known a priori and, therefore,
ð19Þ
Theorem 3.4. Assume ¼ and f is odd (even). If n is odd (even), then nþ2
By the continuity of the functions, we conclude that, between every two consecutive yi s, fðxÞ bn ðxÞ has a zero, i.e., there are y0 < z0 < . . . < yi < zi < yiþ1 < . . . < zn < ynþ1
for some 2 ½; and t is defined as in (10). It can be seen that Theorem 3.3 holds with cn replaced by c~n . There is an important implication to this fact and Theorem 3.3. Assume ¼ and n is odd (or even), depending on whether f is odd (or even). For better error bounds, we shall always think of interpolating f (or f0 ) at the translated zeros of Tnþ2 (or of Tnmþ2 ), in contrast to Tnþ1 (or Tnmþ1 ) for the general f. In general, this will produce a Chebyshev interpolating polynomial of degree n þ 1, but, since it has to be odd (or even), its leading coefficient must be zero! and thus degenerates to a polynomial of degree n.
1
681
fðxÞ bn ðxÞ ¼
n f ðnþ1Þ ðÞ Y ðx zi Þ; ðn þ 1Þ! i¼0
ð22Þ
for some 2 ½; . Such reinterpretation is usually not useful from a numerical point of view, but it is a key step in deriving our main theorems in the next section. For the same reason as explained in Sections 2 and 3, it may also be advantageous to work with f0 ðxÞ for the sake of simple but important mathematical identities. Let hnm ðxÞ be the best polynomial for f0 ðxÞ in the sense that max jf0 ðxÞ hnm ðxÞj ¼ min max jf0 ðxÞ pnm ðxÞj pnm x
x
taken over all possible polynomials pnm ðxÞ of degree n m, and set b~n ðxÞ :¼ fð0Þ þ xm hnm ðxÞ; which approximates fðxÞ over the interval. Of course, then b~n ðxÞ may do worse than3 bn ðxÞ, but the trade off is that now fð0Þ ¼ b~n ð0Þ always! Similarly, we have the following error equation Y f ðnþ1Þ ð~Þ m nm fðxÞ b~n ðxÞ ¼ ðx z~i Þ; x ðn þ 1Þ! i¼0
ð23Þ
for some ~ 2 ½; and z~0 < z~1 < . . . < z~nm . Even though, now, b~n ðxÞ is no longer the best polynomial for fðxÞ, we still call it the best polynomial at places where no confusion may rise. It can be seen that the uniqueness of the best polynomial implies that bn and b~n are odd (even) if f is odd (even) and ¼ . This implies: Theorem 4.2. Assume ¼ and f is odd (even). If n is odd (even), then
pn x
where the min is taken over all possible polynomials pn ðxÞ of degree n. It is known [13], [14, p. 58], [15] that bn ðxÞ exists and is unique and is characterized by the following theorem, traditionally called the Chebyshev Theorem, although he did not prove it (see [14, p. 58] for a detailed history of it). Theorem 4.1 (Chebyshev). Assume f is not a polynomial of degree less than or equal to n. Then, there are at least n þ 2 points y0 < y1 < . . . < ynþ1 at which fðyi Þ bn ðyi Þ ¼ sð1Þi n for 0 i n þ 1, where s ¼ 1, independent of i. In other words, fðxÞ pn ðxÞ travels from one extreme point to another of the opposite sign for at least n þ 1 times.
bn ¼ bnþ1 ;
b~n ¼ b~nþ1 :
Consequently, fðxÞ bn ðxÞ ¼
þ1 f ðnþ2Þ ðÞ nY ðx zi Þ; ðn þ 2Þ! i¼0
ð24Þ
fðxÞ b~n ðxÞ ¼
Y f ðnþ2Þ ð~Þ m nmþ1 ðx z~i Þ; x ðn þ 2Þ! i¼0
ð25Þ
where zi ; ; z~i ; ~ 2 ½; . Proof. Consider the case where n is odd and f is odd only. The other case can be proven in the same way. All the 3. It may be true the other way. See Remark 5.1.
682
IEEE TRANSACTIONS ON COMPUTERS,
best polynomials of f are odd, so bnþ1 is odd. But, since n þ 1 is even, a necessary condition for bnþ1 to be odd is its leading coefficient is zero, i.e., bnþ1 ’s degree is no higher than n. Now, the uniqueness of the best polynomials implies bn ¼ bnþ1 . b~n ¼ b~nþ1 can be proven in the same way. u t
5
VOL. 53,
NO. 6,
JUNE 2004
5.1 Approximations by bn ðxÞ and cn ðxÞ Theorem 5.1. Let bn ðxÞ and cn ðxÞ be defined by (25) and (26). We have 1
cheb ðf ðnþ1Þ Þ; best
ð29Þ
1 cheb ðfÞ ðf ðnþ1Þ Þ: ðfÞ best
MAIN THEOREMS
We shall use BestApproxn ðfÞ :best approximation polynomial; ChebApproxn ðfÞ :Chebyshev interpolating polynomial; of degree n, as building blocks to construct polynomials. Specifically, we will study two sets of approximation polynomials of degree n defined either as bn ðxÞ :¼ BestApproxn ðfÞðxÞ;
ð25Þ
cn ðxÞ :¼ ChebApproxn ðfÞðxÞ;
ð26Þ
ð30Þ
Proof. best cheb holds since bn ðxÞ is the best among all polynomials of degree n. This proves the leftmost inequality in (29). We now turn to the rightmost inequality in (29). Our proof relies on (3) with translated zeros xi s of Tnþ1 , (22), and Theorem 3.2. We have Y n 1 cheb max jf ðnþ1Þ ðxÞj max ðx xi Þ; x ðn þ 1Þ! x i¼0 Y n 1 max ðx zi Þ; best min jf ðnþ1Þ ðxÞj x ðn þ 1Þ! x i¼0
or (when 0 2 ð; Þ of course) b~n ðxÞ :¼ fð0Þ þ xm BestApproxnm ðf0 ÞðxÞ;
ð27Þ
c~n ðxÞ :¼ fð0Þ þ xm ChebApproxnm ðf0 ÞðxÞ;
ð28Þ
which preserve fð0Þ exactly. Define absolute approximation errors best :¼ max jfðxÞ bn ðxÞj; x
cheb :¼ max jfðxÞ cn ðxÞj; x
and relative approximation errors fðxÞ bn ðxÞ ; best :¼ max x fðxÞ fðxÞ cn ðxÞ ; cheb :¼ max x fðxÞ and similarly define errors ~best , ~cheb , ~best , and ~cheb when bn ðxÞ and cn ðxÞ are replaced by b~n ðxÞ and c~n ðxÞ, respectively. For any function hðxÞ on ½; , define min jhðxÞj: ðhÞ :¼ max jhðxÞj x
x
Our theorems below are only useful when ðhÞ is around 1 or at least not too big for the h of interest. The worst case is when h has a zero in ½; and then ðhÞ ¼ þ1. Fortunately, in elementary function computations, the interval ½; is kept relatively narrow by choice and h varies little and stays away from zero and, in fact, often ðhÞ 1:2, as will be detailed in Section 7. Our theorems show the ratio of cheb over best lies in a tiny neighborhood of 1 unless fð0Þ ¼ 0, for which case the ratio of ~cheb over ~best enjoys the same property. This is significant because, as we mentioned in the introduction, it shows that simple Chebyshev interpolation is sufficient. It is also argued that, in the case of b~n or c~n , good approximations to f0 ðxÞ will generally yield good approximations to fðxÞ.
where the second inequality holds because of jwðxÞjjhðxÞj ðmin jwðxÞjÞjhðxÞj x
) max jwðxÞjjhðxÞj min jwðxÞj max jhðxÞj x
x
x
for two arbitrary functions w and h. The rightmost inequality in (29) is a consequence of the above two inequalities and Theorem 3.2. The inequalities in (30) are consequences of those in (29) and cheb cheb cheb ; maxx jfðxÞj minx jfðxÞj best best best : minx jfðxÞj maxx jfðxÞj t u Remark 5.1. For (29) to be of any use at all, ðf ðnþ1Þ Þ must not be too big and, for (30) to be useful, both ðf ðnþ1Þ Þ and ðfÞ must not be too big, which, if < 0 < , excludes the case fð0Þ ¼ 0. Although cheb best always by (29), there is no guarantee that cheb best by (30), except the two will be close under certain circumstances. For example, expðxÞ on ½; with ¼ 25 and n ¼ 5, cheb 244:453 < best 244:447 ; ~cheb 243:453 < ~best 243:447 : Odd (or even) functions are special. By Theorems 3.4 and 4.2, along the same lines of the above proof, we have the following theorem: Theorem 5.2. Let bn ðxÞ and cn ðxÞ be defined by (25) and (26). Suppose ¼ , f is even (odd), and n is even (odd). Then, 1
cheb ðf ðnþ2Þ Þ: best
ð31Þ
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
best0 :¼ max jf0 ðxÞ Bestnm ðf0 ÞðxÞj;
For even f, only and fð0Þ 6¼ 0, 1 cheb ðfÞ ðf ðnþ2Þ Þ: ðfÞ best
683
x
cheb0 :¼ max jf0 ðxÞ Chebnm ðf0 ÞðxÞj;
ð32Þ
5.2 Approximations by b~n ðxÞ and c~n ðxÞ It is natural to develop some version of Theorem 5.1 for b~n ðxÞ and c~n ðxÞ defined by (27) and (28). It is actually very hard to bound one of ~cheb and ~best by the other or one of ~cheb and ~best by the other. The reason for this is that, now, approximations to fðxÞ are obtained indirectly through f0 ðxÞ. First, we will still have max jf0 ðxÞ Bestnm ðf0 ÞðxÞj
x
x
and fðxÞ b~n ðxÞ ; best0 :¼ max x xm f0 ðxÞ fðxÞ c~n ðxÞ : cheb0 :¼ max x xm f ðxÞ 0
Theorem 5.3. Let b~n ðxÞ and c~n ðxÞ be defined by (27) and (28). We have 1
max jf0 ðxÞ Chebnm ðf0 ÞðxÞj;
cheb0 ðf ðnþ1Þ Þ; best0
ð33Þ
x
but this does not imply that one of ~cheb and ~best is always smaller than the other. Second, we also still have nm nm Y Y ðx x~i Þ max ðx z~i Þ; max x x i¼0 i¼0 but this does not imply that one of nm nm mY mY max x ðx x~i Þ and max x ðx z~i Þ x x i¼0 i¼0 is always smaller than the other, either, where the x~i s are the translated zeros of Tnmþ1 . Nevertheless, the following analysis indicates that a good approximation to f0 ðxÞ would yield a good approximation to fðxÞ. Let’s say that we have an approximation hðxÞ to f0 ðxÞ and hðxÞ ¼ f0 ðxÞ½1 þ ðxÞ; where ðxÞ is the relative error (and is bounded by a tiny number). Then, fð0Þ þ xm hðxÞ ¼ fð0Þ þ xm f0 ðxÞð1 þ ðxÞÞ ¼ fð0Þ þ xm f0 ðxÞ þ xm f0 ðxÞ ðxÞ xm f0 ðxÞ ðxÞ : ¼ fðxÞ 1 þ fðxÞ In the case fð0Þ ¼ 0, the induced relative error to fðxÞ is also ðxÞ and, thus, an approximation to f0 ðxÞ with a smaller relative error maxx jðxÞj leads to one with a smaller relative error for fðxÞ as well. However, when fð0Þ 6¼ 0, a smaller maxx jðxÞj does not necessarily imply a smaller m x f0 ðxÞ ðxÞ: max x fðxÞ But, since m m x f0 ðxÞ x f0 ðxÞ max jðxÞj; ðxÞ max max x x fðxÞ fðxÞ x unless the two sides of the above inequality differ dramatically, we still expect that a small relative error maxx jðxÞj would ultimately indicate a good approximation to fðxÞ. Set
1 cheb0 ðf0 Þ ðf ðnþ1Þ Þ: ðf0 Þ best0
ð34Þ
If fð0Þ ¼ 0, then (34) becomes 1 ~cheb ðf0 Þ ðf ðnþ1Þ Þ: ðf0 Þ ~best
ð35Þ
Proof. Equations (33) and (34) can be proven along the same lines as for Theorem 5.1, using (5) and (23). For (35), we notice that ~best ¼ best0 and ~cheb ¼ cheb0 for both the best and the Chebyshev polynomials when fð0Þ ¼ 0. u t As the counterpart to Theorem 5.2, we have: Theorem 5.4. Let b~n ðxÞ and c~n ðxÞ be defined by (27) and (28). Suppose ¼ , f is odd (even), and n is odd (even). Then, cheb0 ðf ðnþ2Þ Þ; best0
ð36Þ
1 cheb0 ðf0 Þ ðf ðnþ2Þ Þ: ðf0 Þ best0
ð37Þ
1
For the case fð0Þ ¼ 0, (37) becomes 1 ~cheb ðf0 Þ ðf ðnþ2Þ Þ: ðf0 Þ ~best
ð38Þ
5.3 Powell’s Theorem Using an entirely different technique, Powell [6] showed that n cheb 1 X ði þ 1=2Þ : ð39Þ n :¼ 1 þ tan 1 n þ 1 i¼0 2ðn þ 1Þ best Powell’s technique is quite elegant and needs no assumption on f being differentiable as ours does. Powell also showed that n grows rather slowly as n grows. In fact, 2:414 n 4:037
for 1 n 25:
Analogously to (30), (39) yields 1 cheb n ðfÞ: ðfÞ best
ð40Þ
684
IEEE TRANSACTIONS ON COMPUTERS,
Since 1:27 log2 ð n Þ 2:02 for 1 n 25, (40) is translated into differences of no more than 2 to 3 bits in accuracy if ðfÞ is nearly 1; see Section 6 for more analysis. Similar conclusions on approximations to f0 ðxÞ can be derived, but we omit the detail.
6
AN ASYMPTOTIC ANALYSIS
We shall now perform an asymptotic analysis to explain the implications of our main theorems in the context of elementary function computations, which is the original motivation for our work. In this application, usually
VOL. 53,
NO. 6,
JUNE 2004
2 f 0 ð0Þ < < 2 f 0 ð0Þ 2 f ðnþ2Þ ð0Þ ; þ ln 2 fð0Þ ln 2 fð0Þ ln 2 f ðnþ1Þ ð0Þ
ðmþ1Þ ðmþ1Þ f f 2 ð0Þ < 2 ð0Þ < 0 ðmÞ ðmÞ ðm þ 1Þ ln 2 f ð0Þ ðm þ 1Þ ln 2 f ð0Þ 2 f ðnþ2Þ ð0Þ þ ; ln 2 f ðnþ1Þ ð0Þ ð49Þ
and, if fð0Þ ¼ 0,
¼ ¼ > 0; < 1 is tiny; e:g:; 25 or smaller: ð41Þ This is especially true for table-based methods. Notice that
ðmþ1Þ ðmþ1Þ f f 2 ð0Þ < < 2 ð0Þ ~ ðmÞ ðmÞ ðm þ 1Þ ln 2 f ð0Þ ðm þ 1Þ ln 2 f ð0Þ 2 f ðnþ2Þ ð0Þ þ : ln 2 f ðnþ1Þ ð0Þ ð50Þ
log2 jðrel: err:Þj ¼ ðNo: of accurate bits in an approx:Þ: We are interested in knowing how big the differences, :¼ log2 cheb log2 best ;
ð42Þ
0 :¼ log2 cheb0 log2 best0 ;
ð43Þ
~ :¼ log2 ~cheb log2 ~best ;
ð44Þ
can get. As we shall see later, these differences are on the order of OðÞ, in general, and Oð2 Þ for even or odd functions. It follows from Theorem 5.1 and Theorem 5.3 that log2 ðfÞ log2 ðfÞ þ log2 ðf ðnþ1Þ Þ;
ð45Þ
log2 ðf0 Þ 0 log2 ðf0 Þ þ log2 ðf ðnþ1Þ Þ;
ð46Þ
In elementary function computations, usually these upper bounds are less than 1 and lower bounds bigger than 1 and, thus, the differences in the numbers of correct bits are under 1 bit! The bounds in (48)-(50) are of order OðÞ. This can be improved to order Oð2 Þ for odd (or even) function f. To see this, it suffices to show, under appropriate conditions, the s in both Theorems 5.2 and 5.4 are in a neighborhood of width Oð2 Þ of 1. For the purpose of showing this idea, we take Theorem 5.4 as an example. Under the conditions of Theorem 5.4, we have f0 ðxÞ ¼
f ðmÞ ð0Þ f ðmþ2Þ ð0Þ þ x2 þ Oðx4 Þ; m! ðm þ 2Þ!
f ðnþ2Þ ðxÞ ¼ f ðnþ2Þ ð0Þ þ x2
and, if fð0Þ ¼ 0, log2 ðf0 Þ ~ log2 ðf0 Þ þ log2 ðf ðnþ1Þ Þ:
ð47Þ
What do these inequalities tell us? Suppose fðxÞ is wellbehaved in ½; , as are most elementary mathematical functions. Since is very small, we have4 fðxÞ ¼ fð0Þ þ xf 0 ð0Þ þ Oðx2 Þ; f0 ðxÞ ¼
f ðmÞ ð0Þ f ðmþ1Þ ð0Þ þx þ Oðx2 Þ; m! ðm þ 1Þ!
f ðnþ1Þ ðxÞ ¼ f ðnþ1Þ ð0Þ þ xf ðnþ2Þ ð0Þ þ Oðx2 Þ; which imply 0 f ð0Þ þ Oð2 Þ; ðfÞ 1 þ 2 fð0Þ 2 f ðmþ1Þ ð0Þ þ Oð2 Þ; ðf0 Þ 1 þ m þ 1 f ðmÞ ð0Þ ðnþ2Þ f ð0Þ ðnþ1Þ ðf Þ 1 þ 2 ðnþ1Þ þ Oð2 Þ: ð0Þ f Therefore, ignoring the terms of order Oð2 Þ or higher, we have that 4. For f0 ðxÞ, we notice that m is the algebraic multiplicity of x ¼ 0 as a zero of fðxÞ fð0Þ.
ð48Þ
f ðnþ4Þ ð0Þ þ Oðx4 Þ; 2!
since, for odd (even) functions, even (odd) derivatives at x ¼ 0 are zero. Therefore, ðmþ2Þ f 22 ð0Þ þ Oð4 Þ; ðf0 Þ 1 þ ðmÞ ðm þ 1Þðm þ 2Þ f ð0Þ ðnþ4Þ f ð0Þ ðf ðnþ2Þ Þ 1 þ 2 ðnþ2Þ þ Oð4 Þ: ð0Þ f
7
APPLICATIONS TO ELEMENTARY FUNCTION COMPUTATIONS
In this section, we shall see the implications of the bounds we developed in the previous sections on the most common elementary functions—exponentials, logarithms, and trigonometric functions. Assume (41).
7.1 Functions ax , ax 1 for a > 1 Here, fðxÞ ¼ ax or ax 1, where a > 1, and, then, f0 ðxÞ ¼ ðax 1Þ=x. This includes all natural exponential functions as special cases. We have f 0 ðxÞ ¼ ax ln a, f ðnþ1Þ ðxÞ ¼ ax ðln aÞnþ1 , and f00 ðxÞ ¼
xax ln a ðax 1Þ : x2
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
685
TABLE 1 Bounds on Differences in Numbers of Accurate Bits to ax
It can be proven that f0 ðxÞ > 0 and f00 ðxÞ > 0, using their Maclaurin series. Therefore, for fðxÞ ¼ ax ðfÞ ¼ a2 ¼ ðf ðnþ1Þ Þ;
ðf0 Þ ¼ f0 ðÞ=f0 ðÞ ¼ a :
We have, by (45) and (46) 2 log2 a 4 log2 a;
ð51Þ
log2 a 0 3 log2 a;
ð52Þ
and, if fðxÞ ¼ ax 1 for which f0 ðxÞ is the same as above, polynomials b~n ðxÞ or c~n ðxÞ should be used, by (47) log2 a ~ 3 log2 a:
ð53Þ
Table 1 displays 4 log2 a as an upper bound on the differences in accurate bits by the two approximations as a and ¼ 2N vary.
7.2 Function loga ð1 þ xÞ for a > 1 Here, fðxÞ ¼ loga ð1 þ xÞ for a > 1. This includes the natural logarithm as a special case. Since fð0Þ ¼ 0, f0 ðxÞ ¼ loga ð1 þ xÞ=x and approximations b~n ðxÞ or c~n ðxÞ should be used. We have f 0 ðxÞ ¼ ð1 þ xÞ1 = ln a, and n! ð1 þ xÞn1 ; f ðnþ1Þ ðxÞ ¼ ð1Þ ln a x=ð1 þ xÞ lnð1 þ xÞ : f00 ðxÞ ¼ x2 ln a n
It can be seen that5 for jxj < 1, f0 ðxÞ > 0 and f00 ðxÞ < 0, using their Maclaurin series. Therefore, f0 ðÞ log ð1 Þ lnð1 Þ ; ¼ a f0 ðÞ loga ð1 þ Þ lnð1 þ Þ 1 þ nþ1 ðf ðnþ1Þ Þ ¼ : 1 ðf0 Þ ¼
We have bounds lnð1 Þ lnð1 Þ ~ log2 log2 lnð1 þ Þ lnð1 þ Þ 1þ : þ ðn þ 1Þ log2 1
7.3 Function sin x Here, fðxÞ ¼ sin x. Assume jxj < =4 or smaller and n is odd. Since fð0Þ ¼ 0, f0 ðxÞ ¼ sin x=x and approximations b~n ðxÞ or c~n ðxÞ should be used. We have f 0 ðxÞ ¼ cos x, f ð2kÞ ðxÞ ¼ ð1Þk sin x, f ð2kþ1Þ ðxÞ ¼ ð1Þk cos x, and f00 ðxÞ ¼
x cos x sin x : x2
It can be seen that for f0 ðxÞ > 0, and that f00 ðxÞ > 0 for x < 0 and f00 ðxÞ < 0 for 0 < x, using their Maclaurin series. Therefore, ðf0 Þ ¼
f0 ð0Þ ; ¼ f0 ðÞ sin
ðf ðnþ2Þ Þ ¼
1 : cos
By (47), we have bounds ~ log2 log2 ðcos Þ: log2 sin sin
ð55Þ
The bounds in (55), as well as those in (59) and (60) for cos x, are independent of n. Table 3 displays largest bounds log2 ðcos Þ; 1 :¼ log2 ð56Þ sin 2 :¼ 2 log2 cos ; 3 :¼ log2
2 log2 cos 2ð1 cos Þ
ð57Þ ð58Þ
in (55), and (59) and (60) below as upper bounds on the difference in accurate bits by the two approximations as ¼ 2N varies for both sin x and cos x.
ð54Þ
Table 2 displays the last quantity in (54) as an upper bound on the differences in accurate bits by the two approximations as n and ¼ 2N vary. More than 1 bit differences may occur only for ¼ 25 with n ¼ 10 or higher. Usually, is made much smaller than this. Current HP-UX implementation [12] and 5. They are actually true for any x > 1.
Intel’s implementation [1] for loga for a ¼ 2; e; 10 on Itanium processors [16] uses 28 .
7.4 Function cos x Now, fðxÞ ¼ cos x and then f0 ðxÞ ¼ ðcos x 1Þ=x2 . Assume jxj < =4 or smaller, and n is even. We have f 0 ðxÞ ¼ sin x, f ð2k1Þ ðxÞ ¼ ð1Þk sin x, f ð2kÞ ðxÞ ¼ ð1Þk cos x, and f00 ðxÞ ¼
x sin x 2 cos x þ 2 : x3
It can be proven that f0 ðxÞ < 0 and xf00 ðxÞ 0; so, f0 ðxÞ is an increasing function for x > 0 and a decreasing one for x < 0. Therefore,
686
IEEE TRANSACTIONS ON COMPUTERS,
VOL. 53,
NO. 6,
JUNE 2004
TABLE 2 Bounds on Differences in Numbers of Accurate Bits to loga x
TABLE 3 Bounds on Differences in Numbers of Accurate Bits to sin x and cos x
ðfÞ ¼
1 ¼ ðf ðnþ2Þ Þ; cos
ðf0 Þ ¼
f0 ð0Þ 2 ¼ : f0 ðÞ 2ð1 cos Þ
min jf ðnþ2Þ ðxÞj ¼ f ðnþ2Þ ð0Þ jxj
¼ ðn þ 2Þ!
We have, by (45) and (46), log2 cos 2 log2 cos ; log2
ð59Þ
2 2 0 log2 log2 cos : ð60Þ 2ð1 cos Þ 2ð1 cos Þ
Table 3 displays the upper bounds in (59) and (60) as 2 and 3 .
7.5 Function arcsin x Here, fðxÞ ¼ arcsin x for jxj 1=2 and, then, f0 ðxÞ ¼ arcsin x=x. Assume n is odd. The Maclaurin series of fðxÞ is [17, p. 81]
n!! ¼ ðn!!Þ2 ; ðn þ 1Þ!! ðn þ 2Þ
minjxj f0 ðxÞ ¼ f0 ð0Þ ¼ 1, and maxjxj f0 ðxÞ ¼ f0 ðÞ. Therefore, ðf ðnþ2Þ Þ ¼
1 2
ðn!!Þ
f ðnþ2Þ ðÞ;
ðf0 Þ ¼
f0 ðÞ arcsin ; ¼ f0 ð0Þ
which imply log2
arcsin arcsin ~ log2 1 f ðnþ2Þ ðÞ: þ log2 ðn!!Þ2
ð61Þ
Table 4 displays the last quantity in (61) as an upper bound 3
arcsin x ¼ x þ
5
2kþ1
1x 1 3x ð2k 1Þ!! x þ þ þ þ : 2 3 24 5 ð2kÞ!! 2k þ 1
Thus, maxjxj jf ðnþ2Þ ðxÞj ¼ f ðnþ2Þ ðÞ and
on the differences in accurate bits by the two approximations as n and ¼ 2N vary.
LI: NEAR OPTIMALITY OF CHEBYSHEV INTERPOLATION FOR ELEMENTARY FUNCTION COMPUTATIONS
687
TABLE 4 Bounds on Differences in Numbers of Accurate Bits to arcsin x
8
CONCLUSIONS
[7]
This paper shows that the Chebyshev interpolating polynomials are just as good as the best polynomials in the context of table-driven methods for elementary function computations and, thus, the common practice of feeding Chebyshev interpolating polynomials into a Remez procedure offers negligible improvement, if any, and thus may be considered unnecessary. Our results improve the earlier results due to Powell [6] which, however, have more general applications.
[11]
ACKNOWLEDGMENTS
[12]
The author is grateful to Prof. W. Kahan of the University of California at Berkeley who brought Powell [6] to his attention and to Professor M.J.D. Powell of the University of Cambridge for sending the paper to him. He thanks the referees for their careful reading of the manuscript and constructive suggestions that improved the presentation. Part of this work was done while the author was on leave at Hewlett-Packard Company. He is grateful for help received from Jim Thomas, Jon Okada, and Peter Markstein of the HP Itanium floating-point and elementary math library team at Cupertino, California. This work was supported in part by the US National Science Foundation CAREER award under Grant No. CCR-9875201.
REFERENCES [1] [2] [3] [4] [5] [6]
J. Harrison, T. Kubaska, S. Story, and P.T.P. Tang, “The Computation of Transcendental Functions on the IA-64 Architecture,” Intel Technology J., no. Q4, pp. 1-7, Nov. 1999. R.-C. Li, S. Boldo, and M. Daumas, “Theorems on Efficient Argument Reductions,” Proc. 16th IEEE Symp. Computer Arithmetic, June 2003. P. Markstein, IA-64 and Elementary Functions: Speed and Precision, 2000. J.-M. Muller, Elementary Functions: Algorithms and Implementation. Boston, Basel, Berlin: Birkha˚user, 1997. K.C. Ng, “Argument Reduction for Huge Arguments: Good to the Last Bit,” SunPro, technical report, 1992, http://www.validlab. com/. M.J.D. Powell, “On the Maximum Errors of Polynomial Approximations Defined by Interpolation and by Least Squares Criteria,” The Computer J., vol. 9, no. 4, pp. 404-407, Feb. 1967.
[8] [9]
[10]
[13] [14] [15] [16]
[17]
P.T.P. Tang, “Table-Driven Implementation of the Exponential Function in IEEE Floating-Point Arithmetic,” ACM Trans. Math. Software, vol. 15, no. 2, pp. 144-157, June 1989. P.T.P. Tang, “Table-Driven Implementation of the Logarithm Function in IEEE Floating-Point Arithmetic,” ACM Trans. Math. Software, vol. 16, no. 4, pp. 378-400, Dec. 1990. P.T.P. Tang, “Table Lookup Algorithms for Elementary Functions and Their Error Analysis,” Proc. 10th IEEE Symp. Computer Arithmetic, P. Kornerup and D.W. Matula, eds., pp. 232-236, June 1991. P.T.P. Tang, “Table-Driven Implementation of the expm1 Function in IEEE Floating-Point Arithmetic,” ACM Trans. Math. Software, vol. 18, no. 2, pp. 211-222, June 1992. E. Remez, “Sur un Proce´de´ Convergent d’Approximations Successives pour De´terminer les Polynoˆmes d’Approximation,” C.R. Acade´mie des Sciences, Paris, vol. 198, 1934. R.-C. Li, P. Markstein, J. Okada, and J. Thomas, “The libm Library and Floating-Point Arithmetic in HP-UX for Itanium II,” June 2000, http://h21007.www2.hp.com/dspp/files/unprotected/ Itanium/FP_White_Paper_v2.pdf. E.W. Cheney, Approximation Theory, second ed. Providence, R.I.: AMS, 1981. R.A. DeVore and G.G. Lorentz, Constructive Approximation. Berlin, Heidelberg, New York: Springer-Verlag, 1993. G.G. Lorentz, Approximation of Functions, second ed. New York: Chelsea Publishing, 1986. INTEL, Intel IA-64 Architecture Software Developer’s Manual. Intel Corp., 2002, vol. 3: Instruction Set Reference, document No. 245319-002, http://developer.intel.com/design/itanium/ manuals/index.htm. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, M. Abramowitz and I.A Stegun, eds., ninth ed. New York: Dover Publications, 1970.
Ren-Cang Li received the BS degree in computational mathematics from Xiamen University, People’s Republic of China, in 1985, the MS degree, also in computational mathematics, from the Chinese Academy of Science in 1988, and the PhD degree in applied mathematics from the University of California at Berkeley in 1995. He is currently an associate professor in the Department of Mathematics, University of Kentucky, Lexington. He was awarded the 1995 Householder Fellowship in Scientific Computing by Oak Ridge National Laboratory, a Friedman memorial prize in Applied Mathematics from the University of California at Berkeley in 1996, and a CAREER award from the US National Science Foundation in 1999. His research interest includes elementary function computations, floating-point support for scientific computing, large and sparse linear systems, eigenvalue problems, and model reduction, and unconventional schemes for differential equations.
