/* Last edited on 2007-02-16 00:09:02 by stolfi */

/* ---------------------------------------------------------------------- */

typedef struct ia_butfly_t 
  { Float xmd; 
    Interval yxlo, yxmd, yxhi;
  } ia_butfly_t;
  /* Describes a butterfly-shaped region of the XY plane.
    The region consists of two trapezoids with vertical 
    bases, spanning a X interval {[xr} provided by the context. 
    The first trapezoid spans the X interval {[xr.lo _ xmd}},
    with vertices

      {(xr.lo, yxlo.lo)}, {(xr.lo, yxlo.hi)}
      {(xmd, yxmd.hi)},   {(xmd, yxmd.lo)}
      
    The second one spans {[xmd _ xr.hi]}, with vertices
    
      {(xmd, yxmd.lo)},   {(xmd, yxmd.hi)}
      {(xr.hi, yxhi.hi)}, {(xr.hi, yxhi.lo)}
      
    The middle abscissa {xmd} must belong to {xr}.
    If {xmd=xr.lo}, the firt trapezoid is empty;
    if {xmd=xr.hi}, the second one is empty. */

/* ---------------------------------------------------------------------- */
            /* Set {a_this} to a trivial butterfly: */
            /* Lo trapezoid is degenerate: */
            Float xlo = x_this.lo, xhi = x_this.hi;
            a_this->tp[0].x = (Interval) { xlo, xlo };
            a_this->tp[0].yxlo = a_this->tp[0].yxhi = y_this;
            /* Hi trapezoid is just a box: */
            a_this->tp[1].x = (Interval) { xlo, xhi };
            a_this->tp[1].yxlo = a_this->tp[1].yxhi = y_this;
            
/* ---------------------------------------------------------------------- */

        Float mrhi_this = zf_compute_max_root_hi (x_this.lo, epsilon, delta);

        k_this = zf_classify_interval(&x_this, &y_this, mrhi_this);

/* ---------------------------------------------------------------------- */

ia_butfly_t zf_butterfly_from_ia_diff(Interval *xr, Float xm, Interval *ymr, Interval *dyr)
  { ia_butfly_t a;
    
    a.xmd = xm;
    a.yxmd = *ymr;

    if (ia_is_full(dyr))
      { a.yxlo = a.yxmd = a.yxhi = ia_full(); }
    else
      { ROUND_UP;
        double dxlo = a.xmd - xr->lo;
        double dxhi = xr->hi - a.xmd;
        a.yxlo.hi = a.yxmd.hi + dyr->lo * (-dxlo);
        a.yxhi.hi = a.yxmd.hi + dyr->hi * (+dxhi);
        ROUND_DOWN;
        a.yxlo.lo = a.yxmd.lo + dyr->hi * (-dxlo);
        a.yxhi.lo = a.yxmd.lo + dyr->lo * (+dxhi);
      }
    return a;
  }

/* ---------------------------------------------------------------------- */

ia_butfly_t zf_butterfly_from_trapezoid(Interval *xr, Interval *xw, Interval *yxlor, Interval *yxhir)
  { ia_butfly_t a;

    demand((xw->lo <= tr->x.lo) && (xw->hi >= tr->x.hi), "X ranges should be nested");

    a.yxlo = *yxlor;
    a.yxhi = *yxhir;

    /* Adjust the vertical intervals {a.yxlo,a.yxhi} as needed to
      account for the clipping of {H} to the X-range {xr}: */

    if (a.yxlo.hi > a.yxhi.hi)
      { a.yxhi.hi = flt_interpolate_hi(tr->x.lo, a.yxlo.hi, xw->hi, a.yxhi.hi, tr->x.hi); }
    else if (a.yxlo.hi < a.yxhi.hi)
      { a.yxlo.hi = flt_interpolate_hi(xw->lo, a.yxlo.hi, tr->x.hi, a.yxhi.hi, tr->x.lo); }
    else
      { /* {a.yxlo.hi, a.yxhi.hi} are  OK */ }

    if (a.yxlo.lo < a.yxhi.lo)
      { a.yxhi.lo = flt_interpolate_lo(tr->x.lo, a.yxlo.lo, xw->hi, a.yxhi.lo, tr->x.hi); }
    else if (a.yxlo.lo > a.yxhi.lo)
      { a.yxlo.lo = flt_interpolate_lo(xw->lo, a.yxlo.lo, tr->x.hi, a.yxhi.lo, tr->x.lo); }
    else
      { /* {a.yxlo.lo, a.yxhi.lo} are OK */ }
      
    /* Provide an empty left trapezoid: */
    a.xmd = tr->x.lo;
    a.yxmd = a.yxlo;
    
    return a;
  }

/* ---------------------------------------------------------------------- */
        

        /* Let {RANGE(aa)} denote the implicit range of an affine form {aa},
          and {FIX(aa,k,v)} the result of substituting {v} for {e_k} in {aa},
          both computed with exact (infinite-precison) arithmetic.
          Let {e_i} be the single noise variable in {xf}.

          Since {xf} depends on only one epsilon {e_i}, the joint
          range {H} of {xf} and {yf} is actually a trapezoid with
          vertical bases. The projection of {H} on the {x}-axis is
          {XF=RANGE(xf)}, and the two bases of {H} are {YXLO =
          RANGE(FIX(yf,i,-1))} and {YXHI = RANGE(FIX(yf,i,+1))}.

          In practice we cannot compute {H}, but only conservative enclosures 
          {xfr} for {XF}, {yxlor} for {YXLO}, and {yxhir} for {YXHI}. */
        
        Interval xfr = aa_implicit_range(xf); /* Conservative enclosure of {XF}. */

        a = ia_trapez_from_trapez(xv, &xfr, &yxlor, &yxhir);

/* ---------------------------------------------------------------------- */

/* ---------------------------------------------------------------------- */

/* ---------------------------------------------------------------------- */