# Accurate Evaluation of Polynomials in Legendre Basis

Authors: Peibing Du; Hao Jiang; Lizhi Cheng
Journal: Journal of Applied Mathematics
 (2014)

---

## 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.


---

## Body

## 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 p
k
(
x
),
(1)
p
0
(
x
)
=
1
p
1
(
x
)
=
x
p
k
+
1
(
x
)
=
2
k
+
1
k
+
1
x
p
k
(
x
)
-
k
k
+
1
p
k
-
1
(
x
)
(
k
>
0
)
.
The polynomial represented in Legendre basis is p
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
), where a
j
∈
R and p
j
(
x
) 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 p
^
(
x
) by the Clenshaw algorithm verifies
(2)
|
p
(
x
)
-
p
^
(
x
)
|
p
(
x
)
≤
cond
(
p
,
x
)
×
o
(
u
)
.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(3)
|
p
(
x
)
-
p
^
(
x
)
|
p
(
x
)
≤
o
(
u
)
+
cond
(
p
,
x
)
×
o
(
u
2
)
.The paper is organized as follows. Section2 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. Leto
p
∈
{
⊕
,
⊖
,
⊗
,
⊘
} represent the floating point computation; then the computation obeys the model
(4)
a
o
p
b
=
(
a
∘
b
)
(
1
+
ε
1
)
=
(
a
∘
b
)
(
1
+
ε
2
)
,
where ∘
∈
{
+
,
-
,
×
,
÷
} and |
ε
1
|
,
|
ε
2
|
≤
u (u is the round-off unit). We also assume that the computed result of a
∈
R in floating point arithmetic is denoted by a
^ and the set of all floating point numbers is denoted by F. The following definition will be used in error analysis (see more details in [17]).Definition 1.
One defines(5)
〈
n
〉
:
=
∏
i
=
1
n
‍
(
1
+
δ
i
)
ρ
i
=
1
+
θ
n
,
where |
δ
i
|
≤
u
,
ρ
i
=
±
1 with ∀
i
=
1,2
,
…
,
n, |
θ
n
|
≤
γ
n
:
=
(
n
u
/
(
1
-
n
u
)
)
=
n
u
+
O
(
u
2
), and n
u
<
1.There are three classic properties which will also be used in error analysis:(i)
〈
k
〉
〈
j
〉
=
〈
k
+
j
〉;
(ii)
γ
k
<
γ
k
+
1;
(iii)
γ
k
+
γ
j
+
γ
k
γ
j
≤
γ
k
+
j.


### 2.2. Accurate Sum and Product

Leta
,
b
∈
F, and no overflow nor underflow occurs. The transformation (
a
,
b
)
→
(
x
,
y
) is regarded as an error-free transformation (EFT) that causes y
∈
F to exist such that a
∘
b
=
x
+
y, x
=
a
o
p
b.Let us show the error-free transformations of the sum and product of two floating point numbers in Algorithms1-3 which are the T
w
o
S
u
m algorithm by Knuth [18] and the T
w
o
P
r
o
d algorithm by Dekker [19], respectively.Algorithm 1:Sum of two floating point numbers.
function[
x
,
y
]
=
T
w
o
S
u
m
(
a
,
b
)
x
=
a
⊕
b
z
=
x
⊖
a
y
=
(
a
⊖
(
x
⊖
z
)
)
⊕
(
b
⊖
z
)Algorithm 2:Split of a floating point number into two parts.
function[
x
,
y
]
=
S
p
l
i
t
(
a
)
c
=
f
a
c
t
o
r
⊗
a (in double precision factor =
2
27
+
1)
x
=
c
⊖
(
c
⊖
a
)
y
=
a
⊖
xAlgorithm 3:Product of two floating point numbers.
function[
x
,
y
]
=
T
w
o
P
r
o
d
(
a
,
b
)
x
=
a
⊗
b
[
a
h
,
a
l
]
=
S
p
l
i
t
(
a
)
[
b
h
,
b
l
]
=
S
p
l
i
t
(
b
)
y
=
a
l
⊗
b
l
⊖
(
(
(
x
⊖
a
h
⊗
b
h
)
⊖
a
l
⊗
b
h
)
⊖
a
h
⊗
b
l
)Algorithms1–3 satisfy the Theorem 2.Theorem 2 (see [11]).
Fora
,
b
∈
F and x
,
y
∈
F, T
w
o
S
u
m and T
w
o
P
r
o
d verify
(6)
[
x
,
y
]
=
T
w
o
S
u
m
(
a
,
b
)
,
x
+
y
=
a
+
b
,
|
y
|
≤
min
⁡
{
u
|
x
|
,
u
|
a
+
b
|
}
,
[
x
,
y
]
=
T
w
o
P
r
o
d
(
a
,
b
)
,
x
+
y
=
a
b
,
|
y
|
≤
min
⁡
{
u
|
x
|
,
u
|
a
b
|
}
.We present the compensated algorithm for the product of three floating point numbers in Algorithm4, which refers to [20].Algorithm 4:Product of three floating point numbers.
function[
x
,
y
]
=
T
h
r
e
e
P
r
o
d
(
a
,
b
,
c
)
[
t
,
h
]
=
T
w
o
P
r
o
d
(
a
,
b
)
[
x
,
e
]
=
T
w
o
P
r
o
d
(
t
,
c
)
y
=
c
⊗
h
⊕
eAccording to Theorem2, we have |
h
|
≤
u
|
a
b
| and t
+
h
=
a
b. Hence, |
e
|
≤
u
|
(
a
b
-
h
)
c
|
≤
u
|
a
b
c
|+u
|
h
c
|
≤
u
|
a
b
c
|
+
u
2
|
a
b
c
|; then |
c
h
+
e
|
≤
|
c
h
|
+
|
e
|
≤
2
u
|
a
b
c
|
+
u
2
|
a
b
c
|. According to Lemma 3.2 in [20], we can propose the error bound of Algorithm 4 in Theorem 3.Theorem 3.
Fora
,
b
,
c
∈
F and x
,
y
,
τ
∈
F, T
h
r
e
e
P
r
o
d verifies
(7)
[
x
,
y
]
=
T
h
r
e
e
P
r
o
d
(
a
,
b
,
c
)
,
x
+
y
=
a
b
c
-
τ
,
|
τ
|
≤
γ
2
γ
6
|
a
b
c
|
,
|
y
|
≤
2
u
|
a
b
c
|
+
14
u
2
|
a
b
c
|
.Proof.
According to Lemma 3.2 in [20], we have
(8)
|
τ
|
≤
γ
2
γ
6
|
a
b
c
|
.
In Algorithm 4, y
+
τ
=
c
h
+
e; then
(9)
|
y
|
=
|
y
+
τ
-
τ
|
≤
|
c
h
+
e
|
+
|
τ
|
≤
2
u
|
a
b
c
|
+
u
2
|
a
b
c
|
+
γ
2
γ
6
|
a
b
c
|
≤
2
u
|
a
b
c
|
+
14
u
2
|
a
b
c
|
.


### 2.3. The Clenshaw Algorithm and Condition Number

The standard general algorithm for the evaluation of polynomial in Legendre basisp
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
) is the Clenshaw algorithm [5]. We recall it in Algorithm 5.Algorithm 5:The Clenshaw algorithm for finite Legendre series.
functionp
(
x
)
=
C
l
e
n
s
h
a
w
(
p
,
x
)
b
n
+
2
=
b
n
+
1
=
0
forj
=
n
:
-
1
:
0
b
j
=
2
j
+
1
j
+
1
x
b
j
+
1
-
j
+
1
j
+
2
b
j
+
2
+
a
j
end
p
(
x
)
=
b
0Barrio 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.
Letp
k
(
x
) be Legendre polynomial. We define the absolute Legendre polynomial p
~
k
(
x
) which is associated with p
k
(
x
) and satisfies
(10)
p
~
0
(
x
)
=
1
p
~
1
(
x
)
=
x
p
~
k
(
x
)
=
2
k
-
1
k
x
p
~
k
-
1
(
x
)
+
k
k
+
1
p
~
k
-
2
(
x
)
(
k
≥
2
)
,
where |
p
k
(
x
)
|
≤
p
~
k
(
x
), ∀
x
≥
0.The absolute Legendre polynomial satisfies the following property.Lemma 5.
Letp
~
k
(
x
) be the absolute Legendre polynomial. Then we have
(11)
p
~
k
(
x
)
p
~
j
(
x
)
≤
p
~
k
+
j
(
x
)
.Proof .
LetA
k
=
(
2
k
+
1
)
/
(
k
+
1
) and k
≥
2. From Definition 4, we have
(12)
p
~
k
(
x
)
=
A
k
-
1
x
p
~
k
-
1
(
x
)
+
k
k
+
1
p
~
k
-
2
(
x
)
=
(
∏
i
=
0
k
-
1
‍
A
i
)
x
k
+
⋯
;
then
(13)
p
~
k
(
x
)
p
~
j
(
x
)
=
(
∏
i
=
0
k
-
1
‍
A
i
)
(
∏
i
=
0
j
-
1
‍
A
i
)
x
k
+
j
+
⋯
.
Thus the equivalent form of (11) is
(14)
(
∏
i
=
0
k
-
1
‍
A
i
)
(
∏
i
=
0
j
-
1
‍
A
i
)
≤
∏
i
=
0
k
+
j
-
1
‍
A
i
.
Since A
k is increasing with k, we obtain
(15)
∏
i
=
0
j
-
1
‍
A
i
≤
∏
i
=
k
k
+
j
-
1
‍
A
i
.
So we finally obtain p
~
k
(
x
)
p
~
j
(
x
)
≤
p
~
k
+
j
(
x
).Following Definition4, we introduce the condition number for the evaluation of polynomials in Legendre basis [21].Definition 6.
Letp
(
x
)
=
∑
k
=
0
n
‍
a
k
p
k
(
x
), where p
k
(
x
) is Legendre polynomial. Let p
~
k
(
x
) be the absolute Legendre polynomial. Then the absolute condition number is
(16)
cond
abs
(
p
,
x
)
=
p
~
(
|
x
|
)
=
∑
k
=
0
n
‍
|
a
k
|
p
~
k
(
|
x
|
)
,
and the relative condition number is
(17)
cond
rel
⁡
(
p
,
x
)
=
p
~
(
|
x
|
)
|
p
(
x
)
|
=
∑
k
=
0
n
‍
|
a
k
|
p
~
k
(
|
x
|
)
|
p
(
x
)
|
.


## 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. Leto
p
∈
{
⊕
,
⊖
,
⊗
,
⊘
} represent the floating point computation; then the computation obeys the model
(4)
a
o
p
b
=
(
a
∘
b
)
(
1
+
ε
1
)
=
(
a
∘
b
)
(
1
+
ε
2
)
,
where ∘
∈
{
+
,
-
,
×
,
÷
} and |
ε
1
|
,
|
ε
2
|
≤
u (u is the round-off unit). We also assume that the computed result of a
∈
R in floating point arithmetic is denoted by a
^ and the set of all floating point numbers is denoted by F. The following definition will be used in error analysis (see more details in [17]).Definition 1.
One defines(5)
〈
n
〉
:
=
∏
i
=
1
n
‍
(
1
+
δ
i
)
ρ
i
=
1
+
θ
n
,
where |
δ
i
|
≤
u
,
ρ
i
=
±
1 with ∀
i
=
1,2
,
…
,
n, |
θ
n
|
≤
γ
n
:
=
(
n
u
/
(
1
-
n
u
)
)
=
n
u
+
O
(
u
2
), and n
u
<
1.There are three classic properties which will also be used in error analysis:(i)
〈
k
〉
〈
j
〉
=
〈
k
+
j
〉;
(ii)
γ
k
<
γ
k
+
1;
(iii)
γ
k
+
γ
j
+
γ
k
γ
j
≤
γ
k
+
j.


## 2.2. Accurate Sum and Product

Leta
,
b
∈
F, and no overflow nor underflow occurs. The transformation (
a
,
b
)
→
(
x
,
y
) is regarded as an error-free transformation (EFT) that causes y
∈
F to exist such that a
∘
b
=
x
+
y, x
=
a
o
p
b.Let us show the error-free transformations of the sum and product of two floating point numbers in Algorithms1-3 which are the T
w
o
S
u
m algorithm by Knuth [18] and the T
w
o
P
r
o
d algorithm by Dekker [19], respectively.Algorithm 1:Sum of two floating point numbers.
function[
x
,
y
]
=
T
w
o
S
u
m
(
a
,
b
)
x
=
a
⊕
b
z
=
x
⊖
a
y
=
(
a
⊖
(
x
⊖
z
)
)
⊕
(
b
⊖
z
)Algorithm 2:Split of a floating point number into two parts.
function[
x
,
y
]
=
S
p
l
i
t
(
a
)
c
=
f
a
c
t
o
r
⊗
a (in double precision factor =
2
27
+
1)
x
=
c
⊖
(
c
⊖
a
)
y
=
a
⊖
xAlgorithm 3:Product of two floating point numbers.
function[
x
,
y
]
=
T
w
o
P
r
o
d
(
a
,
b
)
x
=
a
⊗
b
[
a
h
,
a
l
]
=
S
p
l
i
t
(
a
)
[
b
h
,
b
l
]
=
S
p
l
i
t
(
b
)
y
=
a
l
⊗
b
l
⊖
(
(
(
x
⊖
a
h
⊗
b
h
)
⊖
a
l
⊗
b
h
)
⊖
a
h
⊗
b
l
)Algorithms1–3 satisfy the Theorem 2.Theorem 2 (see [11]).
Fora
,
b
∈
F and x
,
y
∈
F, T
w
o
S
u
m and T
w
o
P
r
o
d verify
(6)
[
x
,
y
]
=
T
w
o
S
u
m
(
a
,
b
)
,
x
+
y
=
a
+
b
,
|
y
|
≤
min
⁡
{
u
|
x
|
,
u
|
a
+
b
|
}
,
[
x
,
y
]
=
T
w
o
P
r
o
d
(
a
,
b
)
,
x
+
y
=
a
b
,
|
y
|
≤
min
⁡
{
u
|
x
|
,
u
|
a
b
|
}
.We present the compensated algorithm for the product of three floating point numbers in Algorithm4, which refers to [20].Algorithm 4:Product of three floating point numbers.
function[
x
,
y
]
=
T
h
r
e
e
P
r
o
d
(
a
,
b
,
c
)
[
t
,
h
]
=
T
w
o
P
r
o
d
(
a
,
b
)
[
x
,
e
]
=
T
w
o
P
r
o
d
(
t
,
c
)
y
=
c
⊗
h
⊕
eAccording to Theorem2, we have |
h
|
≤
u
|
a
b
| and t
+
h
=
a
b. Hence, |
e
|
≤
u
|
(
a
b
-
h
)
c
|
≤
u
|
a
b
c
|+u
|
h
c
|
≤
u
|
a
b
c
|
+
u
2
|
a
b
c
|; then |
c
h
+
e
|
≤
|
c
h
|
+
|
e
|
≤
2
u
|
a
b
c
|
+
u
2
|
a
b
c
|. According to Lemma 3.2 in [20], we can propose the error bound of Algorithm 4 in Theorem 3.Theorem 3.
Fora
,
b
,
c
∈
F and x
,
y
,
τ
∈
F, T
h
r
e
e
P
r
o
d verifies
(7)
[
x
,
y
]
=
T
h
r
e
e
P
r
o
d
(
a
,
b
,
c
)
,
x
+
y
=
a
b
c
-
τ
,
|
τ
|
≤
γ
2
γ
6
|
a
b
c
|
,
|
y
|
≤
2
u
|
a
b
c
|
+
14
u
2
|
a
b
c
|
.Proof.
According to Lemma 3.2 in [20], we have
(8)
|
τ
|
≤
γ
2
γ
6
|
a
b
c
|
.
In Algorithm 4, y
+
τ
=
c
h
+
e; then
(9)
|
y
|
=
|
y
+
τ
-
τ
|
≤
|
c
h
+
e
|
+
|
τ
|
≤
2
u
|
a
b
c
|
+
u
2
|
a
b
c
|
+
γ
2
γ
6
|
a
b
c
|
≤
2
u
|
a
b
c
|
+
14
u
2
|
a
b
c
|
.


## 2.3. The Clenshaw Algorithm and Condition Number

The standard general algorithm for the evaluation of polynomial in Legendre basisp
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
) is the Clenshaw algorithm [5]. We recall it in Algorithm 5.Algorithm 5:The Clenshaw algorithm for finite Legendre series.
functionp
(
x
)
=
C
l
e
n
s
h
a
w
(
p
,
x
)
b
n
+
2
=
b
n
+
1
=
0
forj
=
n
:
-
1
:
0
b
j
=
2
j
+
1
j
+
1
x
b
j
+
1
-
j
+
1
j
+
2
b
j
+
2
+
a
j
end
p
(
x
)
=
b
0Barrio 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.
Letp
k
(
x
) be Legendre polynomial. We define the absolute Legendre polynomial p
~
k
(
x
) which is associated with p
k
(
x
) and satisfies
(10)
p
~
0
(
x
)
=
1
p
~
1
(
x
)
=
x
p
~
k
(
x
)
=
2
k
-
1
k
x
p
~
k
-
1
(
x
)
+
k
k
+
1
p
~
k
-
2
(
x
)
(
k
≥
2
)
,
where |
p
k
(
x
)
|
≤
p
~
k
(
x
), ∀
x
≥
0.The absolute Legendre polynomial satisfies the following property.Lemma 5.
Letp
~
k
(
x
) be the absolute Legendre polynomial. Then we have
(11)
p
~
k
(
x
)
p
~
j
(
x
)
≤
p
~
k
+
j
(
x
)
.Proof .
LetA
k
=
(
2
k
+
1
)
/
(
k
+
1
) and k
≥
2. From Definition 4, we have
(12)
p
~
k
(
x
)
=
A
k
-
1
x
p
~
k
-
1
(
x
)
+
k
k
+
1
p
~
k
-
2
(
x
)
=
(
∏
i
=
0
k
-
1
‍
A
i
)
x
k
+
⋯
;
then
(13)
p
~
k
(
x
)
p
~
j
(
x
)
=
(
∏
i
=
0
k
-
1
‍
A
i
)
(
∏
i
=
0
j
-
1
‍
A
i
)
x
k
+
j
+
⋯
.
Thus the equivalent form of (11) is
(14)
(
∏
i
=
0
k
-
1
‍
A
i
)
(
∏
i
=
0
j
-
1
‍
A
i
)
≤
∏
i
=
0
k
+
j
-
1
‍
A
i
.
Since A
k is increasing with k, we obtain
(15)
∏
i
=
0
j
-
1
‍
A
i
≤
∏
i
=
k
k
+
j
-
1
‍
A
i
.
So we finally obtain p
~
k
(
x
)
p
~
j
(
x
)
≤
p
~
k
+
j
(
x
).Following Definition4, we introduce the condition number for the evaluation of polynomials in Legendre basis [21].Definition 6.
Letp
(
x
)
=
∑
k
=
0
n
‍
a
k
p
k
(
x
), where p
k
(
x
) is Legendre polynomial. Let p
~
k
(
x
) be the absolute Legendre polynomial. Then the absolute condition number is
(16)
cond
abs
(
p
,
x
)
=
p
~
(
|
x
|
)
=
∑
k
=
0
n
‍
|
a
k
|
p
~
k
(
|
x
|
)
,
and the relative condition number is
(17)
cond
rel
⁡
(
p
,
x
)
=
p
~
(
|
x
|
)
|
p
(
x
)
|
=
∑
k
=
0
n
‍
|
a
k
|
p
~
k
(
|
x
|
)
|
p
(
x
)
|
.


## 3. Compensated Algorithm for Evaluating Polynomials

In this section, we exhibit the exact round-off errors generated by the Clenshaw algorithm withEFT. 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)
A
=
A
(
h
)
+
A
(
l
)
+
A
(
m
)
,
where A
(
h
)
,
A
(
l
)
∈
F, A
,
A
(
m
)
∈
R and |
A
(
l
)
|
≤
u
|
A
(
h
)
|
,
|
A
(
m
)
|
≤
u
|
A
(
l
)
|. Let A
k
=
(
2
k
+
1
)
/
(
k
+
1
) and C
k
=
-
(
k
/
(
k
+
1
)
). Then we describe the recurrence relation at jth step of Algorithm 5 for theoretical computation as
(19)
b
j
=
(
A
j
(
h
)
+
A
j
(
l
)
+
A
j
(
m
)
)
x
b
j
+
1
+
(
C
j
+
1
(
h
)
+
C
j
+
1
(
l
)
+
C
j
+
1
(
m
)
)
b
j
+
2
+
(
a
j
(
h
)
+
a
j
(
l
)
+
a
j
(
m
)
)
.
For numerical computation associated with (19), it is
(20)
b
^
j
=
A
j
(
h
)
⊗
x
⊗
b
^
j
+
1
⊕
C
j
+
1
(
h
)
⊗
b
^
j
+
2
⊕
a
j
(
h
)
.Remark 7.
LetA
∈
R. We split a coefficient A like (18); then the representation of splitting is unique and A
^
=
A
(
h
)
⊕
A
(
l
), A
^
=
A
〈
1
〉.LetA
j
=
A
j
(
h
) and C
j
=
C
j
(
h
) in Algorithm 5. Since every elementary floating point operation in Algorithm 5 causes round-off errors in numerical computation, we apply E
F
T and the T
h
r
e
e
P
r
o
d algorithms to take notes of all round-off errors and obtain
(21)
[
r
,
α
j
]
=
T
h
r
e
e
P
r
o
d
(
A
j
(
h
)
,
x
,
b
^
j
+
1
)
,
[
s
,
β
j
]
=
T
w
o
P
r
o
d
(
C
j
+
1
(
h
)
,
b
^
j
+
2
)
,
[
t
,
η
j
]
=
T
w
o
S
u
m
(
r
,
s
)
,
[
b
^
j
,
ξ
j
]
=
T
w
o
S
u
m
(
t
,
a
j
(
h
)
)
.
The sum of the perturbation and the round-off errors of the recurrence relation at jth step is
(22)
w
j
=
(
A
j
(
l
)
+
A
j
(
m
)
)
x
b
^
j
+
1
+
(
C
j
+
1
(
l
)
+
C
j
+
1
(
m
)
)
b
^
j
+
2
+
(
a
j
(
l
)
+
a
j
(
m
)
)
+
α
j
+
τ
j
+
β
j
+
η
j
+
ξ
j
,
where r
+
α
j
+
τ
j
=
A
j
(
h
)
x
b
^
j
+
1, τ
j is defined in Theorem 3. Then we obtain the following theorem.Theorem 8.
Letp
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
) be a Legendre series of degree n, and let x be a floating point value. One assumes that A
j
=
A
j
(
h
) and C
j
=
C
j
(
h
) in Algorithm 5 and b
^
0 is the numerical result; w
j is described in (22) for j
=
0,1
,
…
,
n
-
1. Then one obtains
(23)
∑
j
=
0
n
-
1
‍
w
j
p
j
(
x
)
+
b
^
0
=
∑
j
=
0
n
‍
a
j
p
j
(
x
)
.Proof.
The perturbation of Algorithm5 is
(24)
(
A
j
(
l
)
+
A
j
(
m
)
)
x
b
^
j
+
1
+
(
C
j
+
1
(
l
)
+
C
j
+
1
(
m
)
)
b
^
j
+
2
+
(
a
j
(
l
)
+
a
j
(
m
)
)
.
Let A
j
=
A
j
(
h
) and C
j
=
C
j
(
h
) in Algorithm 5. From Theorems 2–3 and (21), we can easily obtain
(25)
A
j
(
h
)
b
^
j
+
1
x
=
r
+
α
j
+
τ
j
,
C
j
+
1
(
h
)
b
^
j
+
2
=
s
+
β
j
,
r
+
s
=
t
+
η
j
t
+
a
j
(
h
)
=
b
^
j
+
ξ
j
.
Lets
j
=
α
j
+
τ
j
+
β
j
+
η
j
+
ξ
j for j
=
0,1
,
…
,
n
-
1; then
(26)
b
^
j
+
s
j
=
A
j
(
h
)
b
^
j
+
1
x
+
C
j
+
1
(
h
)
b
^
j
+
2
+
a
j
(
h
)
,
j
=
0,1
,
…
,
n
-
1
.
From (18), (19), (22), and (24), we have
(27)
b
^
j
+
w
j
=
A
j
b
^
j
+
1
x
+
C
j
+
1
b
^
j
+
2
+
a
j
,
j
=
0,1
,
…
,
n
-
1
.
Thus(28)
b
^
0
=
∑
j
=
0
n
‍
(
a
j
-
w
j
)
p
j
(
x
)
;
then
(29)
∑
j
=
0
n
-
1
‍
w
j
p
j
(
x
)
+
b
^
0
=
∑
j
=
0
n
‍
a
j
p
j
(
x
)
.To devise a compensated algorithm of Algorithm5 and give its error bound, we need the following lemmas.Lemma 9.
Letp
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
) be a Legendre series of degree n, x a floating point value, and b
^
0 the numerical result of Algorithm 5. Then
(30)
|
C
l
e
n
s
h
a
w
(
p
,
x
)
-
p
(
x
)
|
≤
γ
5
n
-
1
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.Proof.
Applying the standard model of floating point arithmetic, from Algorithm5 and Remark 7, we have
(31)
b
^
n
-
1
=
A
n
-
1
b
^
n
x
〈
4
〉
+
a
n
-
1
〈
1
〉
,
b
^
n
-
2
=
A
n
-
2
b
^
n
-
1
x
〈
5
〉
+
C
n
-
1
b
^
n
〈
4
〉
+
a
n
-
2
〈
1
〉
,
⋮
b
^
0
=
A
0
b
^
1
x
〈
5
〉
+
C
1
b
^
2
〈
4
〉
+
a
0
〈
1
〉
.
Thus we obtain(32)
b
^
0
=
〈
4
+
5
(
n
-
1
)
〉
(
∏
i
=
0
n
-
1
‍
A
i
)
x
n
a
n
+
⋯
+
a
0
〈
1
〉
=
〈
5
n
-
1
〉
(
∏
i
=
0
n
-
1
‍
A
i
)
x
n
a
n
+
⋯
+
a
0
〈
1
〉
.
Then, from Definition1, we get
(33)
|
b
^
0
-
p
(
x
)
|
=
|
(
∏
i
=
0
n
-
1
‍
A
i
)
x
n
a
n
θ
5
n
-
1
+
⋯
+
a
0
θ
1
|
≤
γ
5
n
-
1
(
(
∏
i
=
0
n
-
1
‍
A
i
)
|
x
|
n
|
a
n
|
+
⋯
+
|
a
0
|
)
=
γ
5
n
-
1
∑
j
=
0
n
‍
|
a
j
|
P
~
j
(
|
x
|
)
.Lemma 10.
Letp
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
) be a Legendre series of degree n and x a floating point value. We assume that σ
j
≤
ω
j
(
A
j
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
|
|
b
^
j
+
2
|
+
|
a
j
|
), where ω
j is real numbers; then
(34)
∑
j
=
0
n
-
1
‍
σ
j
p
~
j
(
|
x
|
)
≤
n
ω
j
(
1
+
γ
5
(
n
-
1
)
)
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.Proof.
According to (31), we get
(35)
|
b
^
j
|
≤
(
1
+
γ
5
)
(
A
j
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
|
|
b
^
j
+
2
|
+
|
a
j
|
)
.
Since σ
j
≤
ω
j
(
A
j
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
|
|
b
^
j
+
2
|
+
|
a
j
|
), from Definition 1, we have
(36)
σ
j
≤
ω
j
(
1
+
γ
5
(
n
-
j
-
1
)
)
Σ
k
=
j
n
|
a
k
|
p
~
k
-
j
(
|
x
|
)
;
then
(37)
∑
j
=
0
n
-
1
‍
σ
j
p
~
j
(
|
x
|
)
≤
ω
j
(
1
+
γ
5
(
n
-
1
)
)
×
∑
j
=
0
n
-
1
‍
(
∑
k
=
j
n
‍
|
a
k
|
p
~
k
-
j
(
|
x
|
)
)
p
~
j
(
|
x
|
)
.
According to Lemma 5, we obtain
(38)
∑
j
=
0
n
-
1
‍
σ
j
p
~
j
(
|
x
|
)
≤
ω
j
(
1
+
γ
5
(
n
-
1
)
)
∑
j
=
0
n
-
1
‍
∑
k
=
j
n
‍
|
a
k
|
p
~
k
(
|
x
|
)
≤
n
ω
j
(
1
+
γ
5
(
n
-
1
)
)
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.Among the perturbed and round-off errors in Algorithm5, 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 g
j
(
l
)
=
A
j
(
l
)
x
b
^
j
+
1
+
C
j
+
1
(
l
)
b
^
j
+
2
+
a
j
(
l
). According to A
(
l
)
≤
u
A
(
h
) and Remark 7, we obtain
(39)
|
g
j
(
l
)
|
≤
u
(
1
+
γ
1
)
(
A
j
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
|
|
b
^
j
+
2
|
+
|
a
j
|
)
.
Then let ω
j
=
u
(
1
+
γ
1
) and from Lemma 10, taking into account that (
n
+
1
)
u
(
1
+
γ
5
n
+
1
)
≤
(
5
n
+
1
)
u
(
1
+
γ
5
n
+
1
)
=
γ
5
n
+
1, we have
(40)
|
∑
j
=
0
n
‍
g
j
(
l
)
p
j
(
x
)
|
≤
γ
5
n
+
1
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.γ
5
n
+
1 is O
(
u
), so this coefficient perturbation may influence the accuracy; we need to consider it in our compensated algorithm.Similarly, we letg
j
(
m
)
=
A
j
(
m
)
x
b
^
j
+
1
+
C
j
+
1
(
m
)
b
^
j
+
2
+
a
j
(
m
); then
(41)
|
∑
j
=
0
n
‍
g
j
(
m
)
p
j
(
x
)
|
≤
u
γ
5
n
+
1
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.u
γ
5
n
+
1 is O
(
u
2
), so this coefficient perturbation does not influence the accuracy.Remark 11.
Whenj
=
n, the nth step of Algorithm 5 is b
^
n
=
a
n
(
h
), so that we only need to consider the perturbation of coefficient a
n.Next we deduce the round-off error bound. Lets
j
=
α
j
+
β
j
+
η
j
+
ξ
j. According to Theorems 2 and 3 and Remark 7, we obtain
(42)
|
α
j
|
≤
(
2
u
+
14
u
2
)
|
A
j
(
h
)
x
b
^
j
+
1
|
≤
(
2
u
+
14
u
2
)
(
1
+
γ
1
)
|
A
j
x
b
^
j
+
1
|
;
|
β
j
|
≤
u
|
C
j
+
1
(
h
)
b
^
j
+
2
|
≤
u
(
1
+
γ
1
)
|
C
j
+
1
b
^
j
+
2
|
;
|
η
j
|
≤
u
|
A
j
(
h
)
x
b
^
j
+
1
〈
3
〉
+
C
j
+
1
(
h
)
b
^
j
+
2
〈
2
〉
|
≤
u
(
1
+
γ
3
)
|
A
j
x
b
^
j
+
1
+
C
j
+
1
b
^
j
+
2
|
;
|
ξ
j
|
≤
u
|
A
j
(
h
)
x
b
^
j
+
1
〈
4
〉
+
C
j
+
1
(
h
)
b
^
j
+
2
〈
3
〉
+
a
j
〈
1
〉
|
≤
u
(
1
+
γ
4
)
|
A
j
x
b
^
j
+
1
+
C
j
+
1
b
^
j
+
2
+
a
j
|
;
then
(43)
|
α
j
|
+
|
β
j
|
+
|
η
j
|
+
|
ξ
j
|
≤
(
4
u
+
14
u
2
)
|
1
+
γ
4
|
×
(
A
j
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
|
|
b
^
j
+
2
|
+
|
a
j
|
)
.Letω
j
=
(
4
u
+
14
u
2
)
|
1
+
γ
4
| and using it into Lemma 10, from (
4
+
14
u
)
n
u
(
1
+
γ
5
n
-
1
)
≤
(
5
n
-
1
)
u
(
1
+
γ
5
n
-
1
)
=
γ
5
n
-
1
(
n
≥
2
), we have
(44)
|
∑
j
=
0
n
-
1
‍
(
|
α
j
|
+
|
β
j
|
+
|
η
j
|
+
|
ξ
j
|
)
p
j
(
x
)
|
≤
γ
5
n
-
1
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.γ
5
n
-
1 is O
(
u
), so this round-off error may influence the accuracy, we also need to consider it in our compensated algorithm.From Theorem3 we have known that round-off error |
τ
j
|
≤
γ
2
γ
6
|
A
j
(
h
)
x
b
^
j
+
1
|
≤
γ
2
γ
6
(
1
+
γ
1
)
|
A
j
x
b
^
j
+
1
|. According to Lemma 10 we let ω
j
=
γ
2
γ
6
(
1
+
γ
1
), from γ
2
≤
2
u
(
1
+
γ
1
) and 2
n
u
(
1
+
γ
5
n
-
3
)
≤
γ
5
n
-
3, we have
(45)
|
∑
j
=
0
n
-
1
‍
|
τ
j
|
p
j
(
x
)
|
≤
γ
6
γ
5
n
-
3
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.γ
6
γ
5
n
-
3 is O
(
u
2
), 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 d
i
v
_
d
_
d algorithm in Appendix B (see Algorithm 13), we only need to split the coefficients into two floating point numbers.ApplyingEFT, the T
h
r
e
e
P
r
o
d algorithm and the d
i
v
_
d
_
d algorithm, considering all errors which may influence the numerical result in working precision in Algorithm 5, we obtain the compensated C
l
e
n
s
h
a
w algorithm in Algorithm 6.Algorithm 6:Compensated Clenshaw algorithm.
functionr
e
s
=
C
o
m
p
C
l
e
n
s
h
a
w
(
p
,
x
)
b
^
n
+
2
=
b
^
n
+
1
=
0
ϵ
b
^
n
+
2
=
ϵ
b
^
n
+
1
=
0
forj
=
n
:
-
1
:
0
[
A
j
(
h
)
,
A
j
(
l
)
]
=
d
i
v
_
d
_
d
(
2
j
+
1
,
j
+
1
)
[
C
j
(
h
)
,
C
j
(
l
)
]
=
d
i
v
_
d
_
d
(
-
j
-
1
,
j
+
2
)
[
r
,
α
j
]
=
T
h
r
e
e
P
r
o
d
(
A
j
(
h
)
,
b
^
j
+
1
,
x
)
[
s
,
β
j
]
=
T
w
o
P
r
o
d
(
C
j
+
1
(
h
)
,
b
^
j
+
2
)
[
t
,
η
j
]
=
T
w
o
S
u
m
(
r
,
s
)
[
b
^
j
,
ξ
j
]
=
T
w
o
S
u
m
(
t
,
a
j
)
s
^
j
=
α
j
⊕
β
j
⊕
η
j
⊕
ξ
j
⊕
(
A
j
(
l
)
⊗
x
⊗
b
^
j
+
1
⊕
C
j
+
1
(
l
)
⊗
b
^
j
+
2
⊕
a
j
(
l
)
)
ϵ
b
^
j
=
A
j
(
h
)
⊗
x
⊗
ϵ
b
^
j
+
1
⊕
C
j
+
1
(
h
)
⊗
ϵ
b
^
j
+
2
⊕
s
^
j
end
r
e
s
=
y
⊕
ϵ
b
^
0Here we give the error bound of Algorithm6.Theorem 12.
Letp
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
) be a Legendre series of degree n
(
n
≥
2
) and x a floating point value. The forward error bound of the compensated Clenshaw algorithm is
(46)
|
C
o
m
p
C
l
e
n
s
h
a
w
(
p
,
x
)
-
p
(
x
)
|
≤
u
|
p
(
x
)
|
+
2
γ
5
n
-
2
2
∑
j
=
0
n
‍
a
j
P
j
(
x
)
.Proof.
From Algorithm6, we obtain
(47)
|
C
o
m
p
C
l
e
n
s
h
a
w
(
p
,
x
)
-
p
(
x
)
|
=
|
(
b
^
0
⊕
ϵ
b
^
0
)
-
p
(
x
)
|
=
|
(
1
+
ϵ
)
(
b
^
0
+
ϵ
b
^
0
)
-
p
(
x
)
|
,
where |
ϵ
|
≤
u.
According to Theorem8, we have p
(
x
)
=
b
^
o
+
ϵ
b
o; then
(48)
|
C
o
m
p
C
l
e
n
s
h
a
w
(
p
,
x
)
-
p
(
x
)
|
=
|
(
1
+
ϵ
)
(
P
(
x
)
-
ϵ
b
o
+
ϵ
b
^
0
)
-
p
(
x
)
|
≤
u
|
P
(
x
)
|
+
(
1
+
u
)
|
ϵ
b
o
-
ϵ
b
^
0
|
.
Next we analyze the bound of |
ϵ
b
o
-
ϵ
b
^
0
|.
Letg
j
(
l
)
=
A
j
(
l
)
x
b
^
j
+
1
+
C
j
+
1
(
l
)
b
^
j
+
2
+
a
j
(
l
); then
(49)
g
^
j
(
l
)
=
A
j
(
l
)
x
b
^
j
+
1
〈
4
〉
+
C
j
+
1
(
l
)
b
^
j
+
2
〈
3
〉
+
a
j
(
l
)
〈
1
〉
.
Thus we obtain(50)
|
∑
j
=
0
n
‍
g
j
(
l
)
p
j
(
x
)
-
∑
j
=
0
n
‍
g
^
j
(
l
)
p
j
(
x
)
|
≤
γ
4
∑
j
=
0
n
‍
(
A
j
(
l
)
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
(
l
)
|
|
b
^
j
+
2
|
+
|
a
j
(
l
)
|
)
p
~
j
(
|
x
|
)
.
According to Lemma9, we get
(51)
|
∑
j
=
0
n
‍
g
^
j
(
l
)
p
j
(
x
)
-
⊕
∑
j
=
0
n
‍
g
^
j
(
l
)
⊗
p
j
(
x
)
|
≤
γ
5
n
-
1
∑
j
=
0
n
‍
|
g
^
j
(
l
)
|
p
~
j
(
|
x
|
)
≤
γ
5
n
-
1
(
1
+
γ
4
)
×
∑
j
=
0
n
‍
(
A
j
(
l
)
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
(
l
)
|
|
b
^
j
+
2
|
+
|
a
j
(
l
)
|
)
p
~
j
(
|
x
|
)
.
From (50), (51) and γ
5
n
-
1
(
1
+
γ
4
)
+
γ
4
≤
γ
5
n
+
3, we derive
(52)
|
∑
j
=
0
n
‍
g
j
(
l
)
p
j
(
x
)
-
⊕
∑
j
=
0
n
‍
g
^
j
(
l
)
⊗
p
j
(
x
)
|
≤
|
∑
j
=
0
n
‍
g
j
(
l
)
p
j
(
x
)
-
∑
j
=
0
n
‍
g
^
j
(
l
)
p
j
(
x
)
|
+
|
∑
j
=
0
n
‍
g
^
j
(
l
)
p
j
(
x
)
-
⊕
∑
j
=
0
n
‍
g
^
j
(
l
)
⊗
p
j
(
x
)
|
≤
γ
5
n
+
3
∑
j
=
0
n
‍
(
A
j
(
l
)
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
(
l
)
|
|
b
^
j
+
2
|
+
|
a
j
(
l
)
|
)
p
~
j
(
|
x
|
)
.
According toA
(
l
)
≤
u
A
(
h
) and Remark 7, we obtain
(53)
A
j
(
l
)
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
(
l
)
|
|
b
^
j
+
2
|
+
|
a
j
(
l
)
|
≤
u
(
1
+
γ
1
)
(
A
j
|
x
|
|
b
^
j
+
1
|
+
|
C
j
+
1
|
|
b
^
j
+
2
|
+
|
a
j
|
)
.
From (53) we use ω
j
=
u
(
1
+
γ
1
) into Lemma 10; taking into account that n
u
(
1
+
γ
5
n
-
4
)
≤
(
5
n
+
1
)
u
(
1
+
γ
5
n
+
1
)
=
γ
5
n
+
1 and γ
5
n
+
1
γ
5
n
+
3
≤
γ
5
n
+
2
2, according to (52), we get
(54)
|
∑
j
=
0
n
-
1
‍
g
j
(
l
)
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
‍
g
^
j
(
l
)
⊗
p
j
(
x
)
|
≤
γ
5
n
+
2
2
∑
j
=
0
n
‍
|
a
j
|
p
~
j
(
|
x
|
)
.
Next we lets
j
=
α
j
+
β
j
+
η
j
+
ξ
j; then
(55)
s
^
j
≤
(
1
+
γ
3
)
(
|
α
j
|
+
|
β
j
|
+
|
η
j
|
+
|
ξ
j
|
)
;
thus
(56)
|
∑
j
=
0
n
-
1
‍
s
j
p
j
(
x
)
-
∑
j
=
0
n
-
1
‍
s
^
j
p
j
(
x
)
|
≤
γ
3
∑
j
=
0
n
-
1
‍
(
|
α
j
|
+
|
β
j
|
+
|
η
j
|
+
|
ξ
j
|
)
p
~
j
(
|
x
|
)
.
According to Lemma9, we obtain
(57)
|
∑
j
=
0
n
-
1
‍
s
^
j
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
‍
s
^
j
⊗
p
j
(
x
)
|
≤
γ
5
(
n
-
1
)
-
1
∑
j
=
0
n
-
1
‍
|
s
^
j
|
p
~
j
(
|
x
|
)
≤
γ
5
(
n
-
1
)
-
1
(
1
+
γ
3
)
×
∑
j
=
0
n
-
1
‍
(
|
α
j
|
+
|
β
j
|
+
|
η
j
|
+
|
ξ
j
|
)
p
~
j
(
|
x
|
)
.
From (56), (57), and γ
5
(
n
-
1
)
-
1
(
1
+
γ
3
)
+
γ
3
≤
γ
5
n
-
3, we derive
(58)
|
∑
j
=
0
n
-
1
‍
s
j
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
‍
s
^
j
⊗
p
j
(
x
)
|
≤
|
∑
j
=
0
n
-
1
‍
s
j
p
j
(
x
)
-
∑
j
=
0
n
-
1
‍
s
^
j
p
j
(
x
)
|
+
|
∑
j
=
0
n
-
1
‍
s
^
j
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
‍
s
^
j
⊗
p
j
(
x
)
|
≤
γ
5
n
-
3
∑
j
=
0
n
-
1
‍
(
|
α
j
|
+
|
β
j
|
+
|
η
j
|
+
|
ξ
j
|
)
p
~
j
(
|
x
|
)
.
From (43), we let ω
j
=
(
4
+
14
u
)
u
|
1
+
γ
4
| and using it into Lemma 10, taking into account that (
4
+
14
u
)
n
u
(
1
+
γ
5
n
-
1
)
≤
(
5
n
-
1
)
u
(
1
+
γ
5
n
-
1
)
=
γ
5
n
-
1
(
n
≥
2
) and γ
5
n
-
3
γ
5
n
-
1
≤
γ
5
n
-
2
2, we have
(59)
|
∑
j
=
0
n
-
1
‍
s
j
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
‍
s
^
j
⊗
p
j
(
x
)
|
≤
γ
5
n
-
2
2
∑
j
=
0
n
-
1
‍
|
a
j
|
p
~
j
(
|
x
|
)
.
Combining (54) and (59), we get
(60)
|
ϵ
b
o
-
ϵ
b
^
0
|
=
|
∑
j
=
0
n
-
1
s
j
p
j
(
x
)
‍
+
∑
j
=
0
n
-
1
g
j
(
l
)
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
‍
s
^
j
⊗
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
g
^
j
(
l
)
⊗
p
j
(
x
)
|
≤
|
∑
j
=
0
n
-
1
s
j
p
j
(
x
)
‍
-
⊕
∑
j
=
0
n
-
1
s
^
j
⊗
p
j
(
x
)
‍
|
+
|
∑
j
=
0
n
-
1
g
j
(
l
)
p
j
(
x
)
-
⊕
∑
j
=
0
n
-
1
g
^
j
(
l
)
⊗
p
j
(
x
)
‍
|
≤
2
γ
5
n
+
2
2
∑
j
=
0
n
-
1
|
a
j
|
p
~
j
‍
(
|
x
|
)
.From Definitions4 and 6 and Theorem 12, we easily get the following corollary.Corollary 13.
Letp
(
x
)
=
∑
j
=
0
n
‍
a
j
p
j
(
x
) be a Legendre series of degree n and x a floating point value. The relative error bound of the compensated Clenshaw algorithm is
(61)
|
C
o
m
p
C
l
e
n
s
h
a
w
(
p
,
x
)
-
p
(
x
)
|
|
p
(
x
)
|
≤
u
+
2
γ
5
n
+
2
2
c
o
n
d
r
e
l
(
p
,
x
)
.


## 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 entryx. 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 basisp
(
x
)
=
∑
i
=
0
17
‍
a
i
p
i
(
x
) converted by the polynomials P
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 in the neighborhood of its multiple roots 0.75 and 1. We use the C
l
e
n
s
h
a
w, C
o
m
p
C
l
e
n
s
h
a
w 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 C
o
n
v
e
r
t 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 D
D
C
o
n
v
e
r
t 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 C
o
n
v
e
r
t algorithm and D
D
C
o
n
v
e
r
t 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
].Algorithm 7:The conversion algorithm from power basis to Legendre basis.
functionb
=
C
o
n
v
e
r
t
(
a
)
a
,
b are the vectors of monomial polynomial and orthogonal polynomial’s coefficients, respectively
b
(
n
)
=
0,
b
0
(
n
)
=
a
n,
forj
=
n
-
1
:
-
1
:
0
fork
=
0
:
1
:
n
-
j
b
0
(
j
)
=
a
j
-
C
1
A
1
b
1
(
j
+
1
)
(
k
=
0
)
b
k
(
j
)
=
1
A
k
-
1
b
k
-
1
(
j
+
1
)
-
C
k
+
1
A
k
+
1
b
k
+
1
(
j
+
1
)  (
k
=
1
,
…
,
n
-
j
-
2
)
b
k
(
j
)
=
1
A
k
-
1
b
k
-
1
(
j
+
1
)  (
k
≥
n
-
j
-
1
)
end
end
b
=
b
(
0
)From Figures1-2 we observe that the polynomial evaluated by C
l
e
n
s
h
a
w 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 D
D
C
o
n
v
e
r
t algorithm are smaller than those by the C
o
n
v
e
r
t 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 D
D
C
o
n
v
e
r
t algorithm to get the coefficients; the result obtained by the C
o
m
p
C
l
e
n
s
h
a
w algorithm is almost the same as that by using the Symbolic Toolbox in working precision.Figure 1
Evaluation of polynomial converted fromp
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 by the C
o
n
v
e
r
t algorithm in Legendre basis in the neighborhood of its multiple roots, using the Clenshaw algorithm (up), the C
o
m
p
C
l
e
n
s
h
a
w algorithm (middle), and Symbolic Toolbox (down).Figure 2
Evaluation of polynomial converted fromp
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 by the D
D
C
o
n
v
e
r
t algorithm in Legendre basis in the neighborhood of its multiple roots, using the Clenshaw algorithm (up), the C
o
m
p
C
l
e
n
s
h
a
w 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,x
=
0.75
-
1.0
3
2
i
-
85, for i
= 1 : 40 and x
=
0.75
-
1.1
3
i
-
85 for i
= 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 C
l
e
n
s
h
a
w, C
o
m
p
C
l
e
n
s
h
a
w, and D
D
C
l
e
n
s
h
a
w algorithms in Appendix B (see Algorithm 15) and the Symbolic Toolbox, respectively, so that the relative forward errors can be obtained by |
p
res
(
x
)
-
p
sym
(
x
)
|
/
p
sym
(
x
) 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 u (u
≈
1.16
×
1
0
-
16) when the condition number is less than 1
0
17. And the accuracy of both algorithms is decreasing linearly for the condition number larger than 1
0
17. However, the C
l
e
n
s
h
a
w algorithm cannot yield the working precision; the accuracy of which decreases linearly since the condition number is less than 1
0
17. When the condition number is lager than 1
0
17, the Clenshaw algorithm cannot obtain even one significant bit.Figure 3
Accuracy of evaluation of polynomial converted fromp
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 by the D
D
C
o
n
v
e
r
t algorithm in Legendre basis with respect to the condition number.


### 4.3. Time Performances

We can easily know that theT
w
o
S
u
m, T
w
o
P
r
o
d, and F
a
s
t
T
w
o
S
u
m algorithms in Appendix B (Algorithm 8) require 6, 17, and 3 flops, respectively. Then we obtain the computational cost of the C
l
e
n
s
h
a
w, C
o
m
p
C
l
e
n
s
h
a
w, and D
D
C
l
e
n
s
h
a
w algorithms:(i)
C
l
e
n
s
h
a
w: 5n-2 flops;
(ii)
C
o
m
p
C
l
e
n
s
h
a
w: 79n-29 flops;
(iii)
D
D
C
l
e
n
s
h
a
w: 110n-44 flops.Algorithm 8:EFT of the sum of two floating point numbers (
|
a
|
≥
|
b
|
).
function[
x
,
y
]
=
F
a
s
t
T
w
o
S
u
m
(
a
,
b
)
x
=
a
⊕
b
y
=
(
a
⊖
x
)
⊕
bAlgorithm 9:Addition of a double-double number and a double number.
function[
r
h
,
r
l
]
=
a
d
d
_
d
d
_
d
(
a
h
,
a
l
,
b
)
[
t
h
,
t
l
]
=
T
w
o
S
u
m
(
a
h
,
b
)
t
l
=
a
l
⊕
t
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 10:Addition of a double-double number and a double-double number.
function[
r
h
,
r
l
]
=
a
d
d
_
d
d
_
d
d
(
a
h
,
a
l
,
b
h
,
b
l
)
[
s
h
,
s
l
]
=
T
w
o
S
u
m
(
a
h
,
b
h
)
[
t
h
,
t
l
]
=
T
w
o
S
u
m
(
a
l
,
b
l
)
s
l
=
s
l
⊕
t
h
t
h
=
s
h
⊕
s
l
s
l
=
s
l
⊖
(
t
h
⊖
s
h
)
t
l
=
t
l
⊕
s
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 11:Multiplication of a double-double number by a double number.
function[
r
h
,
r
l
]
=
p
r
o
d
_
d
d
_
d
(
a
h
,
a
l
,
b
)
[
t
h
,
t
l
]
=
T
w
o
P
r
o
d
(
a
h
,
b
)
t
l
=
a
l
⊗
b
⊕
t
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 12:Multiplication of a double-double number by a double-double number.
function[
r
h
,
r
l
]
=
p
r
o
d
_
d
d
_
d
d
(
a
h
,
a
l
,
b
h
,
b
l
)
[
t
h
,
t
l
]
=
T
w
o
P
r
o
d
(
a
h
,
b
h
)
t
l
=
(
a
h
⊗
b
l
)
⊕
(
a
l
⊗
b
h
)
⊕
t
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 13:Division of a double number by a double number.
function[
r
h
,
r
l
]
=
d
i
v
_
d
_
d
(
a
,
b
)
q
h
=
a
⊘
b
[
t
h
1
,
t
l
1
]
=
T
w
o
P
r
o
d
(
q
h
,
b
)
[
t
h
2
,
t
l
2
]
=
T
w
o
S
u
m
(
a
,
-
t
h
1
)
q
l
=
(
t
h
2
⊕
t
l
2
⊖
t
l
1
)
⊘
b
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
q
h
,
q
l
)Considering the previous comparison of the accuracy, we observe thatC
o
m
p
C
l
e
n
s
h
a
w is as accurate as D
D
C
l
e
n
s
h
a
w in double precision, but it only needs about 71.8% of flops counting on average. We also implement C
o
m
p
C
l
e
n
s
h
a
w and D
D
C
l
e
n
s
h
a
w by using Microsoft Visual C++ 2008 on Windows 7. Similar to the statement in [23], we assume that the computing time of these algorithms does not depend on the coefficients of polynomial in Legendre basis nor the argument x. So we generate the tested polynomials with random coefficients and arguments in the interval (
-
1,1
), whose degrees vary from 20 to 10000 by the step 50. The average measured computing time ratio of C
o
m
p
C
l
e
n
s
h
a
w to D
D
C
l
e
n
s
h
a
w in C code is 58.29
%. The reason why the measured computing time ratio is better than the theoretical flop count one can be referred to the analysis in terms of instruction level parallelism (ILP) described in [24, 25].


## 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 basisp
(
x
)
=
∑
i
=
0
17
‍
a
i
p
i
(
x
) converted by the polynomials P
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 in the neighborhood of its multiple roots 0.75 and 1. We use the C
l
e
n
s
h
a
w, C
o
m
p
C
l
e
n
s
h
a
w 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 C
o
n
v
e
r
t 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 D
D
C
o
n
v
e
r
t 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 C
o
n
v
e
r
t algorithm and D
D
C
o
n
v
e
r
t 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
].Algorithm 7:The conversion algorithm from power basis to Legendre basis.
functionb
=
C
o
n
v
e
r
t
(
a
)
a
,
b are the vectors of monomial polynomial and orthogonal polynomial’s coefficients, respectively
b
(
n
)
=
0,
b
0
(
n
)
=
a
n,
forj
=
n
-
1
:
-
1
:
0
fork
=
0
:
1
:
n
-
j
b
0
(
j
)
=
a
j
-
C
1
A
1
b
1
(
j
+
1
)
(
k
=
0
)
b
k
(
j
)
=
1
A
k
-
1
b
k
-
1
(
j
+
1
)
-
C
k
+
1
A
k
+
1
b
k
+
1
(
j
+
1
)  (
k
=
1
,
…
,
n
-
j
-
2
)
b
k
(
j
)
=
1
A
k
-
1
b
k
-
1
(
j
+
1
)  (
k
≥
n
-
j
-
1
)
end
end
b
=
b
(
0
)From Figures1-2 we observe that the polynomial evaluated by C
l
e
n
s
h
a
w 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 D
D
C
o
n
v
e
r
t algorithm are smaller than those by the C
o
n
v
e
r
t 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 D
D
C
o
n
v
e
r
t algorithm to get the coefficients; the result obtained by the C
o
m
p
C
l
e
n
s
h
a
w algorithm is almost the same as that by using the Symbolic Toolbox in working precision.Figure 1
Evaluation of polynomial converted fromp
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 by the C
o
n
v
e
r
t algorithm in Legendre basis in the neighborhood of its multiple roots, using the Clenshaw algorithm (up), the C
o
m
p
C
l
e
n
s
h
a
w algorithm (middle), and Symbolic Toolbox (down).Figure 2
Evaluation of polynomial converted fromp
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 by the D
D
C
o
n
v
e
r
t algorithm in Legendre basis in the neighborhood of its multiple roots, using the Clenshaw algorithm (up), the C
o
m
p
C
l
e
n
s
h
a
w 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,x
=
0.75
-
1.0
3
2
i
-
85, for i
= 1 : 40 and x
=
0.75
-
1.1
3
i
-
85 for i
= 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 C
l
e
n
s
h
a
w, C
o
m
p
C
l
e
n
s
h
a
w, and D
D
C
l
e
n
s
h
a
w algorithms in Appendix B (see Algorithm 15) and the Symbolic Toolbox, respectively, so that the relative forward errors can be obtained by |
p
res
(
x
)
-
p
sym
(
x
)
|
/
p
sym
(
x
) 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 u (u
≈
1.16
×
1
0
-
16) when the condition number is less than 1
0
17. And the accuracy of both algorithms is decreasing linearly for the condition number larger than 1
0
17. However, the C
l
e
n
s
h
a
w algorithm cannot yield the working precision; the accuracy of which decreases linearly since the condition number is less than 1
0
17. When the condition number is lager than 1
0
17, the Clenshaw algorithm cannot obtain even one significant bit.Figure 3
Accuracy of evaluation of polynomial converted fromp
(
x
)
=
(
x
-
0.75
)
7
(
x
-
1
)
11 by the D
D
C
o
n
v
e
r
t algorithm in Legendre basis with respect to the condition number.


## 4.3. Time Performances

We can easily know that theT
w
o
S
u
m, T
w
o
P
r
o
d, and F
a
s
t
T
w
o
S
u
m algorithms in Appendix B (Algorithm 8) require 6, 17, and 3 flops, respectively. Then we obtain the computational cost of the C
l
e
n
s
h
a
w, C
o
m
p
C
l
e
n
s
h
a
w, and D
D
C
l
e
n
s
h
a
w algorithms:(i)
C
l
e
n
s
h
a
w: 5n-2 flops;
(ii)
C
o
m
p
C
l
e
n
s
h
a
w: 79n-29 flops;
(iii)
D
D
C
l
e
n
s
h
a
w: 110n-44 flops.Algorithm 8:EFT of the sum of two floating point numbers (
|
a
|
≥
|
b
|
).
function[
x
,
y
]
=
F
a
s
t
T
w
o
S
u
m
(
a
,
b
)
x
=
a
⊕
b
y
=
(
a
⊖
x
)
⊕
bAlgorithm 9:Addition of a double-double number and a double number.
function[
r
h
,
r
l
]
=
a
d
d
_
d
d
_
d
(
a
h
,
a
l
,
b
)
[
t
h
,
t
l
]
=
T
w
o
S
u
m
(
a
h
,
b
)
t
l
=
a
l
⊕
t
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 10:Addition of a double-double number and a double-double number.
function[
r
h
,
r
l
]
=
a
d
d
_
d
d
_
d
d
(
a
h
,
a
l
,
b
h
,
b
l
)
[
s
h
,
s
l
]
=
T
w
o
S
u
m
(
a
h
,
b
h
)
[
t
h
,
t
l
]
=
T
w
o
S
u
m
(
a
l
,
b
l
)
s
l
=
s
l
⊕
t
h
t
h
=
s
h
⊕
s
l
s
l
=
s
l
⊖
(
t
h
⊖
s
h
)
t
l
=
t
l
⊕
s
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 11:Multiplication of a double-double number by a double number.
function[
r
h
,
r
l
]
=
p
r
o
d
_
d
d
_
d
(
a
h
,
a
l
,
b
)
[
t
h
,
t
l
]
=
T
w
o
P
r
o
d
(
a
h
,
b
)
t
l
=
a
l
⊗
b
⊕
t
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 12:Multiplication of a double-double number by a double-double number.
function[
r
h
,
r
l
]
=
p
r
o
d
_
d
d
_
d
d
(
a
h
,
a
l
,
b
h
,
b
l
)
[
t
h
,
t
l
]
=
T
w
o
P
r
o
d
(
a
h
,
b
h
)
t
l
=
(
a
h
⊗
b
l
)
⊕
(
a
l
⊗
b
h
)
⊕
t
l
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
t
h
,
t
l
)Algorithm 13:Division of a double number by a double number.
function[
r
h
,
r
l
]
=
d
i
v
_
d
_
d
(
a
,
b
)
q
h
=
a
⊘
b
[
t
h
1
,
t
l
1
]
=
T
w
o
P
r
o
d
(
q
h
,
b
)
[
t
h
2
,
t
l
2
]
=
T
w
o
S
u
m
(
a
,
-
t
h
1
)
q
l
=
(
t
h
2
⊕
t
l
2
⊖
t
l
1
)
⊘
b
[
r
h
,
r
l
]
=
F
a
s
t
T
w
o
S
u
m
(
q
h
,
q
l
)Considering the previous comparison of the accuracy, we observe thatC
o
m
p
C
l
e
n
s
h
a
w is as accurate as D
D
C
l
e
n
s
h
a
w in double precision, but it only needs about 71.8% of flops counting on average. We also implement C
o
m
p
C
l
e
n
s
h
a
w and D
D
C
l
e
n
s
h
a
w by using Microsoft Visual C++ 2008 on Windows 7. Similar to the statement in [23], we assume that the computing time of these algorithms does not depend on the coefficients of polynomial in Legendre basis nor the argument x. So we generate the tested polynomials with random coefficients and arguments in the interval (
-
1,1
), whose degrees vary from 20 to 10000 by the step 50. The average measured computing time ratio of C
o
m
p
C
l
e
n
s
h
a
w to D
D
C
l
e
n
s
h
a
w in C code is 58.29
%. The reason why the measured computing time ratio is better than the theoretical flop count one can be referred to the analysis in terms of instruction level parallelism (ILP) described in [24, 25].


## 5. Conclusions

This paper introduces a compensated Clenshaw algorithm for accurate evaluation of the finite Legendre series. TheC
l
e
n
s
h
a
w algorithm is not precise enough for an ill-conditioned problem, especially evaluating a polynomial in the neighborhood of a multiple root. However, this new algorithm can yield a full precision accuracy in working precision, as the same as the original Clenshaw algorithm using double-double arithmetic and rounding into the working precision. Meanwhile this compensated Clenshaw algorithm is more efficient which means that it is much more useful to accurately evaluate the polynomials in Legendre basis for ill-conditioned situations.



---

*Source: 742538-2014-07-23.xml*