Accurate Evaluation of Polynomials in Legendre Basis

Abstract

This paper presents a compensated algorithm for accurate evaluation of a polynomial in Legendre basis. Since the coefficients of the evaluated polynomial are fractions, we propose to store these coefficients in two floating point numbers, such as double-double format, to reduce the effect of the coefficients’ perturbation. The proposed algorithm is obtained by applying error-free transformation to improve the Clenshaw algorithm. It can yield a full working precision accuracy for the ill-conditioned polynomial evaluation. Forward error analysis and numerical experiments illustrate the accuracy and efficiency of the algorithm.

1. Introduction

Legendre polynomial is often used in numerical analysis [1–3], such as approximation theory and quadrature and differential equations. Legendre polynomial satisfies 3-term recurrence relation; that is, for Legendre polynomial 𝑝_𝑘(𝑥), 𝑝₀(𝑥)𝑝=1₁(𝑥)𝑝=𝑥_𝑘+1(𝑥)=2𝑘+1𝑘+1𝑥𝑝_𝑘(𝑥)−𝑘𝑝𝑘+1^𝑘−1(𝑥)(𝑘>0). (1) The polynomial represented in Legendre basis is 𝑝(𝑥)=∑^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥), where 𝑎_𝑗∈ℝ and 𝑝_𝑗(𝑥) is Legendre polynomial.

The Clenshaw algorithm [4, 5] is usually used to evaluate a linear combination of Chebyshev polynomials, but it can apply to any class of functions that can be defined by a three-term recurrence relation. Therefore the Clenshaw algorithm can evaluate a polynomial in Legendre basis. The error analysis of the Clenshaw algorithm was considered in the literatures [6–10]. The relative accuracy bound of the computed values 𝑝(𝑥) by the Clenshaw algorithm verifies 𝑝(𝑥)−𝑝(𝑥)𝑝(𝑥)≤cond𝑝,𝑥×𝑜(𝑢). (2)

For ill-conditioned problems, several researches applied error-free transformations [11] to propose accurate compensated algorithms [12–15] to evaluate the polynomials in monomial, Bernstein, and Chebyshev bases with Horner, de Casteljau, and Clenshaw algorithms, respectively. Some recent applications of high-precision arithmetic were given in [16].

Motivated by them, we apply error-free transformations to analyze the effect of round-off errors and then compensate them to the original result of the Clenshaw algorithm. Since the coefficients of the Legendre polynomial are fractions, the coefficient perturbations in the evaluation may exist when the coefficients are truncated to floating point numbers. We store the coefficients which are not floating point numbers in double-double format, with the double working precision, to get the perturbation. We also compensate the approximate perturbed errors to the original result of the Clenshaw algorithm. Based on the above, we construct a compensated Clenshaw algorithm for the evaluation of a linear combination of Legendre polynomials, which can yield a full working precision accuracy and its relative accuracy bound satisfies 𝑝(𝑥)−𝑝(𝑥)𝑝(𝑥)≤𝑜(𝑢)+cond𝑝,𝑥×𝑜𝑢². (3)

The paper is organized as follows. Section 2 shows some basic notations in error analysis, floating point arithmetic, error-free transformations, compensated algorithm, Clenshaw algorithm, and condition number. Section 3 presents the compensated algorithm and its error bound. Section 4 gives several numerical experiments to illustrate the efficiency and accuracy of the compensated algorithm for polynomial in Legendre basis.

2. Mathematical and Arithmetical Preliminaries

2.1. Basic Notations and Definitions

Throughout this paper, we assume to work with a floating point arithmetic adhering to IEEE-754 floating point standard in rounding to nearest and no overflow nor underflow occurs. Let 𝑜𝑝∈{⊕,⊖,⊗,⊘} represent the floating point computation; then the computation obeys the model 𝑎𝑜𝑝𝑏=(𝑎∘𝑏)1+𝜀₁=(𝑎∘𝑏)1+𝜀₂, (4) where ∘∈{+,−,×,÷} and |𝜀₁|,|𝜀₂|≤𝑢 (𝑢 is the round-off unit). We also assume that the computed result of 𝑎∈ℝ in floating point arithmetic is denoted by 𝑎 and the set of all floating point numbers is denoted by 𝔽. The following definition will be used in error analysis (see more details in [17]).

Definition 1. One defines ⟨𝑛⟩∶=^𝑛_𝑖=11+𝛿_𝑖^𝜌𝑖=1+𝜃_𝑛, (5) where |𝛿_𝑖|≤𝑢,𝜌_𝑖=±1 with ∀𝑖=1,2,…,𝑛, |𝜃_𝑛|≤𝛾_𝑛∶=(𝑛𝑢/(1−𝑛𝑢))=𝑛𝑢+𝑂(𝑢²), and 𝑛𝑢<1.

There are three classic properties which will also be used in error analysis:(i)⟨𝑘⟩⟨𝑗⟩=⟨𝑘+𝑗⟩;(ii)𝛾_𝑘<𝛾_𝑘+1;(iii)𝛾_𝑘+𝛾_𝑗+𝛾_𝑘𝛾_𝑗≤𝛾_𝑘+𝑗.

2.2. Accurate Sum and Product

Let 𝑎,𝑏∈𝔽, and no overflow nor underflow occurs. The transformation (𝑎,𝑏)→(𝑥,𝑦) is regarded as an error-free transformation (EFT) that causes 𝑦∈𝔽 to exist such that 𝑎∘𝑏=𝑥+𝑦, 𝑥=𝑎𝑜𝑝𝑏.

Let us show the error-free transformations of the sum and product of two floating point numbers in Algorithms 1-3 which are the 𝑇𝑤𝑜𝑆𝑢𝑚 algorithm by Knuth [18] and the 𝑇𝑤𝑜𝑃𝑟𝑜𝑑 algorithm by Dekker [19], respectively.

function [𝑥,𝑦]=𝑇𝑤𝑜𝑆𝑢𝑚(𝑎,𝑏) 𝑥=𝑎⊕𝑏 𝑧=𝑥⊖𝑎 𝑦=(𝑎⊖(𝑥⊖𝑧))⊕(𝑏⊖𝑧)

Algorithm 1 

Sum of two floating point numbers.

function [𝑥,𝑦]=𝑆𝑝𝑙𝑖𝑡(𝑎) 𝑐=𝑓𝑎𝑐𝑡𝑜𝑟⊗𝑎 (in double precision factor =227+1) 𝑥=𝑐⊖(𝑐⊖𝑎) 𝑦=𝑎⊖𝑥

Algorithm 2 

Split of a floating point number into two parts.

function [𝑥,𝑦]=𝑇𝑤𝑜𝑃𝑟𝑜𝑑(𝑎,𝑏) 𝑥=𝑎⊗𝑏 [𝑎ℎ,𝑎𝑙]=𝑆𝑝𝑙𝑖𝑡(𝑎) [𝑏ℎ,𝑏𝑙]=𝑆𝑝𝑙𝑖𝑡(𝑏) 𝑦=𝑎𝑙⊗𝑏𝑙⊖(((𝑥⊖𝑎ℎ⊗𝑏ℎ)⊖𝑎𝑙⊗𝑏ℎ)⊖𝑎ℎ⊗𝑏𝑙)

Algorithm 3 

Product of two floating point numbers.

Algorithms 1–3 satisfy the Theorem 2.

Theorem 2 (see [11]). For 𝑎,𝑏∈𝔽 and 𝑥,𝑦∈𝔽, 𝑇𝑤𝑜𝑆𝑢𝑚 and 𝑇𝑤𝑜𝑃𝑟𝑜𝑑 verify ()𝑦{𝑢|𝑥||,()𝑦{𝑢|𝑥||.𝑥,𝑦=𝑇𝑤𝑜𝑆𝑢𝑚𝑎,𝑏,𝑥+𝑦=𝑎+𝑏,≤min,𝑢𝑎+𝑏|}𝑥,𝑦=𝑇𝑤𝑜𝑃𝑟𝑜𝑑𝑎,𝑏,𝑥+𝑦=𝑎𝑏,≤min,𝑢𝑎𝑏|} (6)

We present the compensated algorithm for the product of three floating point numbers in Algorithm 4, which refers to [20].

function [𝑥,𝑦]=𝑇ℎ𝑟𝑒𝑒𝑃𝑟𝑜𝑑(𝑎,𝑏,𝑐) [𝑡,ℎ]=𝑇𝑤𝑜𝑃𝑟𝑜𝑑(𝑎,𝑏) [𝑥,𝑒]=𝑇𝑤𝑜𝑃𝑟𝑜𝑑(𝑡,𝑐) 𝑦=𝑐⊗ℎ⊕𝑒

Algorithm 4 

Product of three floating point numbers.

According to Theorem 2, we have |ℎ|≤𝑢|𝑎𝑏| and 𝑡+ℎ=𝑎𝑏. Hence, |𝑒|≤𝑢|(𝑎𝑏−ℎ)𝑐|≤𝑢|𝑎𝑏𝑐|+𝑢|ℎ𝑐|≤𝑢|𝑎𝑏𝑐|+𝑢²|𝑎𝑏𝑐|; then |𝑐ℎ+𝑒|≤|𝑐ℎ|+|𝑒|≤2𝑢|𝑎𝑏𝑐|+𝑢²|𝑎𝑏𝑐|. According to Lemma 3.2 in [20], we can propose the error bound of Algorithm 4 in Theorem 3.

Theorem 3. For 𝑎,𝑏,𝑐∈𝔽 and 𝑥,𝑦,𝜏∈𝔽, 𝑇ℎ𝑟𝑒𝑒𝑃𝑟𝑜𝑑 verifies ()|𝜏|𝑥,𝑦=𝑇ℎ𝑟𝑒𝑒𝑃𝑟𝑜𝑑𝑎,𝑏,𝑐,𝑥+𝑦=𝑎𝑏𝑐−𝜏,≤𝛾₂𝛾₆|𝑎𝑏𝑐|,𝑦≤2𝑢|𝑎𝑏𝑐|+14𝑢²|𝑎𝑏𝑐|. (7)

Proof. According to Lemma 3.2 in [20], we have |𝜏|≤𝛾₂𝛾₆|𝑎𝑏𝑐|. (8) In Algorithm 4, 𝑦+𝜏=𝑐ℎ+𝑒; then 𝑦=≤||+|𝜏|||𝑦+𝜏−𝜏𝑐ℎ+𝑒≤2𝑢𝑎𝑏𝑐+𝑢²|𝑎𝑏𝑐|+𝛾₂𝛾₆||||𝑎𝑏𝑐≤2𝑢𝑎𝑏𝑐+14𝑢²|𝑎𝑏𝑐|. (9)

2.3. The Clenshaw Algorithm and Condition Number

The standard general algorithm for the evaluation of polynomial in Legendre basis 𝑝(𝑥)=∑^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥) is the Clenshaw algorithm [5]. We recall it in Algorithm 5.

function 𝑝(𝑥)=𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤(𝑝,𝑥) 𝑏𝑛+2=𝑏𝑛+1=0 for 𝑗=𝑛∶−1∶0  𝑏𝑗=2𝑗+1𝑗+1𝑥𝑏𝑗+1−𝑗+1𝑏𝑗+2𝑗+2+𝑎𝑗 end 𝑝(𝑥)=𝑏0

Algorithm 5 

The Clenshaw algorithm for finite Legendre series.

Barrio et al. [21] proposed a general polynomial condition number for any polynomial basis defined by a linear recurrence and used this new condition number to give the error bounds for the Clenshaw algorithm. Based on this general condition number, we propose the absolute Legendre polynomial, which is similar to the absolute polynomial mentioned in [7, 21].

Definition 4. Let 𝑝_𝑘(𝑥) be Legendre polynomial. We define the absolute Legendre polynomial 𝑝_𝑘(𝑥) which is associated with 𝑝_𝑘(𝑥) and satisfies 𝑝₀(𝑥)𝑝=1₁(𝑥)𝑝=𝑥_𝑘(𝑥)=2𝑘−1𝑘𝑥𝑝_𝑘−1(𝑥)+𝑘𝑝𝑘+1^𝑘−2(𝑥)(𝑘≥2), (10) where |𝑝_𝑘(𝑥)|≤𝑝_𝑘(𝑥), ∀𝑥≥0.

The absolute Legendre polynomial satisfies the following property.

Lemma 5. Let 𝑝_𝑘(𝑥) be the absolute Legendre polynomial. Then we have 𝑝_𝑘(𝑥)𝑝_𝑗(𝑥)≤𝑝_𝑘+𝑗(𝑥). (11)

Proof . Let 𝐴_𝑘=(2𝑘+1)/(𝑘+1) and 𝑘≥2. From Definition 4, we have 𝑝_𝑘(𝑥)=𝐴_𝑘−1𝑥𝑝_𝑘−1(𝑥)+𝑘𝑝𝑘+1^𝑘−2(𝑥)=^𝑘−1_𝑖=0𝐴_𝑖𝑥^𝑘+⋯; (12) then 𝑝_𝑘(𝑥)𝑝_𝑗(𝑥)=^𝑘−1_𝑖=0𝐴_𝑖^𝑗−1_𝑖=0𝐴_𝑖𝑥^𝑘+𝑗+⋯. (13) Thus the equivalent form of (11) is ^𝑘−1_𝑖=0𝐴_𝑖^𝑗−1_𝑖=0𝐴_𝑖≤^{𝑘+𝑗−1}_𝑖=0𝐴_𝑖. (14) Since 𝐴_𝑘 is increasing with 𝑘, we obtain ^𝑗−1_𝑖=0𝐴_𝑖≤^{𝑘+𝑗−1}_𝑖=𝑘𝐴_𝑖. (15) So we finally obtain 𝑝_𝑘(𝑥)𝑝_𝑗(𝑥)≤𝑝_𝑘+𝑗(𝑥).

Following Definition 4, we introduce the condition number for the evaluation of polynomials in Legendre basis [21].

Definition 6. Let 𝑝(𝑥)=∑^𝑛_𝑘=0𝑎_𝑘𝑝_𝑘(𝑥), where 𝑝_𝑘(𝑥) is Legendre polynomial. Let 𝑝_𝑘(𝑥) be the absolute Legendre polynomial. Then the absolute condition number is cond_abs𝑝,𝑥=𝑝(|𝑥|)=^𝑛_𝑘=0𝑎_𝑘𝑝_𝑘(|𝑥|), (16) and the relative condition number is cond_rel=𝑝𝑥𝑝,𝑥(||)𝑝(𝑥=∑)^𝑛𝑘=0𝑎_𝑘𝑝_𝑘(|𝑥|)𝑝(𝑥.) (17)

3. Compensated Algorithm for Evaluating Polynomials

In this section, we exhibit the exact round-off errors generated by the Clenshaw algorithm with EFT. We also analyze the perturbations generated by truncating the fractions in Algorithm 5. We propose a compensated Clenshaw algorithm to evaluate finite Legendre series and present its error bound in the following.

Firstly, in order to analyze the perturbations, we split each coefficient into three parts as follows: 𝐴=𝐴^(ℎ)+𝐴^(𝑙)+𝐴^(𝑚), (18) where 𝐴^(ℎ),𝐴^(𝑙)∈𝔽, 𝐴,𝐴^(𝑚)∈ℝ and |𝐴^(𝑙)|≤𝑢|𝐴^(ℎ)|,|𝐴^(𝑚)|≤𝑢|𝐴^(𝑙)|. Let 𝐴_𝑘=(2𝑘+1)/(𝑘+1) and 𝐶_𝑘=−(𝑘/(𝑘+1)). Then we describe the recurrence relation at 𝑗th step of Algorithm 5 for theoretical computation as 𝑏_𝑗=𝐴_𝑗^(ℎ)+𝐴_𝑗^(𝑙)+𝐴_𝑗^(𝑚)𝑥𝑏_𝑗+1+𝐶^(ℎ)_𝑗+1+𝐶^(𝑙)_𝑗+1+𝐶^(𝑚)_𝑗+1𝑏_𝑗+2+𝑎_𝑗^(ℎ)+𝑎_𝑗^(𝑙)+𝑎_𝑗^(𝑚). (19) For numerical computation associated with (19), it is 𝑏̂_𝑗=𝐴_𝑗^(ℎ)⊗𝑥⊗𝑏̂_𝑗+1⊕𝐶^(ℎ)_𝑗+1⊗𝑏̂_𝑗+2⊕𝑎_𝑗^(ℎ). (20)

Remark 7. Let 𝐴∈ℝ. We split a coefficient 𝐴 like (18); then the representation of splitting is unique and 𝐴=𝐴^(ℎ)⊕𝐴^(𝑙), 𝐴=𝐴⟨1⟩.

Let 𝐴_𝑗=𝐴_𝑗^(ℎ) and 𝐶_𝑗=𝐶_𝑗^(ℎ) in Algorithm 5. Since every elementary floating point operation in Algorithm 5 causes round-off errors in numerical computation, we apply 𝐸𝐹𝑇 and the 𝑇ℎ𝑟𝑒𝑒𝑃𝑟𝑜𝑑 algorithms to take notes of all round-off errors and obtain 𝑟,𝛼_𝑗=𝑇ℎ𝑟𝑒𝑒𝑃𝑟𝑜𝑑𝐴_𝑗^(ℎ),𝑥,𝑏̂_𝑗+1,𝑠,𝛽_𝑗=𝑇𝑤𝑜𝑃𝑟𝑜𝑑𝐶^(ℎ)_𝑗+1,𝑏̂_𝑗+2,𝑡,𝜂_𝑗(),𝑏̂=𝑇𝑤𝑜𝑆𝑢𝑚𝑟,𝑠_𝑗,𝜉_𝑗=𝑇𝑤𝑜𝑆𝑢𝑚𝑡,𝑎_𝑗^(ℎ). (21) The sum of the perturbation and the round-off errors of the recurrence relation at 𝑗th step is 𝑤_𝑗=𝐴_𝑗^(𝑙)+𝐴_𝑗^(𝑚)𝑥𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1+𝐶^(𝑚)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙)+𝑎_𝑗^(𝑚)+𝛼_𝑗+𝜏_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗, (22) where 𝑟+𝛼_𝑗+𝜏_𝑗=𝐴_𝑗^(ℎ)𝑥𝑏̂_𝑗+1, 𝜏_𝑗 is defined in Theorem 3. Then we obtain the following theorem.

Theorem 8. Let 𝑝(𝑥)=∑^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥) be a Legendre series of degree 𝑛, and let 𝑥 be a floating point value. One assumes that 𝐴_𝑗=𝐴_𝑗^(ℎ) and 𝐶_𝑗=𝐶_𝑗^(ℎ) in Algorithm 5 and 𝑏̂₀ is the numerical result; 𝑤_𝑗 is described in (22) for 𝑗=0,1,…,𝑛−1. Then one obtains ^𝑛−1_𝑗=0𝑤_𝑗𝑝_𝑗(𝑥)+𝑏̂₀=^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥). (23)

Proof. The perturbation of Algorithm 5 is 𝐴_𝑗^(𝑙)+𝐴_𝑗^(𝑚)𝑥𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1+𝐶^(𝑚)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙)+𝑎_𝑗^(𝑚). (24) Let 𝐴_𝑗=𝐴_𝑗^(ℎ) and 𝐶_𝑗=𝐶_𝑗^(ℎ) in Algorithm 5. From Theorems 2–3 and (21), we can easily obtain 𝐴_𝑗^(ℎ)𝑏̂_𝑗+1𝑥=𝑟+𝛼_𝑗+𝜏_𝑗,𝐶^(ℎ)_𝑗+1𝑏̂_𝑗+2=𝑠+𝛽_𝑗,𝑟+𝑠=𝑡+𝜂_𝑗𝑡+𝑎_𝑗^(ℎ)=𝑏̂_𝑗+𝜉_𝑗. (25)
Let 𝑠_𝑗=𝛼_𝑗+𝜏_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗 for 𝑗=0,1,…,𝑛−1; then 𝑏̂_𝑗+𝑠_𝑗=𝐴_𝑗^(ℎ)𝑏̂_𝑗+1𝑥+𝐶^(ℎ)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(ℎ),𝑗=0,1,…,𝑛−1. (26)
From (18), (19), (22), and (24), we have 𝑏̂_𝑗+𝑤_𝑗=𝐴_𝑗𝑏̂_𝑗+1𝑥+𝐶_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗,𝑗=0,1,…,𝑛−1. (27)
Thus 𝑏̂₀=^𝑛_𝑗=0𝑎_𝑗−𝑤_𝑗𝑝_𝑗(𝑥); (28) then ^𝑛−1_𝑗=0𝑤_𝑗𝑝_𝑗(𝑥)+𝑏̂₀=^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥). (29)

To devise a compensated algorithm of Algorithm 5 and give its error bound, we need the following lemmas.

Lemma 9. Let 𝑝(𝑥)=∑^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥) be a Legendre series of degree 𝑛, 𝑥 a floating point value, and 𝑏̂₀ the numerical result of Algorithm 5. Then 𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤𝑝,𝑥−𝑝(𝑥)≤𝛾^𝑛_{5𝑛−1𝑗=0}𝑎_𝑗𝑝_𝑗(|𝑥|). (30)

Proof. Applying the standard model of floating point arithmetic, from Algorithm 5 and Remark 7, we have 𝑏̂_𝑛−1=𝐴_𝑛−1𝑏̂_𝑛𝑥⟨4⟩+𝑎_𝑛−1⟨1⟩,𝑏̂_𝑛−2=𝐴_𝑛−2𝑏̂_𝑛−1𝑥⟨5⟩+𝐶_𝑛−1𝑏̂_𝑛⟨4⟩+𝑎_𝑛−2⟨1⟩,⋮𝑏̂₀=𝐴₀𝑏̂₁𝑥⟨5⟩+𝐶₁𝑏̂₂⟨4⟩+𝑎₀⟨1⟩. (31)
Thus we obtain 𝑏̂₀=⟨4+5(𝑛−1)⟩^𝑛−1_𝑖=0𝐴_𝑖𝑥^𝑛𝑎_𝑛+⋯+𝑎₀⟨1⟩=⟨5𝑛−1⟩^𝑛−1_𝑖=0𝐴_𝑖𝑥^𝑛𝑎_𝑛+⋯+𝑎₀⟨1⟩. (32)
Then, from Definition 1, we get 𝑏̂₀−𝑝(𝑥)=^𝑛−1_𝑖=0𝐴_𝑖𝑥^𝑛𝑎_𝑛𝜃_5𝑛−1+⋯+𝑎₀𝜃₁≤𝛾_5𝑛−1^𝑛−1_𝑖=0𝐴_𝑖|𝑥|^𝑛𝑎_𝑛+⋯+𝑎₀=𝛾^𝑛_{5𝑛−1𝑗=0}𝑎_𝑗𝑃_𝑗(|𝑥|). (33)

Lemma 10. Let 𝑝(𝑥)=∑^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥) be a Legendre series of degree 𝑛 and 𝑥 a floating point value. We assume that 𝜎_𝑗≤𝜔_𝑗(𝐴_𝑗|𝑥||𝑏̂_𝑗+1|+|𝐶_𝑗+1||𝑏̂_𝑗+2|+|𝑎_𝑗|), where 𝜔_𝑗 is real numbers; then ^𝑛−1_𝑗=0𝜎_𝑗𝑝_𝑗(|𝑥|)≤𝑛𝜔_𝑗1+𝛾_5(𝑛−1)^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(|𝑥|). (34)

Proof. According to (31), we get 𝑏̂_𝑗≤1+𝛾₅𝐴_𝑗|𝑥|𝑏̂_𝑗+1+𝐶_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗. (35) Since 𝜎_𝑗≤𝜔_𝑗(𝐴_𝑗|𝑥||𝑏̂_𝑗+1|+|𝐶_𝑗+1||𝑏̂_𝑗+2|+|𝑎_𝑗|), from Definition 1, we have 𝜎_𝑗≤𝜔_𝑗1+𝛾_{5(𝑛−𝑗−1)}Σ^𝑛_𝑘=𝑗𝑎_𝑘𝑝_𝑘−𝑗(|𝑥|); (36) then ^𝑛−1_𝑗=0𝜎_𝑗𝑝_𝑗(|𝑥|)≤𝜔_𝑗1+𝛾_5(𝑛−1)×^𝑛−1_𝑗=0^𝑛_𝑘=𝑗𝑎_𝑘𝑝_𝑘−𝑗(|𝑥|)𝑝_𝑗(|𝑥|). (37) According to Lemma 5, we obtain ^𝑛−1_𝑗=0𝜎_𝑗𝑝_𝑗(|𝑥|)≤𝜔_𝑗1+𝛾_5(𝑛−1)^𝑛−1𝑛_{𝑗=0𝑘=𝑗}𝑎_𝑘𝑝_𝑘(|𝑥|)≤𝑛𝜔_𝑗1+𝛾_5(𝑛−1)^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(|𝑥|). (38)

Among the perturbed and round-off errors in Algorithm 5, we deem that some errors do not influence the numerical result in working precision. Now we can give the perturbed error bounds and the round-off error bounds, respectively. At first we analyze the perturbed error in (24). Let 𝑔_𝑗^(𝑙)=𝐴_𝑗^(𝑙)𝑥𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙). According to 𝐴^(𝑙)≤𝑢𝐴^(ℎ) and Remark 7, we obtain 𝑔_𝑗^(𝑙)≤𝑢1+𝛾₁𝐴_𝑗|𝑥|𝑏̂_𝑗+1+𝐶_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗. (39) Then let 𝜔_𝑗=𝑢(1+𝛾₁) and from Lemma 10, taking into account that (𝑛+1)𝑢(1+𝛾_5𝑛+1)≤(5𝑛+1)𝑢(1+𝛾_5𝑛+1)=𝛾_5𝑛+1, we have ^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)≤𝛾^𝑛_{5𝑛+1𝑗=0}𝑎_𝑗𝑝_𝑗(|𝑥|). (40)

𝛾_5𝑛+1 is 𝑂(𝑢), so this coefficient perturbation may influence the accuracy; we need to consider it in our compensated algorithm.

Similarly, we let 𝑔_𝑗^(𝑚)=𝐴_𝑗^(𝑚)𝑥𝑏̂_𝑗+1+𝐶^(𝑚)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑚); then ^𝑛_𝑗=0𝑔_𝑗^(𝑚)𝑝_𝑗(𝑥)≤𝑢𝛾^𝑛_{5𝑛+1𝑗=0}𝑎_𝑗𝑝_𝑗(|𝑥|). (41)

𝑢𝛾_5𝑛+1 is 𝑂(𝑢²), so this coefficient perturbation does not influence the accuracy.

Remark 11. When 𝑗=𝑛, the 𝑛th step of Algorithm 5 is 𝑏̂_𝑛=𝑎_𝑛^(ℎ), so that we only need to consider the perturbation of coefficient 𝑎_𝑛.

Next we deduce the round-off error bound. Let 𝑠_𝑗=𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗. According to Theorems 2 and 3 and Remark 7, we obtain 𝛼_𝑗≤2𝑢+14𝑢²𝐴_𝑗^(ℎ)𝑥𝑏̂_𝑗+1≤2𝑢+14𝑢²1+𝛾₁𝐴_𝑗𝑥𝑏̂_𝑗+1;𝛽_𝑗≤𝑢𝐶^(ℎ)_𝑗+1𝑏̂_𝑗+2≤𝑢1+𝛾₁𝐶_𝑗+1𝑏̂_𝑗+2;𝜂_𝑗≤𝑢𝐴_𝑗^(ℎ)𝑥𝑏̂_𝑗+1⟨3⟩+𝐶^(ℎ)_𝑗+1𝑏̂_𝑗+2⟨2⟩≤𝑢1+𝛾₃𝐴_𝑗𝑥𝑏̂_𝑗+1+𝐶_𝑗+1𝑏̂_𝑗+2;𝜉_𝑗≤𝑢𝐴_𝑗^(ℎ)𝑥𝑏̂_𝑗+1⟨4⟩+𝐶^(ℎ)_𝑗+1𝑏̂_𝑗+2⟨3⟩+𝑎_𝑗⟨1⟩≤𝑢1+𝛾₄𝐴_𝑗𝑥𝑏̂_𝑗+1+𝐶_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗; (42) then 𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗≤4𝑢+14𝑢²1+𝛾₄×𝐴_𝑗|𝑥|𝑏̂_𝑗+1+𝐶_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗. (43)

Let 𝜔_𝑗=(4𝑢+14𝑢²)|1+𝛾₄| and using it into Lemma 10, from (4+14𝑢)𝑛𝑢(1+𝛾_5𝑛−1)≤(5𝑛−1)𝑢(1+𝛾_5𝑛−1) = 𝛾_5𝑛−1(𝑛≥2), we have ^𝑛−1_𝑗=0𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗𝑝_𝑗(𝑥)≤𝛾^𝑛_{5𝑛−1𝑗=0}𝑎_𝑗𝑝_𝑗(|𝑥|). (44)

𝛾_5𝑛−1 is 𝑂(𝑢), so this round-off error may influence the accuracy, we also need to consider it in our compensated algorithm.

From Theorem 3 we have known that round-off error |𝜏_𝑗|≤𝛾₂𝛾₆|𝐴_𝑗^(ℎ)𝑥𝑏̂_𝑗+1|≤𝛾₂𝛾₆(1+𝛾₁)|𝐴_𝑗𝑥𝑏̂_𝑗+1|. According to Lemma 10 we let 𝜔_𝑗=𝛾₂𝛾₆(1+𝛾₁), from 𝛾₂≤2𝑢(1+𝛾₁) and 2𝑛𝑢(1+𝛾_5𝑛−3)≤𝛾_5𝑛−3, we have ^𝑛−1_𝑗=0𝜏_𝑗𝑝_𝑗(𝑥)≤𝛾₆𝛾^𝑛_{5𝑛−3𝑗=0}𝑎_𝑗𝑝_𝑗(|𝑥|). (45)

𝛾₆𝛾_5𝑛−3 is 𝑂(𝑢²), so this round-off error does not influence the accuracy.

Observing the error bounds we described above, the perturbation generated by the third part of coefficients in (18) does not influence the accuracy. Thanks to the 𝑑𝑖𝑣_𝑑_𝑑 algorithm in Appendix B (see Algorithm 13), we only need to split the coefficients into two floating point numbers.

Applying EFT, the 𝑇ℎ𝑟𝑒𝑒𝑃𝑟𝑜𝑑 algorithm and the 𝑑𝑖𝑣_𝑑_𝑑 algorithm, considering all errors which may influence the numerical result in working precision in Algorithm 5, we obtain the compensated 𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithm in Algorithm 6.

function 𝑟𝑒𝑠=𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤(𝑝,𝑥) 𝑏̂𝑛+2=𝑏̂𝑛+1=0 𝜖𝑏̂𝑛+2=𝜖𝑏̂𝑛+1=0 for 𝑗=𝑛∶−1∶0  [𝐴𝑗(ℎ),𝐴𝑗(𝑙)]=𝑑𝑖𝑣_𝑑_𝑑(2𝑗+1,𝑗+1)  [𝐶𝑗(ℎ),𝐶𝑗(𝑙)]=𝑑𝑖𝑣_𝑑_𝑑(−𝑗−1,𝑗+2)  [𝑟,𝛼𝑗]=𝑇ℎ𝑟𝑒𝑒𝑃𝑟𝑜𝑑(𝐴𝑗(ℎ),𝑏̂𝑗+1,𝑥)  [𝑠,𝛽𝑗]=𝑇𝑤𝑜𝑃𝑟𝑜𝑑(𝐶(ℎ)𝑗+1,𝑏̂𝑗+2)  [𝑡,𝜂𝑗]=𝑇𝑤𝑜𝑆𝑢𝑚(𝑟,𝑠)  [𝑏̂𝑗,𝜉𝑗]=𝑇𝑤𝑜𝑆𝑢𝑚(𝑡,𝑎𝑗)  𝑠𝑗̂=𝛼𝑗⊕𝛽𝑗⊕𝜂𝑗⊕𝜉𝑗⊕𝐴𝑗(𝑙)⊗𝑥⊗𝑏̂𝑗+1⊕𝐶(𝑙)𝑗+1⊗𝑏̂𝑗+2⊕𝑎𝑗(𝑙)  𝜖𝑏̂𝑗=𝐴𝑗(ℎ)⊗𝑥⊗𝜖𝑏̂𝑗+1⊕𝐶(ℎ)𝑗+1⊗𝜖𝑏̂𝑗+2⊕𝑠𝑗̂ end 𝑟𝑒𝑠=𝑦⊕𝜖𝑏̂0

Algorithm 6 

Compensated Clenshaw algorithm.

Here we give the error bound of Algorithm 6.

Theorem 12. Let 𝑝(𝑥)=∑^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥) be a Legendre series of degree 𝑛(𝑛≥2) and 𝑥 a floating point value. The forward error bound of the compensated Clenshaw algorithm is (𝑥𝑝(𝑥𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤𝑝,𝑥−𝑝)≤𝑢)+2𝛾^2𝑛_{5𝑛−2𝑗=0}𝑎_𝑗𝑃_𝑗(𝑥). (46)

Proof. From Algorithm 6, we obtain (𝑥=𝑏̂𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤𝑝,𝑥−𝑝)₀⊕𝜖𝑏̂₀(𝑥=)𝑏̂−𝑝)(1+𝜖₀+𝜖𝑏̂₀−𝑝(𝑥), (47) where |𝜖|≤𝑢.
According to Theorem 8, we have 𝑝(𝑥)=𝑏̂_𝑜+𝜖𝑏_𝑜; then (𝑥=)𝑃(𝑥)𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤𝑝,𝑥−𝑝)(1+𝜖−𝜖𝑏_𝑜+𝜖𝑏̂₀(𝑥|𝑃(𝑥+(−𝑝)≤𝑢)|1+𝑢)𝜖𝑏_𝑜−𝜖𝑏̂₀. (48) Next we analyze the bound of |𝜖𝑏_𝑜−𝜖𝑏̂₀|.
Let 𝑔_𝑗^(𝑙)=𝐴_𝑗^(𝑙)𝑥𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙); then 𝑔_𝑗^(𝑙)=𝐴_𝑗^(𝑙)𝑥𝑏̂_𝑗+1⟨4⟩+𝐶^(𝑙)_𝑗+1𝑏̂_𝑗+2⟨3⟩+𝑎_𝑗^(𝑙)⟨1⟩. (49)
Thus we obtain ^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)≤𝛾₄^𝑛_𝑗=0𝐴_𝑗^(𝑙)|𝑥|𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙)𝑝_𝑗(|𝑥|). (50)
According to Lemma 9, we get ^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−⊕^𝑛_𝑗=0𝑔_𝑗^(𝑙)⊗𝑝_𝑗(𝑥)≤𝛾^𝑛_{5𝑛−1𝑗=0}𝑔_𝑗^(𝑙)𝑝_𝑗(|𝑥|)≤𝛾_5𝑛−11+𝛾₄×^𝑛_𝑗=0𝐴_𝑗^(𝑙)|𝑥|𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙)𝑝_𝑗(|𝑥|). (51)
From (50), (51) and 𝛾_5𝑛−1(1+𝛾₄)+𝛾₄≤𝛾_5𝑛+3, we derive ^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−⊕^𝑛_𝑗=0𝑔_𝑗^(𝑙)⊗𝑝_𝑗(𝑥)≤^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)+^𝑛_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−⊕^𝑛_𝑗=0𝑔_𝑗^(𝑙)⊗𝑝_𝑗(𝑥)≤𝛾^𝑛_{5𝑛+3𝑗=0}𝐴_𝑗^(𝑙)|𝑥|𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙)𝑝_𝑗(|𝑥|). (52)
According to 𝐴^(𝑙)≤𝑢𝐴^(ℎ) and Remark 7, we obtain 𝐴_𝑗^(𝑙)|𝑥|𝑏̂_𝑗+1+𝐶^(𝑙)_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗^(𝑙)≤𝑢1+𝛾₁𝐴_𝑗|𝑥|𝑏̂_𝑗+1+𝐶_𝑗+1𝑏̂_𝑗+2+𝑎_𝑗. (53) From (53) we use 𝜔_𝑗=𝑢(1+𝛾₁) into Lemma 10; taking into account that 𝑛𝑢(1+𝛾_5𝑛−4)≤(5𝑛+1)𝑢(1+𝛾_5𝑛+1)=𝛾_5𝑛+1 and 𝛾_5𝑛+1𝛾_5𝑛+3≤𝛾²_5𝑛+2, according to (52), we get ^𝑛−1_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑔_𝑗^(𝑙)⊗𝑝_𝑗(𝑥)≤𝛾^2𝑛_{5𝑛+2𝑗=0}𝑎_𝑗𝑝_𝑗(|𝑥|). (54)
Next we let 𝑠_𝑗=𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗; then 𝑠_𝑗̂≤1+𝛾₃𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗; (55) thus ^𝑛−1_𝑗=0𝑠_𝑗𝑝_𝑗(𝑥)−^𝑛−1_𝑗=0𝑠_𝑗̂𝑝_𝑗(𝑥)≤𝛾₃^𝑛−1_𝑗=0𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗𝑝_𝑗(|𝑥|). (56)
According to Lemma 9, we obtain 𝑠_𝑗̂^𝑛−1_𝑗=0𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑠_𝑗̂⊗𝑝_𝑗(𝑥)≤𝛾_{5(𝑛−1)−1}^𝑛−1_𝑗=0𝑠_𝑗̂𝑝_𝑗(|𝑥|)≤𝛾_{5(𝑛−1)−1}1+𝛾₃×^𝑛−1_𝑗=0𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗𝑝_𝑗(|𝑥|). (57)
From (56), (57), and 𝛾_{5(𝑛−1)−1}(1+𝛾₃)+𝛾₃≤𝛾_5𝑛−3, we derive ^𝑛−1_𝑗=0𝑠_𝑗𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑠_𝑗̂⊗𝑝_𝑗(𝑥)≤^𝑛−1_𝑗=0𝑠_𝑗𝑝_𝑗(𝑥)−^𝑛−1_𝑗=0𝑠_𝑗̂𝑝_𝑗(𝑥)+^𝑛−1_𝑗=0𝑠_𝑗̂𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑠_𝑗̂⊗𝑝_𝑗(𝑥)≤𝛾_5𝑛−3^𝑛−1_𝑗=0𝛼_𝑗+𝛽_𝑗+𝜂_𝑗+𝜉_𝑗𝑝_𝑗(|𝑥|). (58)
From (43), we let 𝜔_𝑗=(4+14𝑢)𝑢|1+𝛾₄| and using it into Lemma 10, taking into account that (4+14𝑢)𝑛𝑢(1+𝛾_5𝑛−1)≤(5𝑛−1)𝑢(1+𝛾_5𝑛−1)=𝛾_5𝑛−1(𝑛≥2) and 𝛾_5𝑛−3𝛾_5𝑛−1≤𝛾²_5𝑛−2, we have ^𝑛−1_𝑗=0𝑠_𝑗𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑠_𝑗̂⊗𝑝_𝑗(𝑥)≤𝛾²_5𝑛−2^𝑛−1_𝑗=0𝑎_𝑗𝑝_𝑗(|𝑥|). (59)
Combining (54) and (59), we get 𝜖𝑏_𝑜−𝜖𝑏̂₀=^𝑛−1_𝑗=0𝑠_𝑗𝑝_𝑗(𝑥)+^𝑛−1_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑠_𝑗̂⊗𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑔_𝑗^(𝑙)⊗𝑝_𝑗(𝑥)≤^𝑛−1_𝑗=0𝑠_𝑗𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑠_𝑗̂⊗𝑝_𝑗(𝑥)+^𝑛−1_𝑗=0𝑔_𝑗^(𝑙)𝑝_𝑗(𝑥)−⊕^𝑛−1_𝑗=0𝑔_𝑗^(𝑙)⊗𝑝_𝑗(𝑥)≤2𝛾²_5𝑛+2^𝑛−1_𝑗=0𝑎_𝑗𝑝_𝑗(|𝑥|). (60)

From Definitions 4 and 6 and Theorem 12, we easily get the following corollary.

Corollary 13. Let 𝑝(𝑥)=∑^𝑛_𝑗=0𝑎_𝑗𝑝_𝑗(𝑥) be a Legendre series of degree 𝑛 and 𝑥 a floating point value. The relative error bound of the compensated Clenshaw algorithm is 𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤𝑝,𝑥−𝑝(𝑥)𝑝(𝑥)≤𝑢+2𝛾²_5𝑛+2𝑐𝑜𝑛𝑑_𝑟𝑒𝑙𝑝,𝑥. (61)

4. Numerical Results

All our experiments are performed using IEEE-754 double precision as working precision. Here, we consider the polynomials in Legendre basis with real coefficients and floating point entry 𝑥. All the programs about accuracy measurements have been written in MATLAB R2012b and that about timing measurements have been written in C code on a 2.53-GHz Intel Core i5 laptop.

4.1. Evaluation of the Polynomial in Legendre Basis

In order to construct an ill-conditioned polynomial, we consider the evaluation of the polynomial in Legendre basis 𝑝(𝑥)=∑¹⁷_𝑖=0𝑎_𝑖𝑝_𝑖(𝑥) converted by the polynomials 𝑃(𝑥)=(𝑥−0.75)⁷(𝑥−1)¹¹ in the neighborhood of its multiple roots 0.75 and 1. We use the 𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤, 𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithms and the Symbolic Toolbox to evaluate the polynomial in Legendre basis. In order to observe the perturbation of polynomial coefficients clearly, we propose the 𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm in Appendix A (see Algorithm 7) to obtain the coefficients of the Legendre series. To decrease the perturbation of the coefficients we also propose the 𝐷𝐷𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm in Appendix B (see Algorithm 16) to store the coefficients in double-double format. That is, the coefficients of the polynomial evaluated, which are obtained by the 𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm and 𝐷𝐷𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm, are in double format and double-double format, respectively. In this experiment we evaluate the polynomials for 400 equally spaced points in the intervals [0.68,1.15], [0.7485,0.7515], and [0.993,1.007].

function 𝑏=𝐶𝑜𝑛𝑣𝑒𝑟𝑡(𝑎) 𝑎,𝑏 are the vectors of monomial polynomial and orthogonal polynomial’s coefficients, respectively 𝑏(𝑛)=0, 𝑏0(𝑛)=𝑎𝑛, for 𝑗=𝑛−1∶−1∶0  for 𝑘=0∶1∶𝑛−𝑗   𝑏0(𝑗)=𝑎𝑗−𝐶1𝐴1𝑏1(𝑗+1) (𝑘=0)   𝑏𝑘(𝑗)=1𝐴𝑘−1𝑏(𝑗+1)𝑘−1−𝐶𝑘+1𝐴𝑘+1𝑏(𝑗+1)𝑘+1  (𝑘=1,…,𝑛−𝑗−2)   𝑏𝑘(𝑗)=1𝐴𝑘−1𝑏(𝑗+1)𝑘−1  (𝑘≥𝑛−𝑗−1)  end end 𝑏=𝑏(0)

Algorithm 7 

The conversion algorithm from power basis to Legendre basis.

From Figures 1-2 we observe that the polynomial evaluated by 𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithm (on the top figure) is oscillating, and the compensated algorithm is more smooth drawing. The polynomials we evaluated by Symbolic Toolbox (on the bottom of Figures 1-2) are different because the perturbations of coefficients obtained by the 𝐷𝐷𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm are smaller than those by the 𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm. We can see that the accuracy of the polynomials evaluated by the compensated algorithm is the same with evaluated by Symbolic Toolbox in Figure 1. The polynomials in Legendre basis evaluated by the compensated algorithm are much more smooth drawing and just a little oscillation in the intervals [0.7485,0.7515] in Figure 2. In fact, if we use the Symbolic Toolbox to get the polynomial coefficients, the oscillation will be smaller than it is in Figure 2. However, this method is expensive. We just need to use the 𝐷𝐷𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm to get the coefficients; the result obtained by the 𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithm is almost the same as that by using the Symbolic Toolbox in working precision.

Figure 1 

Evaluation of polynomial converted from 𝑝(𝑥)=(𝑥−0.75)7(𝑥−1)11 by the 𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm in Legendre basis in the neighborhood of its multiple roots, using the Clenshaw algorithm (up), the 𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithm (middle), and Symbolic Toolbox (down).

Figure 2 

Evaluation of polynomial converted from 𝑝(𝑥)=(𝑥−0.75)7(𝑥−1)11 by the 𝐷𝐷𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm in Legendre basis in the neighborhood of its multiple roots, using the Clenshaw algorithm (up), the 𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithm (middle), and Symbolic Toolbox (down).

4.2. Accuracy of the Compensated Algorithm

The closer to the root, the larger the condition number. Thus, in this experiment, the evaluation is for 120 points near the root 0.75, that is, 𝑥=0.75−1.03^2𝑖−85, for 𝑖= 1 : 40 and 𝑥=0.75−1.13^𝑖−85 for 𝑖= 1 : 80. We compare the compensated algorithm with multiple precision library. Since the working precision is double precision, we choose the double-double arithmetic [22] (see Appendix B) which is the most efficient way to yield a full precision accuracy of evaluating the polynomial in Legendre basis to compare with the compensated algorithm. We evaluate the polynomials by the 𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤, 𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤, and 𝐷𝐷𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithms in Appendix B (see Algorithm 15) and the Symbolic Toolbox, respectively, so that the relative forward errors can be obtained by |𝑝_res(𝑥)−𝑝_sym(𝑥)|/𝑝_sym(𝑥) and the relative error bounds are described from Corollaries 13 and A.2 in Appendix A. Then we propose the relative forward errors of evaluation of the polynomial in Legendre basis in Figure 3. As we can see, the relative errors of the compensated algorithm and double-double arithmetic are both smaller than 𝑢 (𝑢≈1.16×10⁻¹⁶) when the condition number is less than 10¹⁷. And the accuracy of both algorithms is decreasing linearly for the condition number larger than 10¹⁷. However, the 𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithm cannot yield the working precision; the accuracy of which decreases linearly since the condition number is less than 10¹⁷. When the condition number is lager than 10¹⁷, the Clenshaw algorithm cannot obtain even one significant bit.

Figure 3 

Accuracy of evaluation of polynomial converted from 𝑝(𝑥)=(𝑥−0.75)7(𝑥−1)11 by the 𝐷𝐷𝐶𝑜𝑛𝑣𝑒𝑟𝑡 algorithm in Legendre basis with respect to the condition number.

4.3. Time Performances

We can easily know that the 𝑇𝑤𝑜𝑆𝑢𝑚, 𝑇𝑤𝑜𝑃𝑟𝑜𝑑, and 𝐹𝑎𝑠𝑡𝑇𝑤𝑜𝑆𝑢𝑚 algorithms in Appendix B (Algorithm 8) require 6, 17, and 3 flops, respectively. Then we obtain the computational cost of the 𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤, 𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤, and 𝐷𝐷𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤 algorithms:(i)𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤: 5n-2 flops;(ii)𝐶𝑜𝑚𝑝𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤: 79n-29 flops;(iii)𝐷𝐷𝐶𝑙𝑒𝑛𝑠ℎ𝑎𝑤: 110n-44 flops.

function [𝑥,𝑦]=𝐹𝑎𝑠𝑡𝑇𝑤𝑜𝑆𝑢𝑚(𝑎,𝑏) 𝑥=𝑎⊕𝑏 𝑦=(𝑎⊖𝑥)⊕𝑏

Algorithm 8 

EFT of the sum of two floating point numbers (|𝑎|≥|