/* From irtinter.c of 94-02-10 00:00 */

void irt_check_seg_eval_aa(
    pcode_proc_t *proc,   /* The function's pseudo-code */
    h3_point_t *org,      /* Start of ray */
    h3_point_t *dst,      /* End of ray */
    AAP aat,              /* The affine form of the $t$ parameter */
    IntervalPair *ip,     /* Bounding trapezoid for f(aat) */
    AAP f,                /* The affine form for f(aat)  */
    AAP *regs,            /* Evaluation registers.  */
    AAP *stack            /* Evaluation stack  */
  )
  { 
    #define NCHECK 35
    
    MemP frame = aa_top();
    AATerm eps[1];
    int i, j;
    
    affirm((aat->nterms == 1), "irt_check_seg_eval_aa: bad aat");
    eps[0].id = ((AATermP)(aat + 1))->id;
    
    fprintf(stderr, "\n"); 
    for (j=-NCHECK; j<=NCHECK; j++)
      { /* 
          The following assumes that $t(e)$ is the exact value of $t$
          corresponding to $eps_k = e$, where $eps_k$ is the (unique) 
          noise var that appears in $aat$.
        */
        Float e; /* A sample value for $eps_k$ */
        Interval f_cmp; /* Range for $f(t(e))$ computed as f(aa_fix_eps(aat, e)) */
        Interval f_fix; /* Range for $f(t(e))$ computed as aa_fix_eps(f(aat), e) */
        Interval f_trp; /* Range for $f(t(e))$ obtained by interpolating $ip$ */
        AAP ct; /* An affine form for $t(e)$ */
        AAP cr; /* An affine form for $1 - t(e)$. */
        
        ROUND_NEAR; e = ((Float) j)/((Float) NCHECK);
        eps[0].coef = e;
        
        ct = aa_fix_eps(aat, 1, eps);
        cr = aa_affine (ct, -1.0, 1.0, 1.0, 0.0);
        f_fix = aa_range(aa_fix_eps(f, 1, eps));
        for (i=0; i<4; i++) 
          { regs[i] = aa_affine_2 (cr, org->c[i], ct, dst->c[i], 1.0, 0.0, 0.0); }
        aa_eval (regs, stack, proc->code);
        f_cmp = aa_range(stack[0]);
        { Float alo, ahi, blo, bhi;
          if (ip->a.lo >= Zero)
            { ROUND_DOWN; alo = ((One - e)/Two)*ip->a.lo;}
          else
            { ROUND_UP;   alo = -(((One - e)/Two)*(-ip->a.lo)); }
          if (ip->b.lo >= Zero)
            { ROUND_DOWN; blo = ((One + e)/Two)*ip->b.lo;}
          else
            { ROUND_UP;   blo = -(((One + e)/Two)*(-ip->b.lo)); }
          
          if (ip->a.hi >= Zero)
            { ROUND_UP;   ahi = ((One - e)/Two)*ip->a.hi;}
          else
            { ROUND_DOWN; ahi = -(((One - e)/Two)*(-ip->a.hi)); }
          if (ip->b.hi >= Zero)
            { ROUND_UP;   bhi = ((One + e)/Two)*ip->b.hi;}
          else
            { ROUND_DOWN; bhi = -(((One + e)/Two)*(-ip->b.hi)); }
            
          ROUND_DOWN; f_trp.lo = alo + blo;
          ROUND_UP;   f_trp.hi = ahi + bhi;
	}

        fprintf(stderr, "      ");
        fprintf(stderr, "  eps = "); ROUND_NEAR; flt_print(stderr, e);
        fprintf(stderr, "  t = "); ia_print(stderr, aa_range(ct));
        fprintf(stderr, "\n");
        fprintf(stderr, "      ");
        fprintf(stderr, "  f_cmp = "); ia_print(stderr, f_cmp);
        fprintf(stderr, "\n");
        fprintf(stderr, "      ");
        fprintf(stderr, "  f_fix = "); ia_print(stderr, f_fix);
        fprintf(stderr, "\n");
        fprintf(stderr, "      ");
        fprintf(stderr, "  f_trp = "); ia_print(stderr, f_trp);
        fprintf(stderr, "\n");
        fprintf(stderr, "\n");
        if ((f_cmp.lo > f_fix.hi) || (f_cmp.hi < f_fix.lo))
          { error("irt_check_seg_eval_aa: inconsistent f_cmp, f_fix"); }
        if ((f_trp.lo > f_fix.hi) || (f_trp.hi < f_fix.lo))
          { error("irt_check_seg_eval_aa: inconsistent f_trp, f_fix"); }
        if ((f_cmp.lo > f_trp.hi) || (f_cmp.hi < f_trp.lo))
          { error("irt_check_seg_eval_aa: inconsistent f_cmp, f_trp"); }
        aa_flush(frame);
      }
    fprintf(stderr, "\n"); 
      
    #undef NCHECK
  }

void irt_trapezoid_from_aa (AAP aaf, AAP aat, IntervalPair *ip)
  { 
    if (aat->nterms == 0)
      { ip->a = aa_range(aaf); 
        ip->b = ip->a; 
      }
    else
      { AATerm eps[1];
        eps[0].id = ((AATermP)(aat + 1))->id;
        eps[0].coef = -One;
        ip->a = aa_range(aa_fix_eps(aaf, 1, eps));
        eps[0].coef = One;
        ip->b = aa_range(aa_fix_eps(aaf, 1, eps));
      }
  }

/* From irtinter.c of 94-02-11 18:00 */

void irt_check_reported_root(Interval *xv, Interval *yv, zf_type tv, void *data);
  /* 
    Does several consistency checks on a root reported by $zerofind$.
    The latter implicitly asserts that $f(x)$ is in $yv$ for $x$ in $xv$, and that
    $tv$ is the type of $yv$.  This procedure computes $f(xv)$,
    and checks whether the answer is consistent with those assertions. */


int irt_process_root (Interval *xv, Interval *yv, zf_type tv, void *data)
  {
    cmp_int_data *cd = (cmp_int_data *) data;
    int sv;

    affirm(cd->hit->lo > cd->hit->hi, "irt_process_root: zerofind didn't stop!");

    if (tv == zf_type_positive)
      { sv = +1; }
    else if (tv == zf_type_negative)
      { sv = -1; }
    else if (tv == zf_type_root)
      { sv = 0; }
    else
      { error("irt_process_root: invalid interval type"); }
      
    /* Consistency check:*/
    if ((tv == zf_type_positive) || (tv == zf_type_negative))
      { Interval sv;
        irt_check_reported_root(xv, yv, tv, data); 
        sv = (Interval){xv->lo, xv->lo};
        irt_check_reported_root(&sv, yv, tv, data); 
        sv = (Interval){xv->hi, xv->hi};
        irt_check_reported_root(&sv, yv, tv, data); 
      }
      
    if (sv != 0)
      { *(cd->slo) = sv; return(0); }
    else if ((xv->lo > Zero) && (xv->hi < One))
      { *(cd->hit) = *xv; return(1); }
    else 
      { return(0); }
  }

void irt_check_reported_root(Interval *xv, Interval *yv, zf_type tv, void *data)
  { IntervalPair fp = irt_eval_f_on_sub_segment (xv, data);
    int ierr = 0;
    if ((fp.a.hi < yv->lo) || (fp.a.lo > yv->hi) || (fp.b.hi < yv->lo) || (fp.b.lo > yv->hi))
      { fprintf(stderr, "*** irt_process_root: inconsistent yv, f(xv)\n"); ierr = 1; }
    if (
      ((tv == zf_type_positive) && (yv->hi <= Zero)) ||
      ((tv == zf_type_negative) && (yv->lo >= Zero)) ) 
      { fprintf(stderr, "*** irt_process_root: inconsistent tv, yv\n"); ierr = 1; }
    if (
      ((tv == zf_type_positive) && (fp.a.hi <= Zero) && (fp.b.hi <= Zero)) ||
      ((tv == zf_type_negative) && (fp.a.lo >= Zero) && (fp.b.lo >= Zero)) ) 
      { fprintf(stderr, "*** irt_process_root: inconsistent tv, f(xv)\n"); ierr = 1; }
    if (ierr) { error("aborted"); }
  }