void bz_patch_compute_bbox ( bz_patch_t *P, interval_array_t *B );
  /* 
    Stores in each element of the array {B} an interval that is guaranteed to
    contain the corresponding element of {P(x)}, when the argument point {x}
    ranges over the unit cube {U^P.d}. The array {B} must have the {P.r} indices
    and size vector {bz_patch_get_rsizes(P)}. */

void bz_patch_eval_box ( bz_patch_t *P, interval_t B[], bz_patch_t *f );
  /* 
    ??? Check this ??? 
    
    Stores in {*f} a Bézier patch that describes the restriction
    of patch {P} to of the given axis-aligned box {B[0..d-1]},
    reparametrized to the unit cube {U^d} of {R^d}. The patch
    {*f} must be allocated by the caller, with the same dimensions
    and degree as {P}. */

typedef void bz_patch_visit_t(bz_patch_t *P, interval_t B[]);
  /* 
    ??? Check this ??? 
    
    Type of a sub-patch visiting function. It receives a patch {P} and
    its corresponding domain box {B[0..P.d-1]}. The box {B}
    may be null. */

void bz_patch_enum_faces
  ( bz_patch_t *P, 
    interval_t B[], 
    bz_patch_ddim_t df
  );
  /* 
    ??? Check this ??? 
    
    Enumerates all faces of {P} with dimension {d}, in order of
    increasing index, and the corresponding faces of the box
    {B[0..df-1]}. Each face will be a patch {bf} with the same
    domain dimension as {P}, but with degree 0 (constant) along each
    axis perpendicular to the face. Each bace of {B} is an array of
    {P.d} intervals {boxf[0..P.d-1]} which are trivial along those
    same axes. Calls {proc(bf, boxt)} for each face {bf} and its
    corrsponding subdomain {boxf}. The {B} may be null, in which
    case each {boxf} will be null. */
 
double bz_patch_try_flatten_face
  ( bz_patch_t *P, 
    box_signed_dir_t dir[], 
    double tol
  );
  /* 
    ??? Check this ??? 
    
    Checks a face {f} of patch {P} described by {dir[]}; if it can be
    approximated by a single multiaffine (degree 1) patch, with error
    at most {tol} in every coordinate, makes it flat. The face is
    selected by its signature vector {dir[i]}, as explained above. In
    any case, returns the final maximum deviation of face {f} from
    flatness.

    The decision to flatten a face depends only on the Bézier coefficients
    of that face.  Thus, if two patches {b1,b2} share a face, they
    will continue to do so after being flattened by this procedure. */ 

void bz_patch_multiaffine_approx ( bz_patch_t *P, bz_patch_t *T );
  /* 
    ??? Check this ??? 
    
    Stores in {*T} a multiaffine approximation (i.e., a Bézier patch
    of degree {1}) for a given Bézier patch {*P} of arbitrary degree.
    The patch {T} must be allocated by the caller, with the same
    dimensions as {P} and degree 1 in each domain axis.
    
    The approximation has the same corners as the original; thus, if
    two cells share an entire face, their approximations will also
    share that face. (However, that is not necessarily true of cells
    that share only part of a face.) */

double bz_patch_multiaffine_error ( bz_patch_t *P, bz_patch_t *T );
  /* 
    ??? Check this ??? 
    
    Compares the general Bézier patch {P} with the multiaffine
    (degree 1) patch {T}, and returns an upper bound to the difference
    {|P(x) - T(x)|}, for any {x} in the unit cube {U}. Requires 
    {T->d == P->d} and {T->r == P->r}. */

void bz_patch_affine_approx ( bz_patch_t *P, double c[], double D[], double *err );
  /* 
    ??? Check this ??? 
    
    Stores in {c[0..d-1]} and {D[0..d^2-1]} an affine approximation
    {h(r)} for the Bézier patch {*P}. Also stores in {err} an upper
    bound for the approximation error {|P(x) - h(x)|}, for any {x} in
    {U}.
    
    In this approximation, an argument point {x} is mapped to point
    {h(x) = c + 2*(x-d)*D}, where {d = (1/2,.. 1/2)} is the center of
    the unit cube {U}. Note that the approximation may not preserve
    the corners of {P}.
    
    The matrix {D} is stored linearized by rows. Both {c} and {D} must
    be allocated by the caller. */

void bz_patch_diff ( bz_patch_t *P, ix_axis_t i, bz_patch_t *T );
  /* 
    ??? Check this ??? 
    
    Stores in {*T} the derivative of the Bézier patch {*P} along
    coordinate axis {i}; namely, a Bézier patch whose value
    {T(x)} is an array with the same shape and size as {P(x)},
    such that every elements is the derivative of the corresponding
    element of {P(x)} with respect to argument coordinate {x[i]}.
    
    The patch {*T} must be allocated by the caller, with the same number of
    indices as {P}. Its degree vector {tg} must be equal to {P}'s degree vector
    {Pg}, except that {tg[i]} must be {max(0, Pg[i]-1)}. */

void bz_patch_grad ( bz_patch_t *P, bz_patch_t *T );
  /* 
    ??? Check this ??? 
    
    Stores in {*T} the gradient of the Bézier patch {*P}; namely, a Bézier patch from {R^d} to
    {R^{r}}, such that {T(x)[i]} is the derivative of {P(x)[i]} with
    respect to {x[i]}.
    
    The patch {*T} must be allocated by the caller, with the same
    dimensions as {P}. Its degree {T->g} must be equal to {P->g},
    except that {T->g[i] = max(0, P->g[i]-1)}. */

void bz_patch_invert
  ( double *p,
    bz_patch_t *P,
    double tol,
    double x[]
  );
  /* Given a Bézier patch {P} and a point {p[0..d-1]} of {R^d},
    the procedure returns the coordinates {x[0..d-1]} 
    of the point {P^{-1}(p)}.  
    
    More precisely, the procedure finds a point {x} in {R^d} such that
    {|P(x) - p| \leq tol}. Note that {x} may lie outside the reference
    cube {U}, or even be undefined, if {p} lies outside the actual
    cell {P(U)}. */

void bz_patch_split
  ( bz_patch_t *P, 
    ix_axis_t a,
    double ratio,
    bz_patch_t *bLO,
    bz_patch_t *bHI
  );
  /* Returns descriptions of the two halves of a Bézier patch {P},
    that result from bissecting the unit {d}-cube {U} by a plane
    perpendicular to axis {a}, and reparametrizing each piece over the
    whole cube.
    
    The parameter {ratio} defines the position on the splitting plane:
    {0} means the face of the cube facing towards {-oo}, {1} means the
    face facing towards {+oo}, {0.5} means exact bisection. */ 

void bz_patch_enum_sub
  ( bz_patch_t *P, 
    interval_t B[],
    bz_patch_visit_t *proc,
    int minRank,
    int maxRank,
    double tol
  );
  /* 
    Enumerates a set of sub-patches of a Bézier patch {bz} with domain
    box {P[0..P.d-1]}. 
    
    The Bézier patch {bz} must be non-null NULL.
    
    The domain rectangle {P[0..P.d-1]} is subdivided recursively in
    half along each axis in turn, cyclically, at least {minRank}
    times. The subdivision continues until {maxRank} splits, or until
    the values of the patch, restricted to the sub-box, deviate less
    than {tol} from the multiaffine patch with the same corners. The
    function {proc} is called for each leaf (unsplit) sub-box. */

bz_patch_t bz_patch_from_box ( bz_patch_ddim_t d, interval_t B[] );
  /* Creates a Bézier patch of uniform degree 1 (multi-linear)
    with domain dimension {d} and range dimension {d},
    whose range is the box {B[0..d-1]}. */