In applications, it is well known that high smoothness, small support and high vanishing moments are the three most important properties of a biorthogonal wavelet. In this paper, we shall investigate the mutual relations among these three properties. A characterization of $L_p\, (1\le p \le \infty)$ smoothness of multivariate refinable functions is presented. It is well known that there is a close relation between a fundamental refinable function and a biorthogonal wavelet. We shall demonstrate that any fundamental refinable function, whose mask is supported on $[1-2r,2r-1]^s$ for some positive integer $r$ and satisfies the sum rules of optimal order $2r$, has $L_p$ smoothness not exceeding that of the univariate fundamental refinable function with the mask $b_r$. Here the sequence $b_r$ on $\bZ$ is the unique univariate interpolatory refinement mask which is supported on $[1-2r,2r-1]$ and satisfies the sum rules of order $2r$. Based on a similar idea, we shall prove that any orthogonal scaling function, whose mask is supported on $[0,2r-1]^s$ for some positive integer $r$ and satisfies the sum rules of optimal order $r$, has $L_p$ smoothness not exceeding that of the univariate Daubechies orthogonal scaling function whose mask is supported on $[0,2r-1]$. We also demonstrate that a similar result holds true for biorthogonal wavelets. Examples are provided to illustrate the general theory. Finally, a general CBC (Construction By Cosets) algorithm is presented to construct all the dual refinement masks of any given interpolatory refinement mask with the dual masks satisfying arbitrary order of sum rules. Thus, for any scaling function which is fundamental, this algorithm can be employed to generate a dual scaling function with arbitrary approximation order. This CBC algorithm can be easily implemented. As a particular application of the general CBC algorithm, a TCBC (Triangle Construction By Cosets) algorithm is proposed. For any positive integer $k$ and any interpolatory refinement mask $a$ such that $a$ is symmetric about all the coordinate axes, such TCBC algorithm provides us a dual mask of $a$ such that the dual mask satisfies the sum rules of order $2k$ and is also symmetric about all the coordinate axes. As an application of this TCBC algorithm, a family of optimal bivariate biorthogonal wavelets is presented with the scaling function being a spline function.