/* Last edited on 2005-09-24 22:35:12 by stolfi */ ---------------------------------------------------------------------- typedef struct zf_trapezoid_t { Interval yxlo, yxhi; } zf_trapezoid_t; zf_kind_t zf_enum_zeros ( void f (Interval *xr, Interval *yr, zf_trapezoid_t *a, void *data), int report (Interval *xr, Interval *yr, zf_kind_t kind, void *data), void *data, Interval x_dom, Float epsilon, Float delta ) { int nst = 0; /* Count of stacked intervals */ zf_stack_entry_t stack[STACKSIZE]; /* Stacked intervals */ zf_stack_entry_t *top = &(stack[0])-1; /* The top-of-stack entry */ Interval x_this; /* The current interval */ Interval y_this; /* The range of defined values of F in {x_this} */ zf_trapezoid_t a_this; /* A trapezoid containing {F(x)} for {x} in {x_this} */ zf_kind_t t_this; /* The kind of {x_this}. */ int l_this; /* Generation number of {x_this} */ Interval x_prev; /* Previous interval */ Interval y_prev; /* Its range of values */ zf_kind_t t_prev; /* Its kind */ x_this = x_dom; y_this = a_this.yxlo = a_this.yxhi = ia_full(); l_this = 0; t_this = zf_kind_mixed; x_prev.lo = One; x_prev.hi = Zero; /* No previous interval yet. */ while (x_this.lo <= x_this.hi) { /* Current interval invariant: /\ ((x \in x_this) /\ (F(x) != undefined)) => /\ F(x) \in y_this /\ (t_this == zf_kind_positive) => F(x) > 0 /\ (t_this == zf_kind_negative) => F(x) < 0 /\ t_this != zf_kind_undefined /\ t_this != zf_kind_root */ /* Stack invariant: /\ top == &(stack[nst-1]) /\_{i=1}^{nst-1} stack[i].xr.hi == stack[i-1].xr.lo /\ (nst > 0 => xv.hi == stack[nst-1].xr.lo) /\_{i=1}^{nst-1} stack[i].level >= stack[i-1].level /\ (nst > 0 => l_this >= stack[nst-1].level) */ /* filter invariant: /\ (x_prev.lo <= x_prev.hi) => /\ (x_prev.hi == x_this.lo) /\ (t_prev != zf_kind_undefined) => (y_prev.lo <= y_prev.hi) /\ (x_prev.lo > x_prev.hi) => /\ (x_dom.lo == x_this.lo)) */ if (t_this == zf_kind_mixed) { /* Try to resolve the kind of {x_this} by evaluating {F} on it: */ /* Also recompute {y_this} in the process. */ Float max_root_hi; Interval yt; max_root_hi = zf_compute_max_root_hi (x_this.lo, epsilon, delta); /* Avoid evaluating {F} on excessively small intervals; */ /* if {x_this} is too small, grow it by eating up the stack: */ while ((x_this.hi < max_root_hi) && (nst > 0)) { affirm (top->xr.lo == x_this.hi, "zerofind: stack not contiguous"); y_this = ia_join(y_this, top->yr); if (max_root_hi >= top->xr.hi) { x_this.hi = top->xr.hi; nst--; top--; } else { x_this.hi = max_root_hi; top->xr.lo = x_this.hi; top->level = l_this; } } #if DEBUG fprintf(stderr, "\n"); fprintf(stderr, " x_this = "); ia_print(stderr, x_this); fprintf(stderr, "\n"); fprintf(stderr, " y_this = "); ia_print(stderr, y_this); fprintf(stderr, "\n"); fprintf(stderr, " l_this = %3d", l_this); fprintf(stderr, " bits = %6.2f\n", zf_bits_gained(x_this.hi - x_this.lo, x_dom.hi - x_dom.lo) ); #endif f(&x_this, &yt, &a_this, data); #if DEBUG fprintf(stderr, " f = "); ia_print(stderr, y_this); fprintf(stderr, "\n"); fprintf(stderr, " f.yxlo = "); ia_print(stderr, a_this.yxlo); fprintf(stderr, "\n"); fprintf(stderr, " f.yxhi = "); ia_print(stderr, a_this.yxhi); fprintf(stderr, "\n"); fprintf(stderr, "\n"); #endif yt = ia_meet(yt, ia_join(a_this.yxlo, a_this.yxhi)); y_this = ia_meet(y_this, yt); t_this = zf_classify_interval(&x_this, &y_this, max_root_hi); } if (t_this != zf_kind_mixed) { /* Output this interval */ if (x_prev.lo <= x_prev.hi) { if (t_prev == t_this) { /* Merge with previous interval: */ x_prev.hi = x_this.hi; if (t_prev != zf_kind_undefined) { if (y_this.lo < y_prev.lo) y_prev.lo = y_this.lo; if (y_this.hi > y_prev.hi) y_prev.hi = y_this.hi; } } else { /* Report previous interval: */ int stop = report (&x_prev, &y_prev, t_prev, data); if (stop) return (t_this); x_prev = x_this; y_prev = y_this; t_prev = t_this; } } else { /* Save this interval: */ x_prev = x_this; y_prev = y_this; t_prev = t_this; } } else { /* Split x_this into smaller pieces, and stack them: */ zf_split_trapezoid ( &x_this, &y_this, &a_this.yxlo, &a_this.yxhi, l_this, &top, &nst ); } /* Unstack the next interval: */ if (nst == 0) { /* All done */ x_this.lo = One; x_this.hi = Zero; } else { x_this = top->xr; y_this = top->yr; /* Trivial (box-like) trapezoid: */ a_this.yxlo = a_this.yxhi = y_this; l_this = top->level; t_this = top->kind; nst--; top--; } } /* Flush last sub-interval: */ if (x_prev.lo <= x_prev.hi) { (void) report(&x_prev, &y_prev, t_prev, data); } return(zf_kind_undefined); } ---------------------------------------------------------------------- Interval ia_div_affine( Interval x, Float zeta, Float gamma, Float delta ) { Interval z; delta = FABS(delta); if (ia_ISFULL(x) || (FABS(zeta) <= Zero) || (FABS(gamma) >= PlusInfinity) || (delta >= PlusInfinity) ) return (ia_FULL); if (zeta < Zero) { ROUND_DOWN; z.lo = x.hi / zeta + gamma - delta; ROUND_UP; z.hi = x.lo / zeta + gamma + delta; } else { ROUND_DOWN; z.lo = x.lo / zeta + gamma - delta; ROUND_UP; z.hi = x.hi / zeta + gamma + delta; } ia_NORMFULL(z); return (z); } zerofind.ho:: zerofind.h ia.ho $(INC)/flt.ho zerofind.o:: zerofind.c zerofind.ho ia.ho $(INC)/flt.ho $(INC)/js.ho IntervalPair zerofind_join_trapezoids( Interval xv1, IntervalPair yv1, Interval xv2, IntervalPair yv2 ); /* Given two trapezoids $xvi$, $yvi$ for consecutive segments of the graph of a curve $f$, returns a single trapezoid for the same curve over the union of $xv1$ and $xv2$. */ IntervalPair zerofind_join_trapezoids( Interval xv1, IntervalPair yv1, Interval xv2, IntervalPair yv2 ) { IntervalPair yvr; /* Could be improved... */ yvr.a.lo = FMIN(FMIN(yv1.a.lo, yv1.b.lo), FMIN(yv2.a.lo, yv2.b.lo)); yvr.a.hi = FMAX(FMAX(yv1.a.hi, yv1.b.hi), FMAX(yv2.a.hi, yv2.b.hi)); yvr.b = yvr.a; return(yvr); }