/* See {zf2.h} */
/* Last edited on 2020-10-01 20:43:43 by jstolfi */

#define _GNU_SOURCE
#include <stdio.h>
#include <assert.h>
#include <math.h>

#include <bool.h>
#include <affirm.h>
#include <ia.h>
#include <flt.h>

#include <zf2.h>

/*** PROTOTYPES FOR INTERNAL PROCEDURES ***/

float zf2_bits_gained(Float curw, Float orgw);
  /* Number of information bits gained when reducing an interval
    of width {orgw} to one of width {curw}. */
  
double zf2_log2 (double x);  

/*** IMPLEMENTATIONS ***/

#define DEBUG (0)

bool_t zf2_enum_zeros_grid
  ( zf2_eval_func_t *eval,
    zf2_report_func_t *report,
    Interval xd,
    Interval yd,
    int nx,
    int ny
  )
  {
    /* Loop on grid rows: */
    Interval yr = (Interval){ .lo = yd.lo, .hi = NAN }; 
    for (int iy = 0; iy < ny; iy++)
      { /* Scan a cell row that starts at {yr.lo}. */
        /* Set the upper limit of the {y} interval: */
        ROUND_NEAR;
        yr.hi = (Float)(yd.lo + ((double)iy + 1)/((double)ny));
        /* Loop on grid cells in the row: */
        Interval xr = (Interval){ .lo = xd.lo, .hi = NAN };
        for (int ix = 0; ix < nx; ix++)
          { /* Set the upper limit of the {x} interval: */
            ROUND_UP;
            xr.hi = (Float)(xd.lo + ((double)ix + 1)/((double)nx));
            /* Evaluate the function in the cell: */
            Interval zr;
            eval(&xr, &yr, &zr);
            
            /* Classify the cell: */
            zf2_kind_t kr = zf2_classify_cell(&xr, &yr, &zr);
            
            /* Report it: */
            bool_t stop = report(&xr, &yr, &zr, kr);
            if (stop) { return TRUE; }

            /* Advance to next cell in row: */
            xr.lo = xr.hi; xr.hi = NAN;
          }
        /* Advance to next row of cells: */
        yr.lo = yr.hi; yr.hi = NAN;
    }    
  }

bool_t zf2_enum_zeros_dyadic
  ( zf2_eval_func_t *eval,
    zf2_report_func_t *report,
    Interval xd,
    Interval yd,
    int max_split
  )
  {
    auto bool_t zf2_enum_zeros_dyadic_aux(Interval *xr, Interval *yr, int iax, int mxsp);
      /* Finds the zero set of {F} in the sub-box {xr\times yr}.
        If the box needs to be subdivided, splits it in half
        by a plane perpendicular to axis {iax}. The axis
        alternates at each recursion level.  The recursive splitting
        is applied at most {mxsp} levels on each axis. */
      
    bool_t result = zf2_enum_zeros_dyadic_aux(&xd, &yd, 0, max_split);
    return result;
    
    /* INTERNAL IMPLEMENTATIONS */
    
    bool_t zf2_enum_zeros_dyadic_aux(Interval *xr, Interval *yr, int iax, int mxsp) 
      { 
        /* Evaluate the function in the cell: */
        Interval zr;
        eval(xr, yr, &zr);
        /* Classify the cell: */
        zf2_kind_t kr = zf2_classify_cell(&xr, &yr, &zr);
        if (kr != zf2_kind_mixed)
          { 
            /* Report it: */
            bool_t stop = report(xr, yr, &zr, kr);
            if (stop) { return TRUE; }
          }
        else
          { /* Must split and recurse: */
            if (iax == 0)
              { /* Split in {x}. */
                Float xm = ia_mid(*xr);
                Interval xr0 = (Interval){ .lo = xr.lo; .hi = xm };
                Interval xr1 = (Interval){ .lo = xr.lo; .hi = xm };
                /* Next split will be in {y}, same splitting limit: */
                stop = zf2_enum_zeros_dyadic_aux(&xr0, yr, 1, mxsp);
                if (stop) { return TRUE; }
                stop = zf2_enum_zeros_dyadic_aux(&xr1, yr, 1, mxsp);
                if (stop) { return TRUE; }
              }
            else
              { /* Split in {y}: */
                Float ym = ia_mid(*yr);
                Interval yr0 = (Interval){ .lo = yr.lo; .hi = ym };
                Interval yr1 = (Interval){ .lo = yr.lo; .hi = ym };
                /* Next split will be in {x}, splitting limit is one less: */
                stop = zf2_enum_zeros_dyadic_aux(xr, &yr0, 0, mxsp-1);
                if (stop) { return TRUE; }
                stop = zf2_enum_zeros_dyadic_aux(xr, &yr1, 0, mxsp-1);
                if (stop) { return TRUE; }
              }
          }
        return FALSE;
      }  
  }
  

zf2_kind_t zf2_classify_box
  ( Interval *xr,
    Interval *yr,
    Interval *zr,
    Float epsilon,
    Float delta
  )
  { if (yr->lo > yr->hi)
      { yr->lo = One; yr->hi = Zero;
        return(zf2_kind_undefined); 
      }
    else 
      { if (yr->lo > Zero)
          { return(zf2_kind_positive); }
        else if (yr->hi < Zero)
          { return(zf2_kind_negative); }
        else if ((yr->lo == Zero) && (yr->hi == Zero))
          { return(zf2_kind_root); }
        else
          { /* Decide if it is small enough: */
            Float xm;
            ROUND_NEAR;
            xm = Half * xr->lo + Half * xr->hi;
            if 
              ( (xr->hi <= max_root_hi) ||
                (xm <= xr->lo) || (xm >= xr->hi)
              )
              { return(zf2_kind_root); }
            else
              { return(zf2_kind_mixed); }
          }
      }
  }
  
float zf2_bits_gained(Float curw, Float orgw)
  { float curbits, orgbits;
    affirm (curw <= orgw, "zf2_bits_gained: bad arguments");
    ROUND_NEAR;
    if (orgw == Zero) 
      { return(Zero); }
    else
      { orgbits = -(float)zf2_log2((double)orgw); }
    if (curw == Zero) 
      { curbits = (float)zf2_log2((double)MaxFloat); }
    else
      { curbits = -(float)zf2_log2((double)curw); }
    return(curbits - orgbits);
  }

double zf2_log2 (double x)
  { return (double)(1.4426950408889634074L * log(x)); }