Keywords

1 Introduction

Tensor product surfaces are bivariate polynomials in tensor product form. In monomial basis, tensor product polynomials are expressed in the following form,

$$p(x,y)=\sum _{i=0}^n\sum _{j=0}^mc_{i,j}x^iy^j.$$

In Computer Aided Geometric Design (CAGD), tensor product surfaces are usually represented in Bézier form [1]

$$p(x,y)=\sum _{i=0}^n\sum _{j=0}^mc_{i,j}B_i^n(x)B_i^m(y), \quad (x,y)\in [0,1]\times [0,1],$$

where \(B_i^k(t)\) is the Bernstein polynomial of degree k as

$$B_i^k(t)=\begin{pmatrix}k\\ i\end{pmatrix}(1-t)^{k-i}t^i, \quad t\in [0,1],\quad i=0,1,\dots ,k.$$

The de Casteljau algorithm [2, 3] is the usual polynomial evaluation algorithm in CAGD. Nevertheless, evaluating a polynomial of degree n, the de Casteljau algorithm needs \(\mathcal {O}(n^2)\) operations, in contrast to the \(\mathcal {O}(n)\) operations of the Volk and Schumaker (VS) algorithm [4]. The VS basis \(z_n:=(z_0^n(t),z_1^n(t),\dots ,z_n^n(t)) (t\in [0,1])\) is given by \(z_i^n(t)=t^i(1-t)^{n-i}\). Otherwise, the VS algorithm consist of Horner algorithm. For evaluating tensor product surfaces, de Casteljau and VS algorithms are more stable and accurate than Horner algorithm [1]. And these three algorithms satisfy the relative accuracy bound

$$\frac{|p(x,y)-\widehat{p}(x,y)|}{|p(x,y)|}\le \mathcal {O}(u)\times {cond}(p,x,y),$$

where \(\widehat{p}(x,y)\) is the computed result, u is the unit roundoff and cond(p, x, y) is the condition number of p(x, y).

From 2005 to 2009, Graillat et al. proposed compensated Horner scheme for univariate polynomials in [5,6,7]. From 2010 to 2013, Jiang et al. presented compensated de Casteljau algorithms to evaluate univariate polynomials and its first order derivative in Bernstein basis in [8], to evaluate bivariate polynomials in Bernstein-Bézier form in [9], and to evaluate Bézier tensor product surfaces in [10]. From 2014 to 2017, Du et al. improved Clenshaw-Smith algorithm [11] for Legendre polynomial series with real number coefficients, bivariate compensated Horner algorithm [12] for tensor product polynomials and the quotient-difference algorithm [13] which is a double nested algorithm. All these algorithms can yield a full precision accuracy in double precision as applying double-double library [14].

This paper presents new compensated VS algorithms, which have less computational cost than compensated de Casteljau algorithm, to evaluate tensor product polynomial surfaces by applying error-free transformations which is exhaustively studied in [15,16,17]. The relative accuracy bound of our proposed compensated algorithms is satisfied

$$\frac{|p(x,y)-\widehat{p}(x,y)|}{|p(x,y)|}\le {u}+\mathcal {O}(u^2)\times {cond}(p,x,y),$$

where \(\widehat{p}(x,y)\) is computed by the compensated algorithms.

The rest of the paper is organized as follows. Section 2 introduces basic notation in error analysis, error-free transformations and condition numbers are also given. Section 3 presents the new compensated VS tensor product algorithm and its error analysis. Finally all the error bounds are compared in numerical experiments in Sect. 4.

2 Preliminary

2.1 Basic Notations

We assume to work with a floating-point arithmetic adhering to IEEE-754 floating-point standard rounding to nearest. In our analysis we assume that there is no computational overflow or underflow. Let \(op\in \{\oplus ,\ominus ,\otimes ,\oslash \}\) represents a floating-point computation, and the evaluation of an expression in floating-point arithmetic is denoted \(fl(\cdot )\), then its computation obeys the model

$$\begin{aligned} a~op~b=(a\circ {b})(1+\varepsilon _1)=(a\circ {b})/(1+\varepsilon _2), \end{aligned}$$
(1)

where \(a,b\in \mathbb {F}\) (the set of floating-point numbers), \(\circ \in \{+,-,\times ,\div \}\) and \(|\varepsilon _1|,|\varepsilon _2|\le {u}\) (u is the round-off unit of the computer). We also assume that if \(a\circ {b}=x\) for \(x\in {\mathbb {R}}\), then the computed result in floating-point arithmetic is denoted by \(\widehat{x}=a~op~b\), and its perturbation is \(\epsilon x\), i.e.

$$\begin{aligned} \widehat{x}=x+\epsilon {x}. \end{aligned}$$
(2)

The following definition and properties will be used in the forward error analysis (see more details in [18]).

Definition 1

We define

$$\begin{aligned} 1+\theta _n=\prod ^n_{i=1}(1+\delta _i)^{\rho _i}, \end{aligned}$$
(3)

where \(|\delta _i|\le {u},\rho _i=\pm 1\) for \(i=1,2,\ldots ,n\), \(|\theta _n|\le {\gamma _n}:=\dfrac{nu}{1-nu}=nu+\mathcal {O}(u^2)\) and \(nu<1\).

Some basic properties in Definition 1 are given by:

  • \(u+\gamma _k\le \gamma _{k+1},\)

  • \(i\gamma _k<\gamma _{ik},\)

  • \(\gamma _k+\gamma _j+\gamma _k\gamma _j\le {\gamma _{k+j}},\)

  • \(\gamma _i\gamma _j\le \gamma _{i+k}\gamma _{j-k},\) if \(0<k<j-i.\)

2.2 Error-Free Transformations

The development of some families of more stable algorithms, which are called compensated algorithms, is based on the paper [15] on error-free transformations (EFT). For a pair of floating-point numbers \(a,b\in \mathbb {F}\), when no underflow occurs, there exists a floating-point number y satisfying \(a \circ b=x+y\), where \(x=\mathrm{fl}(a \circ b)\) and \(\circ {\in }\{+ , - , \times \}\). Then the transformation \((a,b)\longrightarrow (x,y)\) is regarded as an EFT. For division, the corresponding EFT is constructed using the remainder, so its definition is slightly different (see below). The EFT algorithms of the sum, product and division of two floating-point numbers are the TwoSum algorithm [19], the TwoProd algorithm [20] and the DivRem algorithm [21, 22], respectively.

2.3 Condition Numbers

The condition number of polynomials is with respect to the difficulty of the evaluation algorithm. We assume to evaluate a bivariate polynomial p(x, y) in basis \(u\in \mathcal {U}\) at the point (x, y), then for any \((x,y)\in I\), we have

$$\begin{aligned} \begin{aligned} |p(x,y)-\widehat{p}(x,y)|&=|\sum _{i=0}^n\sum _{j=0}^mc_{i,j}u^n_i(x)u^m_i(y)|\\&\le \sum _{i=0}^n\sum _{j=0}^m|c_{i,j}||u^n_i(x)||u^m_i(y)|. \end{aligned} \end{aligned}$$
(4)

We assume that

$$\begin{aligned} \bar{p}(x,y):=\sum _{i=0}^n\sum _{j=0}^m|c_{i,j}||u^n_i(x)||u^m_i(y)|, \end{aligned}$$
(5)

then the relative condition number is

$$\begin{aligned} cond(p,x,y)=\frac{\bar{p}(x,y)}{|p(x,y)|}. \end{aligned}$$
(6)

In [23], it is known that the condition number in VS basis is as same as in Bernstein basis.

3 The Compensated VS Algorithm for Bézior Tensor Product Surfaces

In this section, we show the VS algorithms, including univariate and bivariate ones. We provide a compensated VSTP algorithm for evaluating Bézior tensor product polynomials. Its forward error bound is also given in the end.

3.1 VS Algorithm

The VS algorithm is a nested-type algorithm for the evaluation of bivariate polynomials of total degree n by Schumaker and Volk [4]. Basically, the VS tensor product algorithm could be represented by the univariate VS algorithm.

figure a

Theorem 1 states the forward error bound of VS algorithm.

Theorem 1

[24]. Let \(p(t)=\sum _{i=0}^nc_iz_i^n(t)\) with floating point coefficients \(c_i\) and a floating point value t. Consider the computed result \(\widehat{p}(t)\) with the VS algorithm and its corresponding theoretical result p(t), if \(4nu<1\) where u is the unit roundoff, then

$$\begin{aligned} |p(t)-\widehat{p}(t)|\le \gamma _{4n}\sum ^n_{i=0}|c_{i}z_i^n(t)|. \end{aligned}$$
(7)
figure b

Similar as Theorem 4 in [10], the forward error bound of the VSTP algorithm is easily performed in Theorem 2.

Theorem 2

Let \(p(x,y)=\sum ^n_{i=0}\sum ^m_{j=0}c_{i,j}z_i^n(x)z_i^m(y)\) with floating point coefficients \(c_{i,j}\) and floating point values x, y. Consider the computed result \(\widehat{p}(x,y)\) of the VSTP algorithm and its corresponding theoretical result p(x, y), if \((4n+4m+1)u<1\) where u is the unit roundoff, then

$$\begin{aligned} |p(x,y)-\widehat{p}(x,y)|\le \gamma _{4(n+m)+1}\bar{p}(x,y), \end{aligned}$$
(8)

where \(\bar{p}(x,y)\) is defined in (5) in VS basis.

3.2 The CompVSTP Algorithm

The CompVS algorithm [23] is proposed by Delgado and Peña, which is as accurate as computing in twice the working precision by VS algorithm. In this section, in order to easily provide the forward error bound of CompVS algorithm, we show a compensated Horner algorithm with double-double precision input in Algorithm 3. A compensated power evaluation algorithm in Algorithm 4 is also given.

In Algorithm 3, assuming input x is real number, and we split x into three parts, i.e. \(x=x^{(h)}+x^{(l)}+x^{(m)}, \)where \(x^{(h)},x^{(l)}\in {\mathbb {F}}\), \(x,x^{(m)}\in {\mathbb {R}}\) and \(|x^{(l)}|\le {u|x^{(h)}|},|x^{(m)}|\le {u|x^{(l)}|}\). Since the perturbation of input \(x^{(m)}\) in Algorithm 3 is \(\mathcal {O}(u^2)\), we just need to consider x in double-double precision. According to Theorem 3.1 in [25], the proof of forward error bound of Algorithm 3 in the following theorem is similar as Theorem 12 in [11].

Theorem 3

If \(p(x)=\sum ^n_{i=0}a_ix^i\) \((n\ge 2)\) with floating point coefficients \(a_i\) and a double-double precision number x. And \(\widehat{\epsilon {b}}_0\) is the computed result err of the CompHorner2 algorithm, \(\epsilon {b}_0\) is corresponding theoretical result of \(\widehat{\epsilon {b}}_0\). Then

$$\begin{aligned} |\epsilon {b}_0-\widehat{\epsilon {b}}_0|\le \gamma _{3n-1}\gamma _{3n}\sum ^n_{i=0}|a_i||x^i|. \end{aligned}$$
(9)

Graillat proposes a compensated power evaluation algorithm [26] as follows.

figure c
figure d

Theorem 4

[26]. If \(p(x)=x^n\) \((n\ge 2)\) with a floating-point number x. And \(\widehat{e}\) is the computed result err of the CompLinpower algorithm, e is corresponding theoretical result of \(\widehat{e}\). Then

$$\begin{aligned} |e-\widehat{e}|\le \gamma _{n}\gamma _{2n}|x^n|. \end{aligned}$$
(10)

In [23], Delgado and Peña present the running error analysis of CompVS algorithm, but they do not propose its forward error analysis. Here, combining Algorithms 3 and 4, we show the CompVS algorithm in the following algorithm which is expressed a little different in [23].

figure e

In Algorithm 5, we can easily obtain that \([q^{(h)},q^{(l)}]\) is the double-double form of \(q={(1-x)}/{x}\) if \(x\ge 1/2\) or \(q={x}/{(1-x)}\) if \(x>1/2\). Then, according to Theorems 1, 3 and 4, the forward error bound of CompVS algorithm is proposed in Theorem 5.

Theorem 5

If \(p(t)=\sum ^n_{i=0}c_iz_i^n(t)\) with floating point coefficients \(c_i\) and a floating point value t. And \(\widehat{\epsilon {b}}_0\) is the computed result err of the CompVS algorithm, \(\epsilon {b}_0\) is corresponding theoretical result of \(\widehat{\epsilon {b}}_0\). Then

$$\begin{aligned} |\epsilon {b}_0-\widehat{\epsilon {b}}_0|\le \gamma _{3n+1}\gamma _{3n+2}\sum ^n_{i=0}|c_iz_i^n(t)|. \end{aligned}$$
(11)

Proof

In Algorithm 5, we assume that \(\widehat{f}+e_1=\sum _{i=1}^n p_iq^i\) and \(\widehat{s}+e_2=x^n\). Then, we can obtain that \(p(t)=(\widehat{f}+e1)(\widehat{s}+e2)\) and assume that \(e=e_1\widehat{s}+e_2\widehat{f}+e_1e_2\). Since \(\widehat{e}=\widehat{e}_1\otimes \widehat{s}\oplus \widehat{e}_2\otimes \widehat{f}\), we have

$$\begin{aligned} \begin{aligned}&|e-\widehat{e}|\\&\le |(1+u)^2[(e_1-\widehat{e}_1)\widehat{s}+(e_2-\widehat{e}_2)\widehat{f}+e_1e_2]-(2u+u^2)e|\\&\le {(2u+u^2)}|e|+(1+u)^2(|e_1-\widehat{e}_1||\widehat{s}|+|e_2-\widehat{e}_2||\widehat{f}|). \end{aligned} \end{aligned}$$
(12)

From Theorem 1, let \(\bar{p}(t)=|c_iz_i^n(t)|\), we obtain that

$$\begin{aligned} |e|\le \gamma _{4n}\bar{p}(t). \end{aligned}$$
(13)

Thus

$$\begin{aligned} (2u+u^2)|e|\le \gamma _2\gamma _{4n+1}\bar{p}(t). \end{aligned}$$
(14)

According to Theorem 3, we have

$$\begin{aligned} (1+u)^2|e_1-\widehat{e}_1||\widehat{s}|\le \gamma _{3n}\gamma _{3n+1}\bar{p}(x)+\mathcal {O}(u^2). \end{aligned}$$
(15)

According to Theorem 4, we have

$$\begin{aligned} (1+u)^2|e_2-\widehat{e}_2||\widehat{f}|\le \gamma _{n+1}\gamma _{2n+1}\bar{p}(x)+\mathcal {O}(u^2). \end{aligned}$$
(16)

From (14), (15) and (16), we can deduce (11).

In fact, \(p(x)=\widehat{p}(x)+\epsilon {b}_0\), where \(\epsilon {b}_0\) is corresponding theoretical error of the computed result \(\widehat{p}(x)\). In order to correct the result by Algorithms 1 and 5 find an approximate value \(\widehat{\epsilon {b}}_0\) of \(\epsilon {b}_0\). Motivated by this principle, we propose to use the CompVS algorithm instead of VS algorithm in Algorithm 2 to improve the accuracy of VSTP algorithm. According to Algorithm 2, we assume that

$$\begin{aligned} b_{i,0}=\widehat{b}_{i,0}+err_{i,0}^{(1)}, \quad 0\le {i}\le {n}, \end{aligned}$$
(17)

where \(err_{i,0}^{(1)}\) is the theoretical error of \(\widehat{b}_{i,0}=\mathtt{VS}(c_{i,:},y)\) and

$$\begin{aligned} b_{i,0}=\sum ^m_{j=0}c_{i,j}z_i^m(y), \end{aligned}$$
(18)

is the exact result for each i. Similarly, we have

$$\begin{aligned} \tilde{a}_0=\widehat{a}_0+err^{(2)}, \end{aligned}$$
(19)

where \(err^{(2)}\) is the theoretical error of \(\widehat{a}_0=\mathtt{VS}(\widehat{b}_{:,0},x)\) and

$$\begin{aligned} \tilde{a}_0=\sum _{i=0}^n\widehat{b}_{i,0}z_i^n(x), \end{aligned}$$
(20)

is the exact result. According to (17)–(20), we can deduce

$$\begin{aligned} \sum _{i=0}^n\sum _{j=0}^mc_{i,j}z_i^n(x)z_i^m(y)=\widehat{a}_0+\sum _{i=0}^nerr_{i,0}^{(1)}z_i^n(x)+err^{(2)}, \end{aligned}$$
(21)

i.e.

$$\begin{aligned} p(x,y)=\widehat{p}(x,y)+\sum _{i=0}^nerr_{i,0}^{(1)}z_i^n(x)+err^{(2)}. \end{aligned}$$
(22)

Using CompVS algorithm, we can easily get the approximation values of \(err_{i,0}^{(1)}\) and \(err^{(2)}\), i.e. \(\widehat{err}_{i,0}^{(1)}\) and \(\widehat{err}^{(2)}\). Thus, we propose the CompVSTP algorithm for evaluating Bézier tensor product polynomials in Algorithm 6.

figure f

From (21) and (22), we assume that \(e_1=\sum _{i=0}^nerr_{i,0}^{(1)}z_i^n(x)\) and \(e_2=err^{(2)}\) so that the real error of the computed result is \(e=e_1+e_2\), i.e. \(p(x,y)=\widehat{p}(x,y)+e\). Firstly, we present the bound of \(|e_1-\widehat{e}_1|\) in Lemma 1.

Lemma 1

From Algorithm 6, we assume that \(e_1=\sum _{i=0}^nerr_{i,0}^{(1)}z_i^n(x)\). Then we have

$$\begin{aligned} |e_1-\widehat{e}_1|\le (\gamma _{3n+1}\gamma _{3n+2}(1+\gamma _{4m})+\gamma _{4n}\gamma _{4m}\bar{p}(x,y), \end{aligned}$$
(23)

where \(\bar{p}(x,y)\) is defined in (5) in VS basis.

Proof

We denote that

$$\begin{aligned} \bar{e}_1=\sum _{i=0}^n\widehat{err}_{i,0}^{(1)}z_i^n(x). \end{aligned}$$
(24)

Hence, we have

$$\begin{aligned} |e_1-\widehat{e}_1|\le |e_1-\bar{e}_1|+|\bar{e}_1-\widehat{e}_1|. \end{aligned}$$
(25)

According to Theorem 5, we have

$$\begin{aligned} |err_{i,0}^{(1)}-\widehat{err}_{i,0}^{(1)}|\le \gamma _{3n+1}\gamma _{3n+2}\sum _{j=0}^m|c_{i,j}z_i^m(y)|, \end{aligned}$$
(26)

thus

$$\begin{aligned} \begin{aligned} |e_1-\bar{e}_1|&=\sum _{i=0}^n|err_{i,0}^{(1)}-\widehat{err}_{i,0}^{(1)}|z_i^n(x)\\&\le \gamma _{3n+1}\gamma _{3n+2}\sum _{i=0}^n\sum _{j=0}^m|c_{i,j}z_i^n(x)z_i^m(y)|. \end{aligned} \end{aligned}$$
(27)

According to Theorem 1, we obtain

$$\begin{aligned} |\bar{e}_1-\widehat{e}_1|\le \gamma _{4m}\sum _{i=0}^n|\widehat{err}_{i,0}^{(1)}z_i^n(x)|. \end{aligned}$$
(28)

Then we have that

$$\begin{aligned} |\widehat{err}_{i,0}^{(1)}|\le |err_{i,0}^{(1)}|+|err_{i,0}^{(1)}-|\widehat{err}_{i,0}^{(1)}|. \end{aligned}$$
(29)

By Theorem 1, we have

$$\begin{aligned} |err_{i,0}^{(1)}|\le \gamma _{4n}\sum _{j=0}^m|c_{i,j}z_i^m(y)|. \end{aligned}$$
(30)

From (26), (29) and (30), we deduce that

$$\begin{aligned} |\widehat{err}_{i,0}^{(1)}|\le (\gamma _{3n+1}\gamma _{3n+2}+\gamma _{4n})\sum _{j=0}^m|c_{i,j}z_i^m(y)|, \end{aligned}$$
(31)

and then from (28) we obtain

$$\begin{aligned} |\bar{e}_1-\widehat{e}_1|\le \gamma _{4m}(\gamma _{3n+1}\gamma _{3n+2}+\gamma _{4n})\bar{p}(x,y). \end{aligned}$$
(32)

Hence, from (25), (27) and (32), we can obtain (23).

Then, we present the bound of \(|e_2-\widehat{e}_2|\) in Lemma 2.

Lemma 2

From Algorithm 6, we assume that \(e_2=err^{(2)}\). Then we have

$$\begin{aligned} |e_2-\widehat{e}_2|\le \gamma _{3m+1}\gamma _{3m+2}(1+\gamma _{4m})\bar{p}(x,y), \end{aligned}$$
(33)

where \(\bar{p}(x,y)\) is defined in (5) in VS basis.

Proof

According to Theorem 5, we have

$$\begin{aligned} |e_2-\widehat{e}_2|\le \gamma _{3m+1}\gamma _{3m+2}\sum _{i=0}^n|\widehat{b}_{i,0}z_i^n(x)|. \end{aligned}$$
(34)

From Theorem 1, we obtain

$$\begin{aligned} |\widehat{b}_{i,0}|\le \sum _{j=0}^m(1+\gamma _{4m})|c_{i,j}z_i^m(y)|. \end{aligned}$$
(35)

Hence, from (34) and (35), we can deduce (33).

Above all, the forward error bound of CompVSTP algorithm is performed in the following theorem.

Theorem 6

Let \(p(x,y)=\sum _{i=0}^n\sum _{j=0}^mc_{i,j}z_i^n(x)z_i^m(y)\) with floating point coefficients \(c_{i,j}\) and floating point values x, y. The forward error bound of Algorithm 6 is

$$\begin{aligned} |CompVSTP(p,x,y)-p(x,y)|\le {u}|p(x,y)|+3(\gamma _{4n+2}^2+\gamma _{4m+2}^2)\bar{p}(x,y), \end{aligned}$$
(36)

where \(\bar{p}(x,y)\) is defined in (5) in VS basis.

Proof

We assume that \(e_1=\sum _{i=0}^nerr_{i,0}^{(1)}x^i\) and \(e_2=err^{(2)}\) so that \(e=e_1+e_2\). From (22), we have

$$\begin{aligned} p(x,y)=\widehat{p}(x,y)+e, \end{aligned}$$
(37)

and from Algorithm 6, we have

$$\begin{aligned} \mathtt{CompVSTP}(p,x,y)=\widehat{p}(x,y)\oplus \widehat{e}. \end{aligned}$$
(38)

Hence

$$\begin{aligned} \begin{aligned} |\mathtt{CompVSTP}(p,x,y)-p(x,y)|&\le |(1+u)(p(x,y)-e+\widehat{e})-p(x,y)|\\&\le {u}|p(x,y)|+(1+u)|e-\widehat{e}|. \end{aligned} \end{aligned}$$
(39)

Since \(\widehat{e}=\widehat{e}_1\oplus \widehat{e}_2\), we have

$$\begin{aligned} \begin{aligned} |e-\widehat{e}|&\le |(1+u)(e_1-\widehat{e}_1+e_2-\widehat{e}_2)-ue|\\&\le {u}|e|+(1+u)(|e_1-\widehat{e}_1|+|e_2-\widehat{e}_2|). \end{aligned} \end{aligned}$$
(40)

From Theorem 2, we obtain that

$$\begin{aligned} |e|\le \gamma _{4(n+m)+1}\bar{p}(x,y). \end{aligned}$$
(41)

Thus

$$\begin{aligned} u(1+u)|e|\le \gamma _1\gamma _{4(n+m+1)}\bar{p}(x,y)\le \gamma _{4n+2}\gamma _{4m+2}\bar{p}(x,y)\le \frac{1}{2}(\gamma _{4n+2}^2+\gamma _{4m+2}^2)\bar{p}(x,y). \end{aligned}$$
(42)

According to Lemma 1, we have

$$\begin{aligned} \begin{aligned} (1+u)^2|e_1-\widehat{e}_1|&\le (2\gamma _{4n+1}^2+\gamma _{4n+1}\gamma _{4m+1})\bar{p}(x,y)\\&\le (\frac{5}{2}\gamma _{4n+1}^2+\frac{1}{2}\gamma _{4m+1}^2)\bar{p}(x,y). \end{aligned} \end{aligned}$$
(43)

According to Lemma 2, we have

$$\begin{aligned} (1+u)^2|e_2-\widehat{e}_2|\le 2\gamma _{4m+1}^2\bar{p}(x,y). \end{aligned}$$
(44)

From (42), (43) and (44), we can deduce (36).

According to the relative condition number defined in (6), we can deduce Corollary 1.

Corollary 1

Let \(p(x,y)=\sum _{i=0}^n\sum _{j=0}^mc_{i,j}z_i^n(x)z_i^m(y)\) with floating point coefficients \(c_{i,j}\) and floating point values x, y. The forward relative error bound of Algorithm 6 is

$$\begin{aligned} \frac{|CompVSTP(p,x,y)-p(x,y)|}{|p(x,y)|}\le {u}+3(\gamma _{4n+2}^2+\gamma _{4m+2}^2)cond(p,x,y). \end{aligned}$$
(45)

4 Numerical Experiments

In this section, we compare CompVSTP algorithm against an implementation of VSTP algorithm that applies the double-double format [14, 27] which we denote as DDVSTP algorithm. In fact, since the working precision is double precision, the double-double arithmetic is the most efficient way to yield a full precision accuracy of evaluating polynomials. Moreover, we also compare CompVSTP algorithm against compensated de Casteljau (CompDCTP) algorithm [10].

Fig. 1.
figure 1

Accuracy of evaluation of ill-conditioned Bézier tensor product polynomials with respect to the condition number

Full size image

All our experiments are performed using IEEE-754 double precision as working precision. All the programs about accuracy measurements have been written in Matlab R2014a on a 1.4-GHz Intel Core i5 Macbook Air. We focus on the relative forward error bounds for ill-conditioned Bézier tensor product polynomials. We use a similar GenPoly algorithm [10, 21] to generate tested polynomials p(x, y). The generated polynomials are \(6\times 7\) degree with condition numbers varying from \(10^4\) to \(10^{36}\), x and y are random numbers in [0, 1] and the inspired computed results of all the tested polynomials are 1. We evaluate the polynomials by the VSTP, CompVSTP, CompDCTP, DDVSTP algorithms and the Symbolic Toolbox, respectively, so that the relative forward errors can be obtained by \((|p_{res}(x,y)-p_{sym}(x,y)|)/|p_{sym}(x,y)|\) and the relative error bounds are described from Corollary 1. Note that the condition number of Bézier tensor product polynomials in Bernstein basis evaluated by CompDCTP algorithm is as same as in VS basis evaluated by CompVSTP algorithm. Then we present the relative forward errors of evaluation of the tested polynomials in Fig. 1. As we can see, the relative errors of CompVSTP, CompDCTP and DDVSTP algorithms are both smaller than u (\(u\approx {1.16\times {10}^{-16}}\)) when the condition number is less than \(10^{16}\). And the accuracy of them is decreasing linearly for the condition number larger than \(10^{16}\). However, the VSTP algorithm can not yield the working precision; the accuracy of which decreases linearly since the condition number is less than \(10^{16}\).

At last, we give the computational cost of VSTP, CompVSTP, CompDCTP and DDVSTP algorithms.

  • VSTP: \((3\text {n}+2)(\text {m}+1)+3\text {m}+2\) flops,

  • CompVSTP: \((50\text {n}+26)(\text {m}+1)+50\text {m}+26+1\) flops,

  • CompDCTP: \((24n^2+24\text {n}+7)(\text {m}+1)+24m^2+24\text {m}+7+1\) flops,

  • DDVSTP: \((68\text {n}+120)(\text {m}+1)+68\text {m}+120\) flops.

CompVSTP and DDVSTP algorithms require almost 17 and 23 times flop than VSTP algorithm, respectively. Meanwhile, CompDCTP algorithm requires \(\mathcal {O}(n^2m)\) flop which is much more than \(\mathcal {O}(nm)\). Hence, CompVSTP algorithm only needs about \(73.5\%\) of flops counting on average of DDVSTP algorithm and needs much less computational cost than CompDCTP algorithm. Meanwhile, CompVSTP algorithm is as accurate as CompDCTP and DDVSTP algorithms.

5 Conclusions and Further Work

In this paper, we present CompVSTP algorithm to evaluate Bézier tensor product polynomials, which are compensated algorithms that obtaining an approximate error to correct the computed results by original algorithm. The proposed algorithm is as accurate as computing in double-double arithmetic which is the most efficient way to yield a full precision accuracy. Moreover, it needs fewer flops than counting on average with double-double arithmetic.

A similar approach can be applied to other problems to obtain compensated algorithms. For example we can consider the evaluation of ill-conditioned tensor product polynomials in orthogonal basis like Chebyshev and Legendre basis. Instead of tensor product surfaces, we can consider triangle surfaces like Bernstein-Bézier form. We can also study compensated algorithms for multivariate polynomials.