#ifndef zf2_H
#define zf2_H

/* General interval-based zero finder for a 2-argument function. */
/* Last edited on 2017-01-02 12:26:02 by jstolfi */

#define _GNU_SOURCE
#include <stdio.h>
#include <bool.h>
#include <flt.h>
#include <ia_butfly.h>
#include <ia.h>

/* 2D ROOT FINDING

  The procedures {zf2_enum_zeros_grid} and {zf2_enum_zeros_dyadic} find
  the root set of a function {z=F(x,y)} for {x} in {xd} and {y} in {yd},
  given a procedure {eval} that bounds the graph of {F} on arbitrary
  sub-boxes of {xd\times yd}.

  More precisely, each procedure decomposes the box (axis-aligned rectangle)
  {xd\times yd} into the union of zero or more /cells/ (sub-boxes)
  of non-zero measure, with pairwise disjoint interiors.  
  
  ENCLOSING PRISMS
  
  Each sub-box is further divided into four triangles with a common
  vertex at the center. 
  
  For each triangle {T}, the procedures {zf2_enum_zeros_grid} and
  {zf2_enum_zeros_dyadic} use the client-given {eval} procedure to
  obtain an /enclosing prism/ {P}.
  
  The prism is defined by the three corners {p[i]} of {T} on the {xy}
  plane, and three intervals {zr[i]} of {z}-coorrdinates. The {xy}
  projection of {P} is the triangle {T}. The vertical edge of the
  prism projecting onto some vertex {p[i]} of {T} spans a given range
  {zr[i]} of {z}-coordinates.
  
  This prism should enclose the graph {z=F(x,y)} where {(x,y)} is in
  {T} (including the boundary) and {F(x,y)} is defined.
  
  If the three ranges {zr[i]} are empty, then the prism is 
  an empty set.  In this case, {F} must be undefined in
  the whole triangle {T}.  Otherwise all three ranges must be
  non-empty.  Note that, in this case, the function {F} may still be
  undefined in part of all of {T}.
  
  TRIANGLE KIND
  
  Each triangle is classified as either "root", "positive",
  "negative", "undefined", or "indeterminate".

  A triangle {t} is classified "positive" or "negative" if {F} is
  known to have that sign in the whole triangle. It is classified as
  "undefined" if the {eval} routine says so (see below). It is
  classified as "root" if {F} is guaranteed to have at least one zero
  in in the box. Finally, it is classified as "indeterminate" if the procedure
  cannot decide for one of the above.

  Each call to {report} is given one bounding prism generated 
  by the procedure, as explained under {zf2_report_func_t}.

  Note that since adjacent sub-boxes share a segment, 
  between a "positive" and a "negative" sub-box there must be
  at least one "root" or "undefined" sub-box.

  If any call to {report} returns TRUE, {zf2_enum_zeros} exits
  immediately with {true}.  If all calls to {report} return
  {false}, {zf2_enum_zeros} returns {false}.

  The function {F} must be coded as a routine {eval} with the
  properties explained under {zf2_eval_func_t}. The arguments {xr,yr}
  given to {eval} define always a sub-box of {xd,yd}. */

typedef enum
  { zf2_kind_undefined = 0,    /* {F} undefined in whole box. */
    zf2_kind_positive = 1,     /* {F} positive in whole box. */
    zf2_kind_negative = 2,     /* {F} negative in whole box. */
    zf2_kind_root = 3,         /* {F} zero in whole box, or box is too small. */
    zf2_kind_indeterminate = 4 /* Procedure cannot decide {F}'s kind as one of the above. */
  } zf2_kind_t;
  /* Classification of a sub-box of the domain. */
  
typedef void zf2_eval_func_t(Interval *xr, Interval *yr, Float xp[], Float yp[], Interval zr[]);
  /* Type of a procedure that can be passed as the {eval} argument
    of {zf2_enum_zeros}, to evaluate the target function {F} on a
    rectangular sub-box {xr \times yr} of the domain. 
    
    Such a procedure should return five points {p[i]=(xp[i],yp[i])} 
    and five corresponding intervals {zr[i]}, for {i} in {0..4} that 
    together define a guaranteed enclosure for the graph of the function {F}
    within that sub-box.
    
    Specifically, the five points will be the approximate center
    {p[0]} of the box {xr\times yr}, and the four corners {p[1..4]} of
    that box, in counterclockwise order, starting with the lowest
    corner. These points will be interpreted as the vertices of four
    triangles {(p[i],p[j],p[k])}, where {(i,j,k)} is {(0,1,2)},
    {(0,2,3)}, {(0,3,4)}, and {(0,4,1)}.
    
    The assumed enclosure of {F} consists of four triangular prisms
    with slanted bases and vertical sides, whose {xy} projections are
    those four triangles. The vertical side that projects on point
    {p[i]} is assumed to span the interval {zr[i]} of {z}-coordinates.
    Those four prisms should contain the graph of {F(x,y)} for all
    points {(x,y)} such that {x} is in {*xr}, {y} is in {*yr}, and
    {F(x,y)} is defined.

    The routine {eval} may set all five {zr} intervals to empty if the
    target function {F} is undefined in the whole box
    {*xr\times*yr}.  Otherwise it should set them all to 
    non-empty intervals.  */

typedef bool_t zf2_report_func_t(Float xp[], Float yp[], Interval zp[], zf2_kind_t kind);
  /* Type of a procedure that can be passed as the {report} argument
    of {zf2_enum_zeros}, to process each enclosing prism 
    returned by the {eval} function (see {zf2_eval_func_t}).
    
    The corners of the {xy} projection of the prism are the 
    points {p[i]=(xp[i],yp[i])} for {i} in {0..2}.  
    
    .  The {xy}the intervals that {xr,yr} that define
  the sub-box, its kind {k}, and an interval {zr} that contains all the
  values {F(x)} where {x} is in {xr}, {y} is in {yr}, and {F(x,y)} is
  defined. For "undefined" intevals, the range {zr} will be empty. For
  "positive" and "negative" intervals, the range {zr} may extend all the
  way to zero, even though the function {F} is guaranteed not to be zero
  in that interval.*/

bool_t zf2_enum_zeros_grid
  ( zf2_eval_func_t *eval,
    zf2_report_func_t *report,
    Interval xd,
    Interval yd,
    int nx,
    int ny
  );
  /* This procedure enumerates the zero set of {F} by dividing the
    domain {xd\times yd} into a regular grid of {nx\times ny} cells
    (sub-boxes) with uniform size (apart from roundoff errors). 
    
    These sub-boxes are then scanned systematically. The {eval}
    procedure is called on each sub-box {br = xr\times yr}. The four
    prism enclosures returned by {eval} are then iclassified, and
    passed to the {report} procedure. */

bool_t zf2_enum_zeros_dyadic
  ( zf2_eval_func_t *eval,
    zf2_report_func_t *report,
    Interval xd,
    Interval yd,
    int max_split
    
  );
  /* This procedure enumerates the zero set of {F} by dividing the domain
    recursively into cells (sub-boxes), in dyadic-tree fashion. 
    
    Namely, starting with the whole domain, each sub-box {br =
    xr\times yr} is submitted to {eval}. Then, if all four prisms
    classify as "positive", "negative", "undefined", or "root", the
    four prisms are passed to {report}. Otherwise (i.e. if one or more
    of the prisms classify as "indeterminate"), the cell {br} is
    divided into two equal sub-boxes, and the procedure is recursively
    applied to each half.
    
    The split direction alternates between the {x} and {y} axes at
    each recursion level. The recursion stops after {max_split} nested
    splittings along both axes. In that case, the prism enclosures of
    the cell {br} are given to {report}, even if some of them are
    "indeterminate". */

/* HANDY ROUTINES */

zf2_kind_t zf2_classify_triangle
  ( Float xp[],
    Float yp[],
    Interval zr[]
  );
  /* Returns the kind of a triangle {T} in the {xy} plane, given its
    corners {p[i]=(xp[i],yp[i])} and the corresponding {z}-ranges
    {zr[i]}, for {i} in {0..2}, that define an enclosing prism of {F}
    in {T}. If the prism is empty, returns {zf2_kind_undefined}. */
  
#endif