void bz_patch_invert
  ( double p[],
    bz_patch_t *P,
    double tol,
    double x[]
  )
  { /* Subdivide the shape by DeCasteljau (possibly unequal), until */
    /* it can be approximated by an affine function, then invert that. */
    affirm(FALSE, "not implemented");
  }

void bz_patch_compute_bbox( bz_patch_t *P, interval_array_t *B )
  { /* Uses the fact that the cell is contained in the hull of the control pts. */
    int tsz, i, j, t;
    int m = P->d, n = P->n;
    double *c = P->c;
    /* Compute total size {tsz} of the control point array: */
    tsz = n; for (i = 0; i < m; i++) { tsz *= (P->g[i]+1); }
    /* Now find the control point range in each range coordinate {j}: */
    for (j = 0; j < n; j++)
      { double lo = INFINITY;
        double hi = -INFINITY;
        for (t = j; t < tsz; t += n)
          { double ci = c[t];
            if (ci < lo) { lo = ci; }
            if (ci > hi) { hi = ci; }
          }
        box[j] = (interval_t){{lo, hi}};
      }
  }
  
void bz_patch_multiaffine_approx(bz_patch_t *P, bz_patch_t *t)
  { /* Just get the corner control points. */
    /* Then subtract from original to bound the error. */
    int m = P->d;
    const bz_degree_t *bg = &(P->g[0]);
    const bz_degree_t *tg = &(t->g[0]);
    bool_t trivial = TRUE;
    demand(t->d == P->d, "domain dimension mismatch");
    demand(t->n == P->n, "range dimension mismatch");
    { int i;
      for (i = 0; i < m; i++)
        { demand(tg[i] == 1, "approximant is not multiaffine"); 
          if (bg[i] != 1) { trivial = FALSE; }
        }
    }
    { int n = P->n;
      int tsz = ipow(2, m);
      if (trivial)
        { int ij;
          double *tc = &(t->c[0]), *bc = &(P->c[0]); 
          /* Copy the control points literally: */
          for (ij = 0; ij < n*tsz; ij++,tc++,bc++) { (*tc) = (*bc); }
        }
      else
        { static int bstep[bz_patch_MAX_DDIM];
          int bsz, tix;
          bz_patch_compute_steps(P, bstep, &bsz);
          /* Copy the corner control points: */
          /* The indices {tix} and {bix} do not include the factor {n}. */
          for (tix = 0; tix < tsz; tix++)
            { /* Compute index {bix} s.t. {P->c[bix]} corresponds to {t->c[tix]}: */
              int bix = 0, tt = tix, i;
              for (i = m-1; i >= 0; i--) 
                { int tei = tt & 1;
                  bix += tei*bg[i]*bstep[i];
                  tt >>= 1;
                }
              { double *tc = &(t->c[n*tix]), *bc = &(P->c[n*bix]); 
                int j;
                for (j = 0; j < n; j++,tc++,bc++) { (*tc) = (*bc); }
              }
            }
        }
    }
  }
  
double bz_patch_multiaffine_error(bz_patch_t *P, bz_patch_t *t)
  { /* Just get the corner control points. */
    /* Then subtract from original to bound the error. */
    const bz_degree_t *bg = &(P->g[0]);
    const bz_degree_t *tg = &(t->g[0]);
    bool_t trivial = TRUE;
    demand(t->d == P->d, "domain dimension mismatch");
    demand(t->n == P->n, "range dimension mismatch");
    { int i;
      for (i = 0; i < P->d; i++)
        { demand(tg[i] == 1, "approximant is not multiaffine"); 
          if (bg[i] != 1) { trivial = FALSE; }
        }
    }
    if (trivial)
      { return 0.0; }
    else
      { int m = P->d;
        int n = P->n;
        static double x[bz_patch_MAX_DDIM];
        static double tx[bz_patch_MAX_DDIM];
        static int bstep[bz_patch_MAX_DDIM];
        int bsz, bix, i;
        double dmax = 0.0;
        bz_patch_compute_steps(P, bstep, &bsz);
        for (bix = 0; bix < bsz; bix++)
          { /* Compute nominal coordinates {x} of control point {P->c[bix]}: */
            int bb = bix;
            for (i = m-1; i >= 0; i--) 
              { int bnumi = bg[i] + 1, bei = bb % bnumi;
                x[i] = ((double)bei)/((double)bg[i]);
                bb /= bnumi;
              }
            /* Evaluate the multiaffine spline {t} at {x}: */
            bz_patch_eval(t, x, tx);
            /* Compare control point {P->c[bix]} with interpolated point {tx}: */
            { double *tc = tx, *bc = &(P->c[n*bix]); 
              int j;
              for (j = 0; j < n; j++,tc++,bc++) 
                { double dj = fabs((*tc) - (*bc));
                  if (dj > dmax) { dmax = dj; }
                }
            }
          }
        return dmax;
      }
  }

void bz_patch_affine_approx(bz_patch_t *P, double c[], double M[], double *err)
  { /* Compute approximate center, mean derivatives. */
    /* Then re-convert to Bézier and subtract from original to bound the error. */
    affirm(FALSE, "not implemented");
  }
  
void bz_patch_grad(bz_patch_t *P, ix_axis_t ax, bz_patch_t *t)
  { affirm(FALSE, "not implemented"); }

void bz_patch_split
  ( bz_patch_t *P, 
    ix_axis_t a, 
    double ratio,
    bz_patch_t *bLO,
    bz_patch_t *bHI
  )
  { int m = P->d, n = P->n;
    const bz_degree_t *g = &(P->g[0]);
    double *c = P->c, *cLO = bLO->c, *cHI = bHI->c;
    
    /* The DeCasteljau interpolation must act on the sequence of {g[a]+1}
      points {c'[t] = c[e[0],..e[m-1]]} where all {e[i]} stay fixed
      except {t = e[a]} which varies from {0} to {g[a]}. 
      
      Recall that the index {r} of control point {c[e[0],..e[m-1]]} 
      in the linearized control point array is
      
        {r = SUM{ G[k]*e[k] :  k = 0..d-1}}
        
      where {G[k] = PROD{g[j]+1 : j = k+1..d-1}}.
      
      Breaking this sum at the term {k = a} we get {r = u + v + w}, where
      
        {u = SUM{ G[k]*e[k] :  k = 0..a-1}}
        {v = G[a]*e[a]}
        {w = SUM{ G[k]*e[k] :  k = a+1..d-1}}
        
      Extracting out the common factors, we get
      
        {u = G[a-1] * SUM{ PROD{g[j]+1 : j = k+1..a-1}*e[k] :  k = 0..a-1}}
        
      As we vary the {e[i]} from {0} to {g[i]}, the terms {u,v,w} vary 
      from 0 (inclusive) to limits {limu,limv,limw} (exclusive)
      with steps {du,dv,dw}, where 
      
        {limw = G[a]}           {dw = 1}
        {limv = G[a-1]}         {dv = G[a] = limw}
        {limu = G[m-1]}         {du = G[a-1] = limv}
        
      Keep in mind that {r = u + v + w} is still a *point* index, and
      each point has {n} coordinates. So {r} must be multiplied by
      {n}, and added to the component index {i} which ranges
      over {0..n-1}. */
    
    int u, v, v1, w, i, k, limu, limv, limw, du, dv, dw;
    double s = 1.0 - ratio, t = ratio;
    
    /* fprintf(stderr, "+ pz_patch_split\n"); */ 
    /* Compute {limu,limv,limw} (already times {n}): */
    limw = n; for(k = a+1; k < m; k++) { limw*= (g[k]+1); }
    limv = limw * (g[a]+1);
    limu = limv; for(k = 0; k < a; k++) { limu*= (g[k]+1); }
    /* Compute the increments {du,dv,dw} (already times {n}): */
    du = limv; dv = limw; dw = n;
    
    /* Now do it: */
    /* fprintf(stderr, "= pz_patch_split cLO = %d cHI = %d\n", (int)cLO, (int)cHI); */ 
    /* fprintf(stderr, "limu = %6d  limv = %6d  limw = %6d\n", limu, limv, limw); */
    /* fprintf(stderr, "du =   %6d  dv =   %6d  dw =   %6d\n", du, dv, dw);  */
    for (u = 0; u < limu; u += du)
      { for (w = 0; w < limw; w += dw)
          { for (i = 0; i < n; i++)
              { int uwi = u + w + i;
                /* DeCasteljau algorithm along {v} index: */
                for (v = 0; v < limv; v += dv) 
                  { cHI[uwi + v] = c[uwi + v]; }
                for (v = 0; v < limv; v += dv)
                  { cLO[uwi + v] = cHI[uwi];
                    for (v1 = 0; v1 < limv - v - dv; v1 += dv)
                      { cHI[uwi + v1] = s*cHI[uwi + v1] + t*cHI[uwi + v1 + dv]; }
                  }
              }
          }
      }
    /* fprintf(stderr, "- pz_patch_split\n"); */ 
  }

double bz_patch_try_flatten_face(
    bz_patch_t *P, 
    box_signed_dir_t dir[], 
    double tol
  )
  { /* Data for Bézier patch {P}: */
    int bm = P->d;
    bz_degree_t *bg = &(P->g[0]);
    /* Data for its face {F}: */
    int tm;
    bz_degree_t tg[bz_patch_MAX_DDIM];
    double err;
    int trivial = TRUE;

    /* Get signature, dimension, and degree vector of face: */
    { int ib;
      for (ib = 0, tm = 0; ib < bm; ib++)
        { if (dir[ib] == 00) 
            { tg[tm] = bg[ib]; tm++;
              if (bg[ib] != 1) { trivial = FALSE; } 
            } 
        }
    }
    if (trivial)
      { /* Face is multi-affine, hence already flat: */
        return 0.0;
      } 
    else
      { /* Could be done more efficiently... */
        bz_patch_t t = bz_patch_new(tm, P->n, tg);
        bz_patch_t a = bz_patch_uniform_new(tm, P->n, 1);
        fprintf(stderr, "+ bz_try_flatten_face\n");
        bz_patch_debug("P = \n", P, "\n");
        bz_patch_get_face(P, dir, &t);
        bz_patch_debug("orig t = \n", &t, "\n");
        bz_patch_multiaffine_approx(&t, &a);
        bz_patch_debug("flat a = \n", &a, "\n");
        err = bz_patch_multiaffine_error(&t, &a);
        if (err <= tol)
          { /* Flatten face {f} of {P}, propagate correction to superfaces: */
            fprintf(stderr, "flattening the face:\n");
            bz_patch_raise_degree(&a, &t);
            bz_patch_debug("flat t = \n", &t, "\n");
            bz_patch_set_face(P, dir, &t);
            bz_patch_debug("modf P = \n", P, "\n");
            err = 0.0; 
          }
        free(t.c);
        free(a.c);
        fprintf(stderr, "- bz_patch_try_flatten_face\n");
        return err;
      }
  }

void bz_patch_enum_sub
  ( bz_patch_t *P, 
    interval_t box[],
    bz_patch_visit_t *proc,
    int minRank,
    int maxRank,
    double tol
  )
  {
  
    bz_patch_t t = bz_patch_uniform_new(P->d, P->n, 1);
    
    auto void recurse(bz_patch_t *P, interval_t box[], int rank);
      /* Performs the task on a sub-patch {P,box} at level {rank}. */
    
    void recurse(bz_patch_t *P, interval_t box[], int rank)
      { /* Decide whether a patch must be split or is a leaf: */
        bool_t leaf;
        if (rank < minRank)
          { leaf = FALSE; }
        else if (rank >= maxRank)
          { leaf = TRUE; }
        else
          { bz_patch_multiaffine_approx(P, &t);
            double err = bz_patch_multiaffine_error(P, &t);
            if (err <= tol)
              { leaf = TRUE; }
            else
              { leaf = FALSE; }
          }
        /* Now act on the decision: */
        if (leaf)
          { proc(P, box); }
        else
          { /* Split along axis {rank % P->d}: */
            ix_axis_t a = rank % P->d;
            bz_patch_t bLO = bz_patch_new(P->d, P->n, P->g);
            bz_patch_t bHI = bz_patch_new(P->d, P->n, P->g);
            interval_t boxLO[P->d], boxHI[P->d];
            box_split(P->d, box, a, 0.5, boxLO, boxHI);
            bz_patch_split(P, a, 0.5, &bLO, &bHI);
            recurse(&bLO, boxLO, rank+1);
            recurse(&bHI, boxHI, rank+1);
            bz_patch_free(&bLO);
            bz_patch_free(&bHI);
          }
      }
      
    recurse(P, box, 0);
    free(&t);
  }
  
void bz_patch_eval_box(bz_patch_t *P, interval_t box[], bz_patch_t *f)
  { assert(P->C.ds.na == P->d + P->r);
    affirm(FALSE, "not implemented");
  }


bz_patch_t bz_patch_from_box(bz_patch_ddim_t d, interval_t box[])
  { 
    /* fprintf(stderr, "+ bz_patch_from_box\n"); */
    bz_patch_t P = bz_patch_uniform_new(d, d, 1);
    int ix; 
    int bsz = ipow(2, d);
    for (ix = 0; ix < bsz; ix++)
      { /* Get labels {e[0..d-1]} of Bezier coeff {P.c[ix]} */
        int t = ix, i;
        for (i = 0; i < d; i++) 
          { int ei = (t % 2); t /= 2;
            P.c[d*ix + i] = box[i].end[ei];
          }
      }
    /* fprintf(stderr, "- bz_patch_from_box\n"); */
    return P;
  }