---
author:
- '**Stanis[ł]{}aw Lewanowicz**${}^{a,\!}$ [^1] , **Pawe[ł]{} Keller**${}^{b}$, **Pawe[ł]{} Woźny**${}^{a}$'
date: |
    [*${}^{a}\!\!\!$ Institute of Computer Science, University of Wroc[ł]{}aw, ul. Joliot-Curie 15, 50-383 Wroc[ł]{}aw, Poland\
    ${}^{b}\!\!\!$ Faculty of Mathematics and Information Science, Warsaw University of Technology,\
    ul. Koszykowa 75, 00-662 Warszawa, Poland*]{}
title: '**Bézier form of dual bivariate Bernstein polynomials**'
---

[**Abstract** Dual Bernstein polynomials of one or two variables have proved to be very useful in obtaining Bézier form of the $L^2$-solution of the problem of best polynomial approximation of Bézier curve or surface. In this connection, the Bézier coefficients of dual Bernstein polynomials are to be evaluated at a reasonable cost. In this paper, a set of recurrence relations satisfied by the Bézier coefficients of dual bivariate Bernstein polynomials is derived and an efficient algorithm for evaluation of these coefficients is proposed. Applications of this result to some approximation problems of Computer Aided Geometric Design (CAGD) are discussed.]{}\
[**Key words** Dual bivariate Bernstein basis $\cdot$ Bézier coefficients $\cdot$ Bivariate Jacobi polynomials $\cdot$ Bivariate Hahn polynomials]{}\
[**Mathematics Subject Classification (2010)** 41A10 $\cdot$ 41A63 $\cdot$ 65D17 $\cdot$ 33C45]{}

Introduction and preliminaries {#S:intro}
==============================

The Bernstein-Bézier curves and surfaces have become the standard in the CAGD context due to their favourable geometric properties (cf. [@Far02]). Degree reduction, which consists in approximating a Bézier curve or surface by another one of a lower degree, and approximation of a rational Bézier curve or surface by a polynomial one have many important applications in geometric modelling, such as data exchange, data compression, and data comparison. There have been many papers relevant to this problem (see, e.g., [@Ahn03; @ALPY04; @CW02; @Eck95; @Hu13; @LBPY02; @LuW06; @LuW08; @Rab05; @RLY06; @RA06; @Sun05; @SL04; @Sha13; @WL09; @WL10]).

The above approximation problems, with or without constraints, are studied for different choices of the error norm. Notice that the Bernstein polynomials do not form an orthogonal base, so obtaining the Bézier form of the best $L_2$-norm solution requires some effort. In a frequently used approach, the main tool applied was transformation between the Bernstein and orthogonal polynomial bases. Such methods are not only expensive, but also may be ill-conditioned (cf. a remark in [@LuW08]).

Recently [@WL09; @WL10], a novel approach to the best $L_2$ approximation problem, using the dual basis associated with the univariate or bivariate Bernstein basis, was proposed. The new methods do not use basis transformation matrices explicitly, hence they do not share the abovementioned limitation. High efficiency of the methods was obtained thanks to the application of recursive properties of the dual Bernstein polynomials.

In the present paper, we derive the recurrence relations satisfied by the Bézier coefficients of the dual bivariate Bernstein polynomials and prove in detail a low-cost algorithm for numerical evaluation of these coefficients; see Section \[S:EE\]. This algorithm (without proof) was successfully used in [@LKW17] as a part of the method for polynomial approximation of rational triangular Bézier surfaces[^2]; see the discussion in Section \[S:Appl\]. Several useful properties of the dual bivariate Bernstein basis polynomials, exploited in the proposed method, are obtained in Appendix A. In Appendix B, some results on the bivariate Jacobi and Hahn orthogonal polynomials are collected.

Below, we introduce notation and definitions that are used in the paper.

For $\bm y:=(y_1,y_2,\ldots,y_d)\in{\mathbb R}^d$, we denote $$|\bm y|:=y_1+y_2+\ldots+y_d,\qquad  \|\bm y\|:=\left(y^2_1+y^2_2+\ldots+y^2_d\right)^{\frac12}.$$ For $n\in{\mathbb N}$ and $\bm c:=(c_1,\,c_2,\,c_3)\in{\mathbb N}^3$ such that ${\mbox{$\left|{\mbox{$\bm c$}}\right|$}}<n$, we define the following sets (cf. Figure \[fig:Fig-bc\]): $$\label{E:TOG}
    \left.\begin{array}{l}
    \Theta_n:=\{\bm k=(k_1,k_2)\in{\mathbb N}^2:\: 0\le|\bm k|\le n\},  \\[1ex]
    \Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n:=\{\bm k=(k_1,k_2)\in{\mathbb N}^2:\:k_1\ge c_1,\,k_2\ge c_2,\,|\bm k|\le n-c_3\},\\[1ex]
    \Gamma^{{\mbox{\scriptsize$\bm{c}$}}}_n:=\Theta_n\setminus\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n.
      \end{array}\;\right\}$$

  ----------------------------- ------------------ ------------------
                                                   
   \[-0.2ex\] $\bm{c}=(0,0,2)$   $\bm{c}=(0,1,2)$   $\bm{c}=(2,1,3)$
  ----------------------------- ------------------ ------------------

\[R:TG\] In applications to approximation problems related to triangular (rational or polynomial) Bézier patches, the set $\Theta_n$ corresponds to the set of control points of a patch, while its subset $\Gamma_n^{{\mbox{${\mbox{\scriptsize$\bm{c}$}}$}}}$ – to the boundary points, where some constraints are to be imposed. See Section \[S:Appl\].

Let $T$ be the standard triangle in ${\mathbb R}^2$, $$\label{E:T}        T:=\{(x_1,\,x_2)\,:\,x_1,\,x_2\ge0,\: x_1+x_2\le1\}.$$ For $n\in{\mathbb N}$ and $\bm k:=(k_1,k_2)\in\Theta_n$, we denote $$\binom{n}{\bm k}:=\frac{n!}{k_1!k_2!(n-|\bm k|)!}.$$ The *shifted factorial* is defined for any $a\in{\mathbb C}$ by $$(a)_0:=1;\quad (a)_k:=a(a+1)\cdots(a+k-1),\quad   k\ge1.$$ The *Bernstein polynomial basis* in $\Pi^2_n$, $n\in{\mathbb N}$, is given by (see, e.g., [@Far86], or [@Far02 §18.4]) $$\label{E:Ber2}    B^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x):=\binom{n}{\bm k}x_1^{k_1}x_2^{k_2}(1-|\bm x|)^{n-|{\mbox{\scriptsize$\bm{k}$}}|},
    \qquad \bm k:=(k_1,k_2)\in\Theta_n,\quad\bm x:=(x_1,x_2).$$

The (unconstrained) *dual bivariate Bernstein basis polynomials* [@LW06], $$\label{E:dBer2}
    D^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm\cdot;{\mbox{$\bm{\alpha}$}})\in\Pi^2_n,  
    \qquad \bm k\in\Theta_n,$$ are defined so that $$\left\langle D^{n}_{{\mbox{\scriptsize$\bm{k}$}}},\,B^n_{{\mbox{\scriptsize$\bm{l}$}}} \right\rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}=\delta_{{\mbox{\scriptsize$\bm{k}$}}, {\mbox{\scriptsize$\bm{l}$}}},
    \qquad \bm k,\bm l\in\Theta_n.$$ Here, $\delta_{{\mbox{\scriptsize$\bm{k}$}}, {\mbox{\scriptsize$\bm{l}$}}}$ equals 1 if $\bm k= \bm l$, and 0 otherwise, while the inner product is defined by $$\label{E:Jinprod}\langle f, g \rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}} :={\mbox{$\displaystyle\int\!\!\!\!\int\limits_{\!\!\!\!T}$}}w_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x)f(\bm x)\,g(\bm x){\,\mbox{{\rm d}$\bm x$}},$$ where the weight function $w_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}$ is given by $$\label{E:w}    w_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x):=A_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}x_1^{\alpha_1}x_2^{\alpha_2}(1-|\bm x|)^{\alpha_3}, 
    \qquad{\mbox{$\bm{\alpha}$}}:=(\alpha_1,\,\alpha_2,\,\alpha_3), \quad\alpha_i>-1,$$ and $A_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}$ is a normalisation factor (see ).

In the sequel, we use the notation $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$ for the *connection coefficients* between bivariate Bernstein and dual bivariate Bernstein bases such that $$\label{E:D2inB2-in}    D^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})=\,
    \sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Theta_n}\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)\,
    B^n_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x),\qquad \bm k\in\Theta_n.$$ Investigating the properties of these coefficients is one of the main goals of the paper. The fast recursive scheme for computing the coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$, $\bm k,\bm l\in\Theta_n$, is formulated in the next section, while the wide mathematical background is given in Section \[SS:dualB2\].

For $n\in{\mathbb N}$ and $\bm c:=(c_1,\,c_2,\,c_3)\in{\mathbb N}^3$ such that ${\mbox{$\left|{\mbox{$\bm c$}}\right|$}}<n$, define the constrained bivariate polynomial space $$\Pi^{{\mbox{\scriptsize$\bm{c}$}},\,2}_n :=\left\{P\in\Pi^2_n\::\: 
    P(\bm x)= x_1^{c_1}x_2^{c_2}(1-|\bm x|)^{c_3}\cdot\,Q(\bm x),
\;\mbox{where}\;Q\in\Pi^2_{n-|{\mbox{\scriptsize$\bm{c}$}}|}\right\}.$$ It can be easily seen that the constrained set of polynomials $B^n_{{\mbox{\scriptsize$\bm{k}$}}}$ with the range of $\bm k$ restricted to $\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n$ forms a basis in this space. We define the *constrained dual bivariate Bernstein basis polynomials*, $$\label{E:constrdBer2}
    D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}}(\bm{\cdot};{\mbox{$\bm{\alpha}$}})\in\Pi^{{\mbox{\scriptsize$\bm{c}$}},\,2}_n,
    \qquad \bm k\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n,$$ so that $$\left\langle D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}},\,B^n_{{\mbox{\scriptsize$\bm{l}$}}}
    \right\rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}=\delta_{{\mbox{\scriptsize$\bm{k}$}}, {\mbox{\scriptsize$\bm{l}$}}} \quad \mbox{for}
    \quad \bm k,\bm l\in\Omega^{{\mbox{${\mbox{\scriptsize$\bm{c}$}}$}}}_n,$$ where the notation of is used. For $\bm c=(0,0,0)$, basis reduces to the unconstrained basis in $\Pi^2_n$. Obviously, for any ${\mbox{$\mathsf Q$}}\in\Pi^{{\mbox{\scriptsize$\bm{c}$}},\,2}_n$, we have $$\label{E:Bezform}
    {\mbox{$\mathsf Q$}}(\bm x)=\sum_{{\mbox{\scriptsize$\bm{k}$}}\in\Omega^{{\mbox{${\mbox{\scriptsize$\bm{c}$}}$}}}_n}\left\langle {\mbox{$\mathsf Q$}}, D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}}\right\rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}\,B^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x).$$

Let $E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},n)$, $\bm l\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n$, denote the Bézier coefficients of the constrained bivariate dual Bernstein polynomial $D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})$, $\bm k\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n$: $$\label{E:Dc2inB2-in}    D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})=\sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n}
    \,E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},n)\,B^n_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x).$$ According to Lemma \[L:Dc2inB2\], the following formula relates the above coefficients with the ”unconstrained” connection coefficients defined in : $$\label{E:Dc2inB2-coef-impl}    E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},n):=U\,V_{{\mbox{\scriptsize$\bm{k}$}}}(n)\,V_{{\mbox{\scriptsize$\bm{l}$}}}(n)\,
    \,e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{c}$}}'}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{c}$}}'}({\mbox{$\bm{\alpha}$}}+2\bm c,n-|\bm c|),\qquad {\bm k}, {\bm l}\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n,$$ where $\bm c':=(c_1,c_2)$, and $U,V_{{\mbox{\scriptsize$\bm{k}$}}}(n)$, $\bm k\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n$, are constants defined in .

Hence, for fast evaluation of the coefficients , it would be sufficient to have an efficient algorithm for computing the quantities $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$ for the full range ${\bm k}, {\bm l}\in\Theta_n$, with arbitrary parameters ${\mbox{$\bm{\alpha}$}}$ and $n$. In Section \[S:EE\], we describe such an algorithm that has the computational complexity order proportional to the number of these quantities.

Computing the Bézier coefficients $e^{\bm k}_{\bm l}({\mbox{$\bm{\alpha}$}},n)$ {#S:EE}
===============================================================================

The coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}\equiv e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$, ${\bm k}, {\bm l}\in\Theta_n$, can be arranged in a block triangular matrix $$\label{E:E-tab} 
    \mathbb E=\left[\,\mathbb E^{{\mbox{\scriptsize$\bm{k}$}}}\,\right]_{{\mbox{\scriptsize$\bm{k}$}}\in\Theta_n}=\left[\begin{array}{lllll}
    \mathbb E^{(0,n)}\\[0.75ex]
    \mathbb E^{(0,n-1)}&\mathbb E^{(1,n-1)}\\[0.75ex]
\multicolumn{3}{c}{\dotfill}\\[0.75ex]
\mathbb E^{(0,1)}&\mathbb E^{(1,1)}&\ldots&\mathbb E^{(n-1,1)}\\[0.75ex]
\mathbb E^{(0,0)}&\mathbb E^{(1,0)}&\ldots&\mathbb E^{(n-1,0)}&\mathbb E^{(n,0)}
\end{array}\right],$$ each block $\mathbb E^{{\mbox{\scriptsize$\bm{k}$}}}$, $\bm k\in\Theta_n$, being a triangular matrix $$\mathbb E^{{\mbox{\scriptsize$\bm{k}$}}}=\left[\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}\,\right]_{{\mbox{\scriptsize$\bm{l}$}}\in\Theta_n}=
\left[\begin{array}{lllll}
e^{{\mbox{\scriptsize$\bm{k}$}}}_{(0,n)}\\[1ex]
e^{{\mbox{\scriptsize$\bm{k}$}}}_{(0,n-1)}&e^{{\mbox{\scriptsize$\bm{k}$}}}_{(1,n-1)}\\[1ex]
\multicolumn{3}{c}{\dotfill}\\[1ex]
e^{{\mbox{\scriptsize$\bm{k}$}}}_{(0,1)}&e^{{\mbox{\scriptsize$\bm{k}$}}}_{(1,1)}&\ldots&e^{{\mbox{\scriptsize$\bm{k}$}}}_{(n-1,1)}\\[1ex]
e^{{\mbox{\scriptsize$\bm{k}$}}}_{(0,0)}&e^{{\mbox{\scriptsize$\bm{k}$}}}_{(1,0)}&\ldots&e^{{\mbox{\scriptsize$\bm{k}$}}}_{(n-1,0)}&e^{{\mbox{\scriptsize$\bm{k}$}}}_{(n,0)}
\end{array}\right].$$

The proposed algorithm is based on some recurrences satisfied by the elements of the matrix . In the sequel, we assume that $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}=0$ if $\bm k\not\in\Theta_n$ or $\bm l\not\in\Theta_n$.

The *first recurrence* is the following (cf. Lemma \[L:e-recur\]): $$\label{E:firstrec}    e^{{\mbox{\scriptsize$\bm{k}$}}+{\mbox{\scriptsize$\bm{v}$}}_2}_{{\mbox{\scriptsize$\bm{l}$}}}=\left([\sigma_1(\bm k)-\sigma_1(\bm l)]e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
    -\sigma_2(\bm k) e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{v}$}}_2}_{{\mbox{\scriptsize$\bm{l}$}}} 
    +\sigma_0(\bm l) e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}+{\mbox{\scriptsize$\bm{v}$}}_2}
    +\sigma_2(\bm l) e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{v}$}}_2}\right)/\sigma_0(\bm k),$$ where $\bm v_2:=(0,1)$, and where for $\bm t:=(t_1,t_2)$, we define $$\sigma_0(\bm t):=(|\bm t|-n)(t_2+\alpha_2+1),\quad  
    \sigma_2(\bm t):=t_2(|\bm t|-\alpha_3-n-1),\quad
    \sigma_1(\bm t):=\sigma_0(\bm t)+\sigma_2(\bm t).$$ Observe that recurrence relates three consecutive blocks of a column of the matrix , shown in the following diagram: $$\begin{array}{l}
    \mathbb E^{{\mbox{\scriptsize$\bm{k}$}}+{\mbox{\scriptsize$\bm{v}$}}_2}\\[0.75ex]
        \mathbb E^{{\mbox{\scriptsize$\bm{k}$}}}\\[0.75ex]
    \mathbb E^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{v}$}}_2}\,.
\end{array}$$ According to the convention, the block $\mathbb E^{(k_1,-1)}$ has only zero elements. Thus, we can compute all the blocks of this column, provided that we have computed the block $\mathbb E^{(k_1,0)}$ using another method.

Now, to *initialise* the computation of the first column, we need a method to compute the corner block $\mathbb  E^{(0,0)}$ in the table . We can use the following formula (cf. Corollary \[C:e00\]): $$\label{E:E00}    e^{(0,0)}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n) =\frac{(-1)^{l_1}(|{\mbox{$\bm{\alpha}$}}|+3)_{n}}{n!(\alpha_1+2)_{l_1}}
    \sum_{i=0}^{n-l_1}C^\ast_i\,h_i(l_2;\alpha_2,\alpha_3,n-l_1),
    \qquad \bm l:=(l_1,l_2)\in\Theta_n,$$ where the coefficients $C^\ast_i$ are given by , and we use the notation $h_i(t;a,b,M)$ for the univariate Hahn polynomials (cf. ). As noticed in Remark \[R:Clenshaw\], we can efficiently evaluate the sum in using the Clenshaw’s algorithm, at the cost of $O(n-l_1)$ operations.

To compute the remaining part of the last row of the matrix , i.e., $$\mathbb E^{(1,0)}\quad\mathbb E^{(2,0)}\quad\ldots\quad\mathbb E^{(n-1,0)}\quad\mathbb E^{(n,0)},$$ we use the *second recurrence* (cf. Lemma \[L:e-recur\]): $$\label{E:secrec}    e^{{\mbox{\scriptsize$\bm{k}$}}+{\mbox{\scriptsize$\bm{v}$}}_1}_{{\mbox{\scriptsize$\bm{l}$}}}
    =\left([\tau_1(\bm k)- \tau_1(\bm l)]e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
    -\tau_2(\bm k)\,e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{v}$}}_1}_{{\mbox{\scriptsize$\bm{l}$}}}
                 +
                 \tau_0(\bm l)\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}+{\mbox{\scriptsize$\bm{v}$}}_1}                  
                 + \tau_2(\bm l)\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{v}$}}_1}\right)/\tau_0(\bm k),$$ where $\bm v_1:=(1,0)$, and for $\bm t:=(t_1,t_2)$, the coefficients $\tau_{j}(\bm t)$ are given by $$\tau_0(\bm t):=(|\bm t|-n)(t_1+\alpha_1+1),\quad    
    \tau_2(\bm t):=t_1(|\bm t|-\alpha_3-n-1), \quad 
    \tau_1(\bm t):=\tau_0(\bm t)+\tau_2(\bm t).$$ The recurrence relates each three consecutive blocks of a row of the matrix , shown in the following diagram: $$\begin{array}{lll}
\mathbb E^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{v}$}}_1}&\quad\mathbb E^{{\mbox{\scriptsize$\bm{k}$}}}&\quad\mathbb E^{{\mbox{\scriptsize$\bm{k}$}}+{\mbox{\scriptsize$\bm{v}$}}_1}\,.
\end{array}$$ Again, by the convention, the blocks $\mathbb E^{(-1,k_2)}$ have only zero elements.

The following algorithm can be applied to compute the complete set of the coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}$ arranged in the block matrix $\mathbb E$ (cf. ). The algorithm requires $O(n^4)$ operations, i.e., its cost is proportional to the total number of the coefficients.

\[A:E-comp\]

[Step 1]{}

:   For $l_1=0,1,\ldots, n-1$,\
    $l_2=0,1,\ldots,n-l_1$,\
    using the Clenshaw algorithm, compute $e^{(0,0)}_{(l_1,l_2)}$ defined by , ,\
    and put $e^{(l_1,l_2)}_{(0,0)}:=e^{(0,0)}_{(l_1,l_2)}$.

[Step 2]{}

:   For $k_1=0,1,\ldots,n-1$,

    $1^o$

    :   for $k_2=0,1,\ldots,n-k_1-1$,\
        $l_1=k_1,k_1+1,\ldots,n$,\
        $l_2=0,1,\ldots,n-l_1$,\
        compute $e^{(k_1,k_2+1)}_{(l_1,l_2)}$ using the recurrence , and put $e^{(l_1,l_2)}_{(k_1,k_2+1)}:=e^{(k_1,k_2+1)}_{(l_1,l_2)}$;

    $2^o$

    :   for $l_1=k_1+1,k_1+2,\ldots,n$,\
        $l_2=0,1,\ldots,n-l_1$,\
        compute $e^{(k_1+1,0)}_{(l_1,l_2)}$ using the recurrence , and put $e^{(l_1,l_2)}_{(k_1+1,0)}:=e^{(k_1+1,0)}_{(l_1,l_2)}$.

**Output**: The set of the coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}}, n)$ for $\bm k,\bm l\in\Theta_n$.

\[R:sym\]

(i)

:   In Algorithm \[A:E-comp\], we made use of the symmetry property $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{{\mbox{\scriptsize$\bm{l}$}}}_{{\mbox{\scriptsize$\bm{k}$}}}({\mbox{$\bm{\alpha}$}},n)$ (cf. ).

(ii)

:   When $\alpha_2=\alpha_3$, the cost of completing the table $\mathbb E$ can be reduced significantly using . First, the lower part of the table , containing the blocks $$\left.\begin{array}{l}
        \mathbb E^{(k_1,\lfloor\frac12(n-k_1)\rfloor)} \\[1ex]
        \vdots\\[1ex]
        \mathbb E^{(k_1,1)}\\[1ex]
        \mathbb E^{(k_1,0)}
        \end{array}\;\right\}\;\mbox{with}\; k_1=0,1,\ldots,n,$$ is computed using the properly adapted Algorithm \[A:E-comp\]. In the second step, the upper part of the table is computed, using the formula $
            e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{\hat{{\mbox{\scriptsize$\bm{k}$}}}}_{\hat{{\mbox{\scriptsize$\bm{l}$}}}}({\mbox{$\bm{\alpha}$}},n),
            $ where $\hat{\bm k}:=(k_1,n-|\bm k|)$ and $\hat{\bm l}:=(l_1,n-|\bm l|)$.

(iii)

:   Similar effect of the cost reduction can be obtained if $\alpha_1=\alpha_2$, or $\alpha_1=\alpha_3$, or $\alpha_1=\alpha_2=\alpha_3$ (cf. –).

(iv)

:   Observe that the complexity order of Algorithm \[A:E-comp\] equals $O(n^4)$, i.e., is proportional to the total number of the coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}}, c,n)$ for $\bm k,\bm l\in\Theta_n$.

Applications {#S:Appl}
============

Rational Bézier curves and surfaces, being the most natural generalization of polynomial Bézier curves and surfaces, are important tools in geometric modelling. However, they are sometimes inconvenient in practical applications. For this reason, several algorithms for approximating a rational Bézier geometric form by a polynomial one have been proposed (see, e.g., [@Hu13; @LWK12; @Rab05; @Sha13]).

The following constrained approximation problem was recently considered in [@LKW17].

\[P:approx\] Given the Bézier coefficients $r_{{\mbox{\scriptsize$\bm{k}$}}}$ and positive weights $\omega_{{\mbox{\scriptsize$\bm{k}$}}}$, $\bm k\in\Theta_n$, of the rational function ${\mbox{$\mathsf R$}}_n$ of degree $n$, $$\label{E:Rn}    {\mbox{$\mathsf R$}}_n(\bm x)    
    :=\frac{\displaystyle \sum_{{\mbox{\scriptsize$\bm{k}$}}\in\Theta_n}\omega_{{\mbox{\scriptsize$\bm{k}$}}}r_{{\mbox{\scriptsize$\bm{k}$}}}B^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x)}
    {\displaystyle \sum_{{\mbox{\scriptsize$\bm{k}$}}\in\Theta_n}\omega_{{\mbox{\scriptsize$\bm{k}$}}}B^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x)}
    = \sum_{{\mbox{\scriptsize$\bm{k}$}}\in\Theta_n}r_{{\mbox{\scriptsize$\bm{k}$}}}Q^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x),\qquad \bm x\in T,$$ where $$\label{E:Q}    Q^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x):=\frac{\displaystyle \omega_{{\mbox{\scriptsize$\bm{k}$}}}B^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x)}
    {\displaystyle \sum_{{\mbox{\scriptsize$\bm{i}$}}\in\Theta_n}\omega_{{\mbox{\scriptsize$\bm{i}$}}}B^n_{{\mbox{\scriptsize$\bm{i}$}}}(\bm x)},$$ find a polynomial of degree $m$, of the form $$\label{E:Pm}    {\mbox{$\mathsf P$}}_m(\bm x)
    := \sum_{{\mbox{\scriptsize$\bm{k}$}}\in\Theta_m}p_{{\mbox{\scriptsize$\bm{k}$}}}B^m_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x),\qquad \bm x\in T,$$ with the coefficients $p_{{\mbox{\scriptsize$\bm{k}$}}}$ satisfying the conditions $$\label{E:gcond}    p_{{\mbox{\scriptsize$\bm{k}$}}}=g_{{\mbox{\scriptsize$\bm{k}$}}}\quad \mbox{for}\quad \bm k\in\Gamma^{{\mbox{\scriptsize$\bm{c}$}}}_m,$$ $g_{{\mbox{\scriptsize$\bm{k}$}}}$ being prescribed numbers, and $\bm c:=(c_1,c_2,c_3)\in{\mathbb N}^3$ being a given parameter vector with $|\bm c|<m$, such that the distance between ${\mbox{$\mathsf R$}}_n$ and ${\mbox{$\mathsf P$}}_m$, $$\label{E:dist}    d({\mbox{$\mathsf R$}}_n,{\mbox{$\mathsf P$}}_m):={\mbox{$\displaystyle\int\!\!\!\!\int\limits_{\!\!\!\!T}$}}w_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x)[{\mbox{$\mathsf R$}}_n(\bm x)-{\mbox{$\mathsf P$}}_m(\bm x)]^2{\,\mbox{{\rm d}$\bm x$}},$$ reaches the minimum.

It has been shown that the solution of Problem \[P:approx\] is the polynomial with the coefficients given by (cf. [@LKW17 Thm 2.2]) $$\label{E:psc}    p_{{\mbox{\scriptsize$\bm{k}$}}}=\sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_m}\binom{m}{\bm l}\,E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},m)
\big(u_{{\mbox{\scriptsize$\bm{l}$}}}-v_{{\mbox{\scriptsize$\bm{l}$}}}\big), \qquad \bm k\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_m,$$ where $$\begin{aligned}
      u_{{\mbox{\scriptsize$\bm{l}$}}}:=& \sum_{{\mbox{\scriptsize$\bm{h}$}}\in\Theta_n} \binom{n}{\bm h}
    {\mbox{$\displaystyle\binom{n+m}{\bm h+\bm l}^{\!\!-1}\!\!$}}\,r_{{\mbox{\scriptsize$\bm{h}$}}}\,I_{{\mbox{\scriptsize$\bm{h}$}},{\mbox{\scriptsize$\bm{l}$}}},\\ 
        \label{E:vl}    v_{{\mbox{\scriptsize$\bm{l}$}}}:=& \frac{1}{(|{\mbox{$\bm{\alpha}$}}|+3)_{2m}}\sum_{{\mbox{\scriptsize$\bm{h}$}}\in\Gamma^{{\mbox{\scriptsize$\bm{c}$}}}_m}\binom{m}{\bm h}                   
           \left(\prod_{i=1}^{3}(\alpha_i+1)_{h_i+l_i}\right)g_{{\mbox{\scriptsize$\bm{h}$}}}             \end{aligned}$$ with $h_3:=m-|\bm h|$, $l_3:=m-|\bm l|$, and $$\label{E:I}    I_{{\mbox{\scriptsize$\bm{h}$}},{\mbox{\scriptsize$\bm{l}$}}}:=  {\mbox{$\displaystyle\int\!\!\!\!\int\limits_{\!\!\!\!T}$}}w_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x)Q^{n}_{{\mbox{\scriptsize$\bm{h}$}}}(\bm x)
    B^m_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x){\,\mbox{{\rm d}$\bm x$}}.$$ The symbol $E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},m)$ has the meaning given in (see also ).

Implementation of formula demands $1^o$ evaluation of all the coefficients $E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},m)$ with $\bm k, \bm l \in\Omega^{{\mbox{${\mbox{\scriptsize$\bm{c}$}}$}}}_m$ and $2^o$ computing all the integrals $I_{{\mbox{\scriptsize$\bm{h}$}},{\mbox{\scriptsize$\bm{l}$}}}$ with $\bm h\in\Theta_n$, $\bm l\in\Omega^{{\mbox{${\mbox{\scriptsize$\bm{c}$}}$}}}_m$.

The first task is accomplished in two steps. (*i*) We compute all the the coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}(\bm\mu,{\mbox{$\bm c$}},M)$, $\bm k, \bm l \in\Omega^{{\mbox{${\mbox{\scriptsize$\bm{c}$}}$}}}_m$, with $M:=m-{\mbox{$\left|{\mbox{$\bm c$}}\right|$}}$ and $\bm\mu:={\mbox{$\bm{\alpha}$}}+2{\mbox{$\bm c$}}$ by Algorithm \[A:E-comp\] and then (*ii*) use the result in formula (with $n$ replaced by $m$) to compute $E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},m)$, $\bm k, \bm l \in\Omega^{{\mbox{${\mbox{\scriptsize$\bm{c}$}}$}}}_m$. Remark that the striking simplicity of Algorithm \[A:E-comp\] was obtained thanks to using properties of dual bivariate Bernstein polynomials, investigated in Appendix A.

As for the second task, observe that in general, the integrals cannot be evaluated exactly, and their number, equal to $(n+m-|\bm c|)(n+m-|\bm c|+1)/2$, may lead to the impression that computation of the coefficients is time-expensive in practice. The problem has been defeated in [@LKW17], where we have proposed a very effective algorithm which allows to compute the complete set of integrals $I_{{\mbox{\scriptsize$\bm{h}$}},{\mbox{\scriptsize$\bm{l}$}}}$, $\bm h\in\Theta_n$, $\bm l\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_m$, in a time about twice as long as the time required to compute a single integral of this type (the algorithm is an extension of the adaptive, high precision method of [@Kel07], and it uses some special properties of the integrals ).

In this way, we have obtained a very efficient method for computing the solution to the Problem \[P:approx\].

 Case of equal weights {#SS:eqw}
----------------------

Notice that in the particular case where all the weights are equal, $\omega_{{\mbox{\scriptsize$\bm{i}$}}}=\omega$, $\bm i\in\Theta_n$, we have $Q^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x)=B^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x)$, $\bm k\in\Theta_n$. Consequently, the rational function reduces to a polynomial of degree $n$, so that Problem \[P:approx\] (with additional assumption $m<n$) is actually the constrained polynomial degree reduction problem that has been discussed [@WL10]. In the cited paper, we have developed a method to evaluate the control points of the reduced surface in terms of quantities $\left\langle D^{(m,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}},\,B^n_{{\mbox{\scriptsize$\bm{l}$}}} \right\rangle_{{\mbox{\scriptsize$\bm{a}$}}}$, $\bm k\in \Omega^{{\mbox{\scriptsize$\bm{c}$}}}_{m}$, $\bm l\in\Theta_n$. Computing these coefficients was the most cumbersome part of the method. The main part of the above collection was computed efficiently thanks to using some recurrence relations. However, a quite large portion of the coefficients demanded some complex computations. By applying the approach of this paper, we may solve the degree reduction problem much faster.

N.B. In the case of equal weights, evaluation of the integrals is a simple task as we have (see, e.g., [@WL10]) $$I_{{\mbox{\scriptsize$\bm{h}$}},{\mbox{\scriptsize$\bm{l}$}}}=\binom{n}{\bm h}\binom{m}{\bm l}
    \frac{
    (\alpha_1+1)_{h_1+l_1}(\alpha_2+1)_{h_2+l_2}(\alpha_3+1)_{n+m-|{\mbox{\scriptsize$\bm{h}$}}|-|{\mbox{\scriptsize$\bm{l}$}}|}}
    {(|{\mbox{$\bm{\alpha}$}}|+3)_{n+m}}.$$

Appendix A: Dual bivariate Bernstein bases {#S:AppA .unnumbered}
==========================================

0 0 0

Unconstrained dual bivariate Bernstein polynomials {#SS:dualB2}
--------------------------------------------------

In this section, we prove some important properties of the coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$ introduced in . Obviously, we have (cf. ) $$\label{E:e-is}    e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=\langle D^n_{{\mbox{\scriptsize$\bm{k}$}}}, D^n_{{\mbox{\scriptsize$\bm{l}$}}}\rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}, \qquad
    \bm k,\bm l\in\Theta_n.$$ However, we need some alternative formulas that would be more useful in obtaining recurrence relations satisfied by these quantities.

In the sequel, we adopt the following convention: given $\bm t:=(t_1,t_2)\in\Theta_n$, we use the notation $$\tilde{\bm t}:=(t_2,t_1),\quad \hat{\bm t}:=(t_1,t_3),\quad \check{\bm t}:=(t_3,t_1),\quad \bm t^\ast:=(t_2,t_3)\quad \mbox{and}\quad {\bm t}^{\circ}:=(t_3,t_2),$$ where $t_3:=n-|\bm t|$.

\[L:D2inB2\] Dual bivariate Bernstein polynomials have the following representation: $$D^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})=\,
    \sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Theta_n}\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)\,
    B^n_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x),\qquad \bm k\in\Theta_n,$$ where $$\label{E:D2inB2-coef}    e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n):=\sum_{0\le i\le q\le n}\,C^2_{q,i}({\mbox{$\bm{\alpha}$}},n)
    H_{q,i}\left(\bm k^\ast;{\mbox{$\bm{\alpha}$}},n\right)H_{q,i}\left(\bm l^\ast;{\mbox{$\bm{\alpha}$}},n\right)$$ with $$C_{q,i}({\mbox{$\bm{\alpha}$}},n):=\binom{q}{i}\frac{1}{q!(-n)_q\, \lambda_{q,i}}.$$ Here $H_{q,i}(\bm t;{\mbox{$\bm{\alpha}$}},n)$ are the bivariate Hahn polynomials defined by , and $\lambda_{q,i}$ is the constant given by .

In [@LW06], it has been proved that the dual bivariate Bernstein polynomials have the following representation: $$\label{E:D2inP2}    D^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})=\,
    \sum_{0\le i\le q\le n}\,b^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{q,i}(n,\bm k)\,
    P_{q,i}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x),\qquad \bm k\in\Theta_n,$$ where $$b^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{ q,i}(n,\bm k):=
    (-1)^{i}C_{q,i}({\mbox{$\bm{\alpha}$}},n)
    H_{q,i}\left(\bm k^\ast;{\mbox{$\bm{\alpha}$}},n\right).$$ Here $P_{q,i}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}$, $q\in{\mathbb N};\; 0\le i\le q$, are the *two-variable Jacobi polynomials* ([@Koo75]; see also [@DX01 p. 86] or ), which are orthonormal on $T$, i.e., $\langle P_{m,l}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}},\, P_{n,k}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}} \rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}$ equals 1 if $(m,l)=(n,k)$, and 0 otherwise. Using expansion for both dual polynomials in , we obtain $$e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)
    =\sum_{0\le i\le q\le n}b^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{q,i}(n,\bm k)b^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{q,i}(n,\bm l),$$ which is equivalent to .

\[R:ekl-sym\] Some symmetry properties of the coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$ should be noticed.

1.  The following equation results easily from : $$\label{E:sym1}    e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{{\mbox{\scriptsize$\bm{l}$}}}_{{\mbox{\scriptsize$\bm{k}$}}}({\mbox{$\bm{\alpha}$}},n).$$

2.  Let $\alpha_2=\alpha_3$. The following equality holds: $$\label{E:sym2}        e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{\hat{{\mbox{\scriptsize$\bm{k}$}}}}_{\hat{{\mbox{\scriptsize$\bm{l}$}}}}({\mbox{$\bm{\alpha}$}},n).$$ This equation can be easily verified using , definition of the bivariate Hahn polynomials, and the identity $
                h_i(n-|\bm t|;\alpha_2,\alpha_2,n-t_1)=(-1)^ih_i(t_2;\alpha_2,\alpha_2,n-t_1)
                $ (cf. ).

3.  The following equation holds: $$\label{E:sym3}    e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}(\tilde{{\mbox{$\bm{\alpha}$}}},n),$$ where $\tilde{{\mbox{$\bm{\alpha}$}}}:=(\alpha_2,\alpha_1,\alpha_3)$. By , equation is equivalent to the equation $$\left\langle   D^n_{{\mbox{\scriptsize$\bm{k}$}}}(\cdot;{\mbox{$\bm{\alpha}$}}),D^n_{{\mbox{\scriptsize$\bm{l}$}}}(\cdot;{\mbox{$\bm{\alpha}$}}) \right \rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}
        = \left\langle D^n_{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}(\cdot;\tilde{{\mbox{$\bm{\alpha}$}}}),D^n_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}(\cdot;\tilde{{\mbox{$\bm{\alpha}$}}})\right\rangle_{{\mbox{\scriptsize$\tilde{\bm{\alpha}}$}}}$$ which can be easily verified using the definition of the inner product $\langle\cdot,\cdot\rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}$.

4.  Let $\alpha_1=\alpha_2$. By , we obtain the equation $$\label{E:sym4}        e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}({\mbox{$\bm{\alpha}$}},n).$$

5.  Let $\alpha_1=\alpha_3$. Then $$\label{E:sym5}        e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{{{\mbox{\scriptsize$\bm{k}$}}}^{\circ}}_{{{\mbox{\scriptsize$\bm{l}$}}}^{\circ}}({\mbox{$\bm{\alpha}$}},n).$$ Transforming $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$, using consecutively , , and again , the result follows.

6.  Let $\alpha_1=\alpha_2=\alpha_3$. By , , and , $$\label{E:sym6}        e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}({\mbox{$\bm{\alpha}$}},n)=e^{\hat{{\mbox{\scriptsize$\bm{k}$}}}}_{\hat{{\mbox{\scriptsize$\bm{l}$}}}}({\mbox{$\bm{\alpha}$}},n)
            =e^{{\mbox{\scriptsize$\bm{k}$}}^\ast}_{{\mbox{\scriptsize$\bm{l}$}}^\ast}({\mbox{$\bm{\alpha}$}},n)=e^{\check{{\mbox{\scriptsize$\bm{k}$}}}}_{\check{{\mbox{\scriptsize$\bm{l}$}}}}({\mbox{$\bm{\alpha}$}},n)
            =e^{{{\mbox{\scriptsize$\bm{k}$}}}^{\circ}}_{{{\mbox{\scriptsize$\bm{l}$}}}^{\circ}}({\mbox{$\bm{\alpha}$}},n).$$

\[L:e-recur\] The coefficients $e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}\equiv e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$ defined in , and given by , satisfy the following bivariate recurrence relations: $$\begin{aligned}
    \label{E:e-receq1}            &\sigma_0(\bm k)\, e^{{\mbox{\scriptsize$\bm{k}$}}+{\mbox{\scriptsize$\bm{v}$}}_2}_{{\mbox{\scriptsize$\bm{l}$}}}-\sigma_1(\bm k)\, e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
            +\sigma_2(\bm k) \,e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{v}$}}_2}_{{\mbox{\scriptsize$\bm{l}$}}}         
            =\sigma_0(\bm l)\, e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}+{\mbox{\scriptsize$\bm{v}$}}_2}-\sigma_1(\bm l) \,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
            +\sigma_2(\bm l)\, e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{v}$}}_2},\\[1ex] 
    \label{E:e-receq2}    & \tau_0(\bm k)\,e^{{\mbox{\scriptsize$\bm{k}$}}+{\mbox{\scriptsize$\bm{v}$}}_1}_{{\mbox{\scriptsize$\bm{l}$}}}
                  -\tau_1(\bm k)\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
                  +\tau_2(\bm k)\,e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{v}$}}_1}_{{\mbox{\scriptsize$\bm{l}$}}}  
                =
                  \tau_0(\bm l)\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}+{\mbox{\scriptsize$\bm{v}$}}_1}
                  - \tau_1(\bm l)\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
                  + \tau_2(\bm l)\,e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{v}$}}_1},      
                  \end{aligned}$$ where $\bm v_1:=(1,0)$, $\bm v_2:=(0,1)$, and for $\bm t:=(t_1,t_2)$, we define $$\label{E:sigtau}    \sigma_i(\bm t):=\varphi_i({\mbox{$\bm{\alpha}$}},\bm t),\quad \tau_i(\bm t):=\varphi_i(\tilde{\mbox{$\bm{\alpha}$}},\tilde{\bm t}),\qquad i=0,1,2,$$ with $\tilde{\mbox{$\bm{\alpha}$}}:=(\alpha_2,\alpha_1,\alpha_3)$, and $$\label{E:varphi}    \begin{array}{c}
    \varphi_0({\mbox{$\bm{\alpha}$}},\bm t):=(|\bm t|-n)(t_2+\alpha_2+1),\quad    
    \varphi_2({\mbox{$\bm{\alpha}$}},\bm t):=t_2(|\bm t|-\alpha_3-n-1),\\[1ex]
    \varphi_1({\mbox{$\bm{\alpha}$}},\bm t):=\varphi_0({\mbox{$\bm{\alpha}$}},\bm t)+\varphi_2({\mbox{$\bm{\alpha}$}},\bm t).
    \end{array}$$

First, we prove the recurrence . By , we have $$H_{q,i}\left(\bm t^\ast;{\mbox{$\bm{\alpha}$}},n\right)=h_i(t_2;\alpha_2,\alpha_3,n-t_1)\,
    h_{q-i}(n-t_1-i;\alpha_2+\alpha_3+2i+1,\alpha_1,n-i).$$ Let ${\mbox{$\mathcal D$}}^{n-t_1}_{t_2}$ be the difference operator defined according to . Then, by , $${\mbox{$\mathcal D$}}^{n-t_1}_{t_2}\,h_i(t_2;\alpha_2,\alpha_3,n-t_1)
            =i(i+\alpha_2+\alpha_3+1)\,h_i(t_2;\alpha_2,\alpha_3,n-t_1).$$ Hence, we obtain $$\begin{aligned}
                {\mbox{$\mathcal D$}}^{n-t_1}_{t_2}\,H_{q,i}\left(\bm t^\ast;{\mbox{$\bm{\alpha}$}},n\right)&=
            {\mbox{$\mathcal D$}}^{n-t_1}_{t_2}\,h_i(t_2;\alpha_2,\alpha_3,n-t_1)\cdot 
            h_{q-i}(n-t_1-i;\alpha_2+\alpha_3+2i+1,\alpha_1,n-i)\\
            &=i(i+\alpha_2+\alpha_3+1) \,H_{q,i}\left(\bm t^\ast;{\mbox{$\bm{\alpha}$}},n\right).
      \end{aligned}$$ Consequently, having in mind the form , we obtain the equation $${\mbox{$\mathcal D$}}^{n-k_1}_{k_2}e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n) = {\mbox{$\mathcal D$}}^{n-l_1}_{l_2}e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)$$ which can be simplified to the form .

To prove , let us define the difference operator ${\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{{\mbox{\scriptsize$\bm{t}$}}}$ by $${\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{{\mbox{\scriptsize$\bm{t}$}}}\,F(\bm t)=\varphi_0({\mbox{$\bm{\alpha}$}},\bm t)F(\bm t+\bm v_2)
        -\varphi_1({\mbox{$\bm{\alpha}$}},\bm t)F(\bm t)
        +\varphi_2({\mbox{$\bm{\alpha}$}},\bm t)F(\bm t-\bm v_2),$$ where we use the notation of . The recurrence can be written as $${\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{{\mbox{\scriptsize$\bm{k}$}}}e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n)=
        {\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n).$$ Substituting $\tilde{\bm k}$, $\tilde{\bm l}$ and ${\mbox{$\tilde{{\mbox{$\bm{\alpha}$}}}$}}$ in place of $\bm k$, $\bm l$ and ${\mbox{$\bm{\alpha}$}}$, respectively, gives $$\label{E:Rt}                {\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{\scriptsize$\tilde{\bm{\alpha}}$}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}e^{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}({\mbox{$\tilde{{\mbox{$\bm{\alpha}$}}}$}},n)=
            {\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{\scriptsize$\tilde{\bm{\alpha}}$}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}e^{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}({\mbox{$\tilde{{\mbox{$\bm{\alpha}$}}}$}},n).$$ Now, notice that by we have $$\begin{aligned}
                  & 
                  {\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{\scriptsize$\tilde{\bm{\alpha}}$}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}e^{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}(\tilde{{\mbox{$\bm{\alpha}$}}},n)
                  =\varphi_0(\tilde{{\mbox{$\bm{\alpha}$}}},\tilde{\bm k})e^{{\mbox{\scriptsize$\bm{k}$}}+{\mbox{\scriptsize$\bm{v}$}}_1}_{{\mbox{\scriptsize$\bm{l}$}}}
                  -\varphi_1(\tilde{\mbox{$\bm{\alpha}$}},\tilde{\bm k})e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
                  +\varphi_2(\tilde{\mbox{$\bm{\alpha}$}},\tilde{\bm k})e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{v}$}}_1}_{{\mbox{\scriptsize$\bm{l}$}}}, \\
                  &{\mbox{$\mathcal R$}}^{{\mbox{\scriptsize${\mbox{\scriptsize$\tilde{\bm{\alpha}}$}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}e^{{\mbox{\scriptsize$\tilde{\bm{k}}$}}}_{{\mbox{\scriptsize$\tilde{\bm{l}}$}}}(\tilde{{\mbox{$\bm{\alpha}$}}},n)
                  =\varphi_0(\tilde{\mbox{$\bm{\alpha}$}},\tilde{\bm l})e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}+{\mbox{\scriptsize$\bm{v}$}}_1}
                  -\varphi_1(\tilde{\mbox{$\bm{\alpha}$}},\tilde{\bm l})e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}
                  +\varphi_2(\tilde{\mbox{$\bm{\alpha}$}},\tilde{\bm l})e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{v}$}}_1}.
            \end{aligned}$$ Using this in , equation follows.

\[L:e-k20\] For $\bm k=(k_1,0),\,\bm l=(l_1,l_2)\in\Theta_n$, the following formula holds: $$\label{E:e-k20}    e^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n) =
     G_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}\sum_{i=0}^{n-\max(k_1,l_1)}C_i\,h_i(l_2;\alpha_2,\alpha_3,n-l_1),$$ where $ G_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}} :={\Gamma(\eta+1)\Gamma(\alpha_1+1)}/{\Gamma(|{\mbox{$\bm{\alpha}$}}|+3)}$, $\eta:=\alpha_2+\alpha_3+1$, symbol $h_i(t;a,b,M)$ is defined in , and $$\label{E:Ci}                \left.\begin{array}{ll}
                C_0:=c_{k_1,l_1}(n,\eta,\alpha_1),\\[1ex]
                C_i:=\dfrac{(2i+\eta)(k_1-n)_i(\eta+1)_{i-1}}{i!(-n)_i^2(\alpha_3+1)_i}
                \,c_{k_1,l_1}(n-i,2i+\eta,\alpha_1),&\quad i\ge1,   
                \end{array}\quad\right\}$$ the symbol $c_{h,j}(m,\alpha,\beta)$ denoting the $j$*th* Bézier coefficient of the univariate dual Bernstein polynomial $D^{m}_{h}(x;\alpha,\beta)$ (cf. [@LW11b]).

We give a sketch of the proof. By [@LW06 Thm 3.3], we have for $\bm k:=(k_1,k_2)\in\Theta_n$, $$D^n_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})
    =G_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}\sum_{i=0}^{n-k_1}f_i(n,\bm k)
    (1-x_1)^iR^{(\alpha_3,\alpha_2)}_i\left(x_2/(1-x_1)\right)D^{n-i}_{k_1}(x_1;2i+\eta,\alpha_1),$$ where $R^{(\alpha_3,\alpha_2)}_i$ are the univariate Jacobi polynomials (cf. ), $D^N_j(t;\mu,\nu)$ are univariate dual Bernstein polynomials (see, e.g., [@LW11b]), and $$\begin{array}{ll}
        f_0(n,\bm k):=1,\\[0.5ex]
        f_i(n,\bm k):=(-1)^i\dfrac{(2i+\eta)(\eta+1)_{i-1}}
                           {(-n)_i(\alpha_2+1)_i(\alpha_3+1)_i}\,
                           \,h_i(k_2;\alpha_2,\alpha_3,n-k_1),&\quad i\ge1.
                           \end{array}$$ We use the above formula in , then reduce the integration over the triangle $T$ to evaluating two one-dimensional integrals, and apply the orthogonality property of the polynomials $R^{(\alpha_3,\alpha_2)}_i$ (cf. [@KS98 §1.8]). Letting $\bm k:=(k_1,0)$, the formula follows.

\[C:e00\] The following formula holds: $$e^{(0,0)}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},n) =\frac{(-1)^{l_1}(|{\mbox{$\bm{\alpha}$}}|+3)_{n}}{n!(\alpha_1+2)_{l_1}}
    \sum_{i=0}^{n-l_1}C^\ast_i\,h_i(l_2;\alpha_2,\alpha_3,n-l_1),
    \qquad \bm l:=(l_1,l_2)\in\Theta_n,$$ where $$\label{E:Cast}                \left.\begin{array}{ll}
                C^\ast_0:=\dfrac{(\alpha_1+2)_{n}}
                {(\alpha_2+\alpha_3+2)_{n-l_1}},\\[2.5ex]
                \displaystyle C^\ast_i:=(-1)^{i}\dfrac{(2i+\alpha_2+\alpha_3+1)(\alpha_1+2)_{n-i}(|{\mbox{$\bm{\alpha}$}}|+n+3)_{i}}
                {i!(\alpha_3+1)_i(\alpha_2+\alpha_3+i+1)_{n-l_1+1}},&\quad i\ge1.   
                \end{array}\quad\right\}$$

The result follows by putting $k_1=0$ in , , and using the explicit form for $c_{0,l_1}(n-i,2i+\alpha_2+\alpha_3+1,\alpha_1)$, given in [@LW11b Eq. (2.11)].

Constrained dual bivariate Bernstein polynomials {#SS:cdualA2}
------------------------------------------------

The constrained dual bivariate Bernstein basis polynomials can be expressed in terms of the unconstrained dual bivariate Bernstein polynomials with shifted degree and parameters. Namely, we have the following result.

\[L:constrDual2\] For $\bm k\in \Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n$, the following formula holds: $$\label{E:Dc2inD2}    D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})=
    U\,V_{{\mbox{\scriptsize$\bm{k}$}}}(n)\,
    x_1^{c_1}x_2^{c_2}(1-|\bm x|)^{c_3}\,
    D^{n-|{\mbox{\scriptsize$\bm{c}$}}|}_{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{c}$}}'}(\bm x;{\mbox{$\bm{\alpha}$}}+2\bm c),$$ where $\bm c':=(c_1,c_2)$, and $$\label{E:UV}     U:=
    (|{\mbox{$\bm{\alpha}$}}|+3)_{2|{\mbox{\scriptsize$\bm{c}$}}|}\prod_{i=1}^{3}
    (\alpha_i+1)_{2c_i}^{-1},\qquad
    V_{{\mbox{\scriptsize$\bm{k}$}}}(n):= \binom{n-|\bm c|}{\bm k-\bm c'}{\mbox{$\displaystyle\binom{n}{\bm k}^{\!\!-1}\!\!$}}.$$

\[L:Dc2inB2\] The constrained dual bivariate Bernstein polynomials have the Bézier representation $$D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})=\sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_n}
    \,E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},n)\,B^n_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x),$$ where $$E^{{\mbox{\scriptsize$\bm{k}$}}}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}},{\mbox{$\bm c$}},n):=U\,V_{{\mbox{\scriptsize$\bm{k}$}}}(n)\,V_{{\mbox{\scriptsize$\bm{l}$}}}(n)\,
    \,e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{c}$}}'}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{c}$}}'}({\mbox{$\bm{\alpha}$}}+2\bm c,n-|\bm c|),\qquad \bm c':=(c_1,c_2),$$ notation used being that of and .

By Lemma \[L:D2inB2\], we have $$D^{n-|{\mbox{\scriptsize$\bm{c}$}}|}_{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{c}$}}'}(\bm x;{\mbox{$\bm{\alpha}$}}+2\bm c)
        =\sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Theta_{n-|{\mbox{\tiny$\bm{c}$}}|}}    
        e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{c}$}}'}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}}+2\bm c,n-|\bm c|)\,
        B^{n-|{\mbox{\scriptsize$\bm{c}$}}|}_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x).$$ Putting this result in and using the equation $$x_1^{c_1}x_2^{c_2}(1-|\bm x|)^{c_3}\cdot B^{n-|{\mbox{\scriptsize$\bm{c}$}}|}_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x)=
        \binom{n-|\bm c|}{\bm l}\binom{n}{\bm l+\bm c'}^{\!-1}\,B^n_{{\mbox{\scriptsize$\bm{l}$}}+{\mbox{\scriptsize$\bm{c}$}}'}(\bm x),$$ we obtain $$\begin{aligned}
    D^{(n,{\mbox{\scriptsize$\bm{c}$}})}_{{\mbox{\scriptsize$\bm{k}$}}}(\bm x;{\mbox{$\bm{\alpha}$}})=&
    U\,V_{{\mbox{\scriptsize$\bm{k}$}}}(n)\,\sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Theta_{n}}    
    \binom{n-|\bm c|}{\bm l}{\mbox{$\displaystyle\binom{n}{\bm l+\bm c'}^{\!\!-1}\!\!$}}
    \,e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{c}$}}'}_{{\mbox{\scriptsize$\bm{l}$}}}({\mbox{$\bm{\alpha}$}}+2\bm c,n-|\bm c|)\,
    B^{n}_{{\mbox{\scriptsize$\bm{l}$}}+{\mbox{\scriptsize$\bm{c}$}}'}(\bm x)\\
    =&
    U\,V_{{\mbox{\scriptsize$\bm{k}$}}}(n)\,\sum_{{\mbox{\scriptsize$\bm{l}$}}\in\Omega^{{\mbox{\scriptsize$\bm{c}$}}}_{n}}
    \binom{n-|\bm c|}{\bm l-\bm c'}{\mbox{$\displaystyle\binom{n}{\bm l}^{\!\!-1}\!\!$}}    
    \,e^{{\mbox{\scriptsize$\bm{k}$}}-{\mbox{\scriptsize$\bm{c}$}}'}_{{\mbox{\scriptsize$\bm{l}$}}-{\mbox{\scriptsize$\bm{c}$}}'}({\mbox{$\bm{\alpha}$}}+2\bm c,n-|\bm c|)\,
    B^{n}_{{\mbox{\scriptsize$\bm{l}$}}}(\bm x).\end{aligned}$$ Hence, the lemma is proved.

Appendix B: Bivariate orthogonal polynomials {#S:AppB .unnumbered}
============================================

0 0 0

The notation $${\,{}_{r}F_{s}\!\left(                \begin{array}{cc}{a_1,\ldots, a_r}\\[0.35ex]{b_1,\ldots,b_s} \end{array}\Big|\,{z} \right)}
        :=\sum_{k=0}^{\infty}\frac{(a_1)_k\cdots (a_r)_k}
                           {k!(b_1)_k\cdots (b_s)_k}\,z^k$$ is used for the *generalised hypergeometric series* (see, e.g., [@AAR99 §2.1]); here $r,\,s\in{{\mathbb Z}_+}$, $z$, $a_1,\ldots, a_r$, $b_1,\ldots,\,b_s$$\in {\mathbb C}$, and $(c)_k$ is the shifted factorial.

Recall that *two-variable Jacobi polynomials* $P_{n,k}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x)$, $n=0,1,\ldots,\; k=0,1,\ldots ,n$, are defined by ([@Koo75]; see also [@DX01 p. 86]) $$\begin{aligned}
\label{E:Jacxy}    P_{n,k}^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x):=& \lambda^{-1}_{n,k}
         \,R_{n-k}^{(2k+\alpha_2+\alpha_3+1,\alpha_1)}(x_1)
         \,(1-x_1)^k\,
         R_k^{(\alpha_3 ,\alpha_2)}\left({x_2}/{(1-x_1)}\right),\end{aligned}$$ where $\bm x:=(x_1,x_2)$, ${\mbox{$\bm{\alpha}$}}:=(\alpha_1,\,\alpha_2,\,\alpha_3)$, $\alpha_i>-1$, $$\label{E:uniJac}R^{(\mu,\nu)}_m(t):=\frac{(\mu+1)_m}{m!}{\,{}_{2}F_{1}\!\left(                \begin{array}{cc}{-m,m+\mu+\nu+1}\\[0.35ex]{\mu+1} \end{array}\Big|\,{1-t} \right)}$$ is the $m$th shifted Jacobi polynomial in one variable [@KS98 §1.8], and $$\label{E:lam} \lambda^{2}_{n,k}\equiv
 \left[ \lambda^{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}_{n,k}\right]^{2}:=
 \frac{(\alpha_1+1)_{n-k}
 (\alpha_2+1)_k(\alpha_3+1 )_k(k+\eta)_{n+1}}
 {k!(n-k)!(2k+\eta)(2n/\sigma+1)(\sigma)_{n+k}}$$ with $\eta:=\alpha_2+\alpha_3+1$, $\sigma:=|{\mbox{$\bm{\alpha}$}}|+2$. Polynomials form the orthonormal set with respect to the inner product $$\langle f, g \rangle_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}} :={\mbox{$\displaystyle\int\!\!\!\!\int\limits_{\!\!\!\!T}$}}w_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x)f(\bm x)\,g(\bm x) {\,\mbox{{\rm d}$\bm x$}},$$ where $T:=\{(x_1,\,x_2)\;:\;x_1,\,x_2\ge0,\; 1-x_1-x_2\ge0\}$, and $w_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}(\bm x):=A_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}\, x_1^{\alpha_1} x_2^{\alpha_2}(1-x_1-x_2)^{\alpha_3}$ with $$\label{E:A}A_{{\mbox{\scriptsize${\mbox{$\bm{\alpha}$}}$}}}:=\Gamma({\mbox{$\left|{\mbox{$\bm{\alpha}$}}\right|$}}+3)\prod_{i=1}^{3}[\Gamma(\alpha_i+1)]^{-1}.$$ *Bivariate Hahn polynomials* are defined by (see, e.g., [@Tra91]) $$\label{E:Hahn2}    H_{q,i}(\bm t;{\mbox{$\bm{\alpha}$}},n):=h_i(t_1;\alpha_2,\alpha_3,t_1+t_2)
    \,h_{q-i}(t_1+t_2-i;\alpha_2+\alpha_3+2i+1,\alpha_1,n-i)\,,$$ where $0\le i\le q\le n$, $n\in{\mathbb N}$, ${\mbox{$\bm{\alpha}$}}=:(\alpha_1,\alpha_2\,\alpha_3)$, $\alpha_i>-1$, $i=1,2,3$, $\bm t:=(t_1,t_2)\in{\mathbb R}^2$, and $$\label{E:Hahn1}    h_l(t;a,b,M):=(a+1)_l(-M)_l{\,{}_{3}F_{2}\!\left(                \begin{array}{cc}{-l,l+a+b+1,-t}\\[0.35ex]{a+1,-M} \end{array}\Big|\,{1} \right)}$$ are the *univariate Hahn polynomials* (see, e.g., [@KS98 §1.5]). The latter polynomials satisfy the recurrence relation $$h_{l+1}(t)=A_l(t,M)\,h_l(t)+B_l(M)\,h_{l-1}(t),
        \qquad l\ge0;\; h_{0}(t)\equiv1;    \;h_{-1}(t)\equiv0,$$ where $h_l(t)\equiv h_l(t;a,b,M)$, $$\label{E:Hahn1-rec-coeffs}    A_l(t,M):=C_l\,(2l+s-1)_2\,t-D_l-E_l,\qquad     
 B_l(M):=-D_l\,E_{l-1},$$ with $s:=a+b+1$, $C_l:=(2l+s+1)/[(l+s)(2l+s-1)]$, $D_l:=C_l\,l(l+M+s)(l+b)$, and $E_l:=(l+a+1)(M-l)$. Moreover, the polynomials satisfy the difference equation $$\label{E:Hahn1deq}
    {\mbox{$\mathcal D$}}^M_t\,h_j(t;a,b,M)=j(j+a+b+1)\,h_j(t;a,b,M),$$ where the difference operator ${\mbox{$\mathcal D$}}^M_t$ is given by $$\label{E:Dop}    {\mbox{$\mathcal D$}}^M_t F(t):=U(t;a,b,M)\,F(t+1)-V(t;a,b,M)\,F(t)+W(t;a,b,M)\,F(t-1)$$ with $U(t;a,b,M):=(t-M)(t+a+1)$, $W(t;a,b,M):=t(t-b-M-1)$ and $V(t;a,b,M):=U(t;a,b,M)+W(t;a,b,M)$.

\[R:Clenshaw\] A linear combination $
        s_N(t):=\sum_{i=0}^{N}\gamma_ih_i(t;a,b,M)
    $ can be summed up using the following *Clenshaw’s algorithm* (see, e.g., [@DB08 Thm 3.2.11]). Compute the sequence $V_0,V_1,\ldots$, $V_{n+2}$ from $$V_i:=\gamma_i+A_i(t;M)V_{i+1}+B_{i+1}(M)V_{i+2},\qquad i=N,N-1,\ldots,0,$$ with $V_{N+1}=V_{N+2}=0$, where the coefficients $A_i(t;M)$ and $B_i(M)$ are defined by . Then, $s_N(t)=V_0$.

[99]{}

Y.J. Ahn, Using Jacobi polynomials for degree reduction of Bézier curves with $C^k$-constraints, Comput. Aided Geom. Design 20 (2003) 423–434.

Y.J. Ahn, B.-G. Lee, Y. Park, J. Yoo, Constrained polynomial degree reduction in the $L_2$-norm equals best weighted Euclidean approximation of Bézier coefficients, Comput. Aided Geom. Design 21 (2004) 181–191.

G.E. Andrews, R. Askey, R. Roy, Special Functions, Cambridge Univ. Press, Cambridge, 1999.

G. Chen, G. Wang, Optimal multi-degree reduction of Bézier curves with constraints of endpoints continuity, Comput. Aided Geom. Design 19 (2002) 365–377.

G. Dahlquist, A. Björck, Numerical Methods in Scientific Computing, vol. I, SIAM, Philadelphia, 2008.

C.F. Dunkl, Y. Xu, Orthogonal Polynomials of Several Variables, Cambridge Univ. Press, Cambridge, 2001.

M. Eck, Least squares degree reduction of Bézier curves, Computer-Aided Design 27 (1995) 845–851.

G. Farin, Triangular Bernstein-Bézier patches, Comput. Aided Geom. Design 3 (1986) 83–127.

G. Farin, Curves and Surfaces for Computer-Aided Geometric Design. A Practical Guide, fifth ed., Academic Press, Boston, 2002.

Q.Q. Hu, An iterative algorithm for polynomial approximation of rational triangular Bézier surfaces, Appl. Math. Comput. 219 (2013) 9308–9316.

P. Keller, A method for indefinite integration of oscillatory and singular functions, Numer. Algor. [46]{} (2007) 219–251.

R. Koekoek, R.F. Swarttouw, The Askey scheme of hypergeometric orthogonal polynomials and its $q$-analogue, Rep. 98-17, Fac. Techn. Math. Informatics, Delft Univ. of Technology, Delft, 1998.

T.H. Koornwinder, Two-variable analogues of the classical orthogonal polynomials, in: R.A. Askey (ed.), Theory and Application of Special Functions, Academic Press, New York, 1975, pp. 435–495.

B.-G. Lee, Y. Park, J. Yoo, Application of Legendre-Bernstein basis transformations to degree elevation and degree reduction, Comput. Aided Geom. Design 19 (2002) 709–718.

S. Lewanowicz, P. Keller, P. Woźny, Constrained approximation of rational triangular Bézier surfaces by polynomial triangular Bézier surfaces. Numer. Algor. **75** (2017), in press; DOI 10.1007/s11075-016-0198-4.

S. Lewanowicz, P. Woźny, Connections between two-variable Bernstein and Jacobi polynomials on the triangle, J. Comput. Appl. Math. 197 (2006) 520–533.

S. Lewanowicz, P. Woźny, Bézier representation of the constrained dual Bernstein polynomials, Appl. Math. Comput. 218 (2011) 4580–4586.

S. Lewanowicz, P. Woźny, P. Keller, Polynomial approximation of rational Bézier curves with constraints, Numer. Algor. 59 (2012) 607–622.

L. Lu, G. Wang, Optimal degree reduction of Bézier curves with $G^2$-continuity, Comput. Aided Geom. Design 23 (2006) 673–683.

L. Lu, G. Wang, Application of Chebyshev II–Bernstein basis transformations to degree reduction of Bézier curves, J. Comput. Appl. Math. 221 (2008) 52–65.

A. Rababah, Distances with rational triangular Bézier surfaces, Appl. Math. Comp. 160 (2005) 379–386.

A. Rababah, B.-G. Lee, J. Yoo, A simple matrix form for degree reduction of Bézier curves using Chebyshev–Bernstein basis transformations, Appl. Math. Comput. 181 (2006) 310–318.

A. Rababah, M. Al-Natour, The weighted dual functionals for the univariate Bernstein basis, Appl. Math. Comput. 186 (2006) 1581–1590.

R. Sharma, Conversion of a rational polynomial in triangular Bernstein-Bézier form to polynomial in triangular Bernstein-Bézier form, Internat. J. Comput. Appl. Math. 8 (2013) 45–52.

H. Sunwoo, Matrix representation for multi-degree reduction of Bézier curves, Comput. Aided Geom. Design 22 (2005) 261–273.

H. Sunwoo, N. Lee, A unified matrix representation for degree reduction of Bézier curves, Comput. Aided Geom. Design 21 (2004) 151–164.

M.V. Tratnik, Some multivariate orthogonal polynomials of the Askey tableau – discrete families, J. Math. Physics 32 (1991) 2337–2342.

P. Woźny, S. Lewanowicz, Multi-degree reduction of Bézier curves with constraints, using dual Bernstein basis polynomials, Comput. Aided Geom. Design 26 (2009) 566–579.

P. Woźny, S. Lewanowicz, Constrained multi-degree reduction of triangular Bézier surfaces using dual Bernstein polynomials, J. Comput. Appl. Math. 235 (2010) 785–804.

[^1]: Corresponding author *Email addresses*: Stanislaw.Lewanowicz@cs.uni.wroc.pl (Stanis³aw Lewanowicz), Pawel.Keller@mini.pw.edu.pl (Pawe³ Keller), Pawel.Wozny@cs.uni.wroc.pl (Pawe³ WoŸny)

[^2]: Notice that for a specific choice of parameters, the rational Bézier surface reduces to a polynomial Bézier surface, so that the above approximation problem, with an additional assumption, is actually the polynomial degree reduction problem discussed in [@WL10]. In the method proposed there, dual bivariate Bernstein polynomials also play the basic role. It should be stressed, however, that the present approach is substantially different and much more efficient than the one proposed in [@WL10].
