/******************************************************************************
* MvarZero.c - tools to compute zero sets of multivariates.		      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, July. 97.					      *
******************************************************************************/

#include "mvar_loc.h"
#include "misc_lib.h"

#define MVAR_NUMER_ZERO_NUM_STEPS 30

typedef struct MvarMVParamDomainStruct {
    CagdRType Min, Max;					      /* The domain. */
} MvarMVParamDomainStruct;

STATIC_DATA int
    GlblApplyNormalConeTest = FALSE;

static MvarMVParamDomainStruct *MVarCopySubdivParamDomains(MvarMVStruct *MV,
							   CagdRType t,
							   int Dir,
							   CagdBType FirstSubMV);
static void MvarGetSubdivParamDomains(MvarMVStruct *MV,
				      CagdRType *Min,
				      CagdRType *Max,
				      int Dir);
static MvarPtStruct *MvarMVsSubdivZeros(MvarMVStruct **MVs,
					MvarConstraintType *Constraints,
					int NumOfMVs,
					int NumOfZeroMVs,
					int ApplyNormalConeTest,
					CagdRType SubdivTol,
					int Depth);
static MvarPtStruct *MvarGenPtMidMvar(MvarMVStruct *MV, int SingleSol);
static MvarPtStruct *MvarMVsNumericZeros(MvarPtStruct **ZeroSet,
					 MvarMVStruct **MVs,
					 MvarConstraintType *Constraints,
					 int NumOfMVs,
					 CagdRType Neighborhood,
					 CagdRType NumericTol);
static CagdRType MvarEvalErrorL1(MvarMVStruct **MVs,
				 CagdRType *Params,
				 int NumOfMVs);

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Sets the use (or not) of the normal cone tests inside the multivariate   M
* zero set solver.							     M
*                                                                            *
* PARAMETERS:                                                                M
*   NormalConeTest:   New setting for normal cone testing usage.             M
*                                                                            *
* RETURN VALUE:                                                              M
*   int:       Old setting for normal cone testing usage.                    M
*                                                                            *
* SEE ALSO:                                                                  M
*   MvarMVsZeros                                                             M
*                                                                            *
* KEYWORDS:                                                                  M
*   MvarMVsZerosNormalConeTest                                               M
*****************************************************************************/
int MvarMVsZerosNormalConeTest(int NormalConeTest)
{
    int OldVal = GlblApplyNormalConeTest;

    GlblApplyNormalConeTest = NormalConeTest;

    return OldVal;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes the simultaneous solution of the given set of NumOfMVs          M
* constraints.  A constraints can be equality or ineqaulity as prescribed    M
* by the Constraints vector.  Only equality constraints are employed in the  M
* numerical improvment stage.						     M
*   All multivariates are assumed to be scalar and be in the same parametric M
* domain size and dimension.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   MVs:          Vector of multivariate constraints.                        M
*   Constraints:  Either an equality or an inequality type of constraint.    M
*   NumOfMVs:     Size of the MVs and Constraints vector.                    M
*   SubdivTol:    Tolerance of the subdivision process.  Tolerance is        M
*		  measured in the parametric space of the multivariates.     M
*   NumericTol:   Numeric tolerance of the numeric stage.  The numeric stage M
*		  is employed only if FABS(NumericTol) < SubdivTol.          M
*		  If NumericTol is negative, points that fail to improve     M
*		  numerically are purged away.				     M
*                                                                            *
* RETURN VALUE:                                                              M
*   MvarPtStruct *:   List of points on the solution set.  Dimension of the  M
*		      points will be the same as the dimensions of all MVs.  M
*                                                                            *
* SEE ALSO:                                                                  M
*   MvarMVsZerosNormalConeTest                                               M
*                                                                            *
* KEYWORDS:                                                                  M
*   MvarMVsZeros                                                             M
*****************************************************************************/
MvarPtStruct *MvarMVsZeros(MvarMVStruct **MVs,
			   MvarConstraintType *Constraints,
			   int NumOfMVs,
			   CagdRType SubdivTol,
			   CagdRType NumericTol)
{
    int i, l, NumOfZeroMVs, CanApplyNormalConeTest;
    CagdRType
	MinScale = 1.0;
    MvarPtStruct *Pt, *ZeroSet,
	*OutSet = NULL;
    MvarMVStruct **LclMVs;

    /* Make sure all multivariates are scalar fields. */
    for (i = 0; i < NumOfMVs; i++) {
	if (MVAR_NUM_OF_PT_COORD(MVs[i]) != 1) {
	    MVAR_FATAL_ERROR(MVAR_ERR_SCALAR_PT_EXPECTED);
	    return NULL;
	}
    }

    /* Make sure the parametric domain is identical at all multivariates. */
    for (l = 0; l < MVs[0] -> Dim; l++) {
	CagdRType Min, Max;

	MvarMVDomain(MVs[0], &Min, &Max, l);

	for (i = 1; i < NumOfMVs; i++) {
	    CagdRType Min2, Max2;

	    MvarMVDomain(MVs[i], &Min2, &Max2, l);
	    if (!APX_EQ(Min, Min2) || !APX_EQ(Max, Max2)) {
		MVAR_FATAL_ERROR(MVAR_ERR_INCONS_DOMAIN);
		return NULL;
	    }
	}
    }

    /* Make sure we have domains on all multivariates. */
    LclMVs = (MvarMVStruct **) IritMalloc(sizeof(MvarMVStruct) * NumOfMVs);
    for (i = 0; i < NumOfMVs; i++) {
	CagdRType Scale;
	STATIC_DATA CagdVType
	    Translate = { 0.0, 0.0, 0.0 };
	CagdBBoxStruct BBox;

	LclMVs[i] = MvarMVCopy(MVs[i]);
	if (MVAR_IS_BEZIER_MV(MVs[i])) {
	    int j;
	    MvarMVParamDomainStruct *Dmns;

	    LclMVs[i] -> PAux = Dmns = (MvarMVParamDomainStruct *)
		IritMalloc(sizeof(MvarMVParamDomainStruct) * MVs[i] -> Dim);
	    for (j = 0; j < MVs[i] -> Dim; j++) {
	        Dmns[j].Min = 0.0;
	        Dmns[j].Max = 1.0;
	    }
	} 
	if (MVAR_IS_RATIONAL_MV(LclMVs[i])) {
	    /* Convert P1 point type to E1, in local copy. */
#ifndef MVAR_MALLOC_STRUCT_ONCE
	    IritFree(LclMVs[i] -> Points[0]);
#endif /* MVAR_MALLOC_STRUCT_ONCE */
	    LclMVs[i] -> Points[0] = NULL;
	    LclMVs[i] -> PType = CAGD_PT_E1_TYPE;
	}

	/* For numeric stability, scale all constraints so that the maximal */
	/* coefficient's absolute value is one.				    */
	if (!CagdPointsHasPoles(LclMVs[i] -> Points,
				MVAR_CTL_MESH_LENGTH(LclMVs[i]))) {
	    MvarMVBBox(LclMVs[i], &BBox);
	    Scale = 1.0 / MAX(FABS(BBox.Min[0]), FABS(BBox.Max[0]));
	    MinScale = MIN(MinScale, Scale);
	    MvarMVTransform(LclMVs[i], Translate, Scale);
	}
    }
    NumericTol *= MinScale;

    /* Make sure all zero constraints are first, and count how many. */
    for (i = 0; i < NumOfMVs; i++) {
        if (Constraints[i] != MVAR_CNSTRNT_ZERO) {
	    int j;

	    for (j = i + 1; j < NumOfMVs; j++) {
	        if (Constraints[j] == MVAR_CNSTRNT_ZERO) {
		    /* Swap constraints i and j. */
		    SWAP(MvarConstraintType, Constraints[i], Constraints[j]);
		    SWAP(MvarMVStruct *, LclMVs[i], LclMVs[j]);
		    break;
		}
	    }
	    if (j >= NumOfMVs)
	        break;
	}
    }
    NumOfZeroMVs = i;

    /* See of normal cone test could be applied and if so prepare for it. */
    if (CanApplyNormalConeTest = GlblApplyNormalConeTest &&
			         NumOfZeroMVs == LclMVs[0] -> Dim) {
        for (i = 0; i < NumOfMVs; i++) {
	    if (MVAR_IS_RATIONAL_MV(LclMVs[i])) {
	        MVAR_FATAL_ERROR(MVAR_ERR_RATIONAL_NO_SUPPORT);
		return NULL;
	    }
	    else {
	        int Size = sizeof(CagdRType) * MVAR_CTL_MESH_LENGTH(LclMVs[i]);
		MvarMVGradientStruct
		    *Grad = MvarMVPrepGradient(LclMVs[i]);

	        /*   Move the gradient vector to this LclMVs[i], in place.   */
	        /*   The LclMVs[i] was a scalar E1 function and now it will  */
	        /* become E(1+Dim).					     */
	        if (LclMVs[i] -> Dim > CAGD_MAX_PT_COORD) {
		    MVAR_FATAL_ERROR(MVAR_ERR_DIM_TOO_HIGH);   
		    return NULL;
		}

		for (l = 0; l < LclMVs[i] -> Dim; l++) {
		    LclMVs[i] -> Points[l + 2] =
					        (CagdRType *) IritMalloc(Size);
		    GEN_COPY(LclMVs[i] -> Points[l + 2],
			     Grad -> MVGrad -> Points[l + 1], Size);

		    LclMVs[i] -> PType = MVAR_MAKE_PT_TYPE(FALSE,
							 LclMVs[i] -> Dim + 1);
		}

		MvarMVFreeGradient(Grad);
	    }
	}
    }
 
    /* Do the subdivision stage up to SubdivTol tolerance. */
    ZeroSet = MvarMVsSubdivZeros(LclMVs, Constraints, NumOfMVs, NumOfZeroMVs,
				 CanApplyNormalConeTest, SubdivTol, 0);

    if (CanApplyNormalConeTest) {
        for (i = 0; i < NumOfMVs; i++) {
	    if (MVAR_IS_RATIONAL_MV(LclMVs[i])) {
	        MVAR_FATAL_ERROR(MVAR_ERR_RATIONAL_NO_SUPPORT);
		return NULL;
	    }
	    else {
	        int NumOfCoord = MVAR_NUM_OF_PT_COORD(LclMVs[i]);

	        /* Restore scalar multivariates, dropping gradient. */
	        for (l = 2; l <= NumOfCoord; l++) {
		    IritFree(LclMVs[i] -> Points[l]);
		    LclMVs[i] -> Points[l] = NULL;
		    LclMVs[i] -> PType = CAGD_PT_E1_TYPE;
		}
	    }
	}
    }

#   ifdef DEBUG
    {
        IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugMvarSubdivZeroSet, FALSE) {
	    for (Pt = ZeroSet; Pt != NULL; Pt = Pt -> Pnext) {
	        fprintf(stderr, "\n\tSubdivided (%d) = ",
			AttrGetIntAttrib(Pt -> Attr, "SingleSol"));
		for (i = 0; i < Pt -> Dim; i++)
		    fprintf(stderr, "%f  ", Pt -> Pt[i]);
	    }
	}
    }
#   endif /* DEBUG */

    /* Do the numeric improvement stage up to NumericTol tolerance. */
    if (NumericTol < SubdivTol) {
	ZeroSet = MvarMVsNumericZeros(&ZeroSet, LclMVs, Constraints,
				      NumOfZeroMVs, SubdivTol, NumericTol);
    }

    for (i = 0; i < NumOfMVs; i++) {
        if (MVAR_IS_BEZIER_MV(LclMVs[i]))
	    IritFree(LclMVs[i] -> PAux);
	MvarMVFree(LclMVs[i]);
    }
    IritFree(LclMVs);

#   ifdef DEBUG
    {
        IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugMvarFilteringZeroSet, FALSE) {
	    for (Pt = ZeroSet; Pt != NULL; Pt = Pt -> Pnext) {
	        fprintf(stderr, "\n\tUnfiltered = ");
		for (i = 0; i < Pt -> Dim; i++)
		    fprintf(stderr, "%f  ", Pt -> Pt[i]);
	    }
	}
    }
#   endif /* DEBUG */

    /* Filter out points that fail on the inequality constraints or         */
    /* identical points in the zero set.				    */
    if (NumericTol < SubdivTol) {
	while (ZeroSet != NULL) {
	    Pt = ZeroSet;
	    ZeroSet = ZeroSet -> Pnext;
	    Pt -> Pnext = NULL;

	    if (AttrGetIntAttrib(Pt -> Attr, "Similar") == TRUE) {
		MvarPtFree(Pt);
	    }
	    else {
		MvarPtStruct *Pt2;
		CagdBType
		    PurgePt = FALSE;

		/* Lets see if we fail any inequality constraint. */
		for (i = NumOfZeroMVs; i < NumOfMVs; i++) {
		    CagdRType
		        *R = MvarMVEval(MVs[i], Pt -> Pt);
  
		    if (MVAR_IS_RATIONAL_MV(MVs[i]))
		        R[1] = R[1] / R[0];

		    switch(Constraints[i]) {
			case MVAR_CNSTRNT_POSITIVE:
		            if (R[1] < 0.0)
			        PurgePt = TRUE;
		            break;
			case MVAR_CNSTRNT_NEGATIVE:
		            if (R[1] > 0.0)
			        PurgePt = TRUE;
		            break;
		    }
		}

		if (PurgePt) {
		    MvarPtFree(Pt);
		}
		else {
		    for (Pt2 = ZeroSet; Pt2 != NULL; Pt2 = Pt2 -> Pnext) {
		        for (i = 0; i < Pt -> Dim; i++) {
			    if (!APX_EQ_EPS(Pt -> Pt[i], Pt2 -> Pt[i],
					    FABS(NumericTol)))
			        break;
			}

			if (i >= Pt -> Dim) {
			    /* Pt and Pt2 are same - mark Pt2 as similar. */
			    AttrSetIntAttrib(&Pt2 -> Attr, "Similar", TRUE);
			}
		    }

		    LIST_PUSH(Pt, OutSet);
		}
	    }
	}
    }
    else
	OutSet = ZeroSet;

#   ifdef DEBUG
    {
        IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugMvarFilteringZeroSet, FALSE) {
	    for (Pt = OutSet; Pt != NULL; Pt = Pt -> Pnext) {
	        fprintf(stderr, "\n\tFiltered = ");
		for (i = 0; i < Pt -> Dim; i++)
		    fprintf(stderr, "%f  ", Pt -> Pt[i]);
	    }
	}
    }
#   endif /* DEBUG */

    return OutSet;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Copies (and update) the domains' vector of a Bezier multivariate.        *
*                                                                            *
* PARAMETERS:                                                                *
*   MV:         A Bezier multivariate to copy domains from.                  *
*   t:          The parameter value where MV was subdivided at.              *
*   Dir:        The direction where MV was subdivided at.                    *
*   FirstSubMV: TRUE if this copy should reflect the first half of the       *
*		subdivided domain, FALSE to reflect the second.		     *
*                                                                            *
* RETURN VALUE:                                                              *
*   MvarMVParamDomainStruct *:                                               *
*****************************************************************************/
static MvarMVParamDomainStruct *MVarCopySubdivParamDomains(MvarMVStruct *MV,
							   CagdRType t,
							   int Dir,
							   CagdBType FirstSubMV)
{
    MvarMVParamDomainStruct
	*RetDmns = IritMalloc(sizeof(MvarMVParamDomainStruct) * MV -> Dim);

    if (MV -> PAux == NULL)
        MVAR_FATAL_ERROR(MVAR_ERR_INCONS_DOMAIN);

    CAGD_GEN_COPY(RetDmns, MV -> PAux,
		  sizeof(MvarMVParamDomainStruct) * MV -> Dim);

    if (FirstSubMV)
	RetDmns[Dir].Max = t;
    else
	RetDmns[Dir].Min = t;

    return RetDmns;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Extracts the subdivision parametric domain from a multivariate.          *
*                                                                            *
* PARAMETERS:                                                                *
*   MV:         Multivariate to extract doamin from.			     *
*   Min, Max:   Subdivision parametric domain to extract.		     *
*   Dir:        Direction of sought domain.				     *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void MvarGetSubdivParamDomains(MvarMVStruct *MV,
				      CagdRType *Min,
				      CagdRType *Max,
				      int Dir)
{
    if (MVAR_IS_BEZIER_MV(MV)) {
        MvarMVParamDomainStruct
	    *Dmns = (MvarMVParamDomainStruct *) MV -> PAux;

	if (Dmns == NULL)
	    MVAR_FATAL_ERROR(MVAR_ERR_INCONS_DOMAIN);

	*Min = Dmns[Dir].Min;
	*Max = Dmns[Dir].Max;
    }
    else
	MvarMVDomain(MV, Min, Max, Dir);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Construct a point of the dimension as the given MV in the middle of its  *
* parametric domain.                                                         *
*                                                                            *
* PARAMETERS:                                                                *
*   MV:         To consurct a point in the middle of its domain.             *
*   SingleSol:  If TRUE, this point is a single solution it MV domain.       *
*                                                                            *
* RETURN VALUE:                                                              *
*   MvarPtStruct *:  The construct point in the middle of MV.                *
*****************************************************************************/
static MvarPtStruct *MvarGenPtMidMvar(MvarMVStruct *MV, int SingleSol)
{
    int l;
    MvarPtStruct
        *Pt = MvarPtNew(MV -> Dim);

    for (l = 0; l < MV -> Dim; l++) {
	CagdRType Min, Max;

	MvarGetSubdivParamDomains(MV, &Min, &Max, l);
	Pt -> Pt[l] = (Min + Max) * 0.5;
    }

    AttrSetIntAttrib(&Pt -> Attr, "SingleSol", SingleSol);

    return Pt;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Approximate a solution to the set of constraints, if any, using the      *
* subdivision of the parametric domains of the MVs.  Stops when the          *
* parametric domain is smaller than SubdivTol in all dimensions and returns  *
* a central point to that small multivariate patch.			     *
*                                                                            *
* PARAMETERS:                                                                *
*   MVs:          Vector of multivariate constraints.                        *
*   Constraints:  Either an equality or an inequality type of constraint.    *
*   NumOfMVs:     Size of the MVs and Constraints vector.                    *
*   NumOfZeroMVs: Number of zero or equality constraints.                    *
*   ApplyNormalConeTest:  TRUE to apply normal cones' single intersection    *
*		  tests.						     *
*   SubdivTol:    Tolerance of the subdivision process.  Tolerance is        *
*		  measured in the parametric space of the multivariates.     *
*   Depth:        Of subdivision recursion.				     *
*                                                                            *
* RETURN VALUE:                                                              *
*   MvarPtStruct *:   List of points on the solution set.  Dimension of the  *
*		      points will be the same as the dimensions of all MVs.  *
*****************************************************************************/
static MvarPtStruct *MvarMVsSubdivZeros(MvarMVStruct **MVs,
					MvarConstraintType *Constraints,
					int NumOfMVs,
					int NumOfZeroMVs,
					int ApplyNormalConeTest,
					CagdRType SubdivTol,
					int Depth)
{
    int i, l;
    CagdRType Min, Max;

    /* Test if the multivariate may satisfy their constraints.  Compute a    */
    /* bound on the set of values that the multivariate can assume and check */
    /* against the prescribed constraint to that multivariate.		     */
    for (i = 0; i < NumOfMVs; i++) {
	CagdBBoxStruct BBox;

	MvarMVBBox(MVs[i], &BBox);

	switch (Constraints[i]) {
	    case MVAR_CNSTRNT_ZERO:
	        if (BBox.Min[0] > 0.0 || BBox.Max[0] < 0.0)
		    return NULL;
		break;
	    case MVAR_CNSTRNT_POSITIVE:
		if (BBox.Max[0] < 0.0)
		    return NULL;
		break;
	    case MVAR_CNSTRNT_NEGATIVE:
		if (BBox.Min[0] > 0.0)
		    return NULL;
		break;
	}
    }

    /* Check the normal cone overlapping criteria, if at enough depth and    */
    /* every second call only.						     */
    if (ApplyNormalConeTest && !MvarMVConesOverlap(MVs, NumOfZeroMVs)) {
#	ifdef DEBUG
        {
	    IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugPrintZeroPts, FALSE) {
	        CagdRType Min, Max;

		fprintf(stderr, "Zero by cones at");
		for (i = 0; i < MVs[0] -> Dim; i++) {
		    MvarGetSubdivParamDomains(MVs[0], &Min, &Max, i);
		    fprintf(stderr, " %f [%f %f]", (Min + Max) * 0.5, Min, Max);
		}
		fprintf(stderr, "\n");
	    }
	}
#       endif /* DEBUG */

        return MvarGenPtMidMvar(MVs[0], TRUE);
    }

    /* If we got here then these patches may satisfy all constraints.        */
    /* Subdivide them along their maximal parametric length dimension.	     */
    MvarGetSubdivParamDomains(MVs[0], &Min, &Max, l = 0);
    for (i = 1; i < MVs[0] -> Dim; i++) {
	CagdRType MinCurr, MaxCurr;

	MvarGetSubdivParamDomains(MVs[0], &MinCurr, &MaxCurr, i);
        if (MaxCurr - MinCurr > Max - Min) {
            l = i;
            Min = MinCurr;
            Max = MaxCurr;
        }
    }

    if (Max - Min > SubdivTol) {
	MvarPtStruct *PtList1, *PtList2;
	MvarMVStruct
	    **MVs1 = (MvarMVStruct **) IritMalloc(NumOfMVs *
						    sizeof(MvarMVStruct *)),
	    **MVs2 = (MvarMVStruct **) IritMalloc(NumOfMVs *
						    sizeof(MvarMVStruct *));
	CagdRType
	    t = (Min + Max) * 0.5;

	/* Lets see if we have a Bspline multivariate with interior knot */
	/* in this direction.  If so pick up the interior knot instead.  */
	for (i = 0; i < NumOfMVs; i++) {
	    if (MVAR_IS_BSPLINE_MV(MVs[i]) &&
		MVs[i] -> Lengths[l] != MVs[i] -> Orders[l]) {
		t = MVs[i] -> KnotVectors[l][(MVs[i] -> Lengths[l] +
						MVs[i] -> Orders[l]) >> 1];
		break;
	    }
	}

	for (i = 0; i < NumOfMVs; i++) {
	    MVs1[i] = MvarMVSubdivAtParam(MVs[i],
				    MVAR_IS_BEZIER_MV(MVs[i]) ? 0.5 : t, l);
	    MVs2[i] = MVs1[i] -> Pnext;
	    MVs1[i] -> Pnext = NULL;

	    if (MVAR_IS_BEZIER_MV(MVs[i])) {
		/* Needs to copy and update the domains kept in PAux. */
		MVs1[i] -> PAux = MVarCopySubdivParamDomains(MVs[i], t, l,
								TRUE);
		MVs2[i] -> PAux = MVarCopySubdivParamDomains(MVs[i], t, l,
								FALSE);
	    }
	}
	PtList1 = MvarMVsSubdivZeros(MVs1, Constraints, NumOfMVs,
				     NumOfZeroMVs, ApplyNormalConeTest,
				     SubdivTol, Depth + 1);
	PtList2 = MvarMVsSubdivZeros(MVs2, Constraints, NumOfMVs,
				     NumOfZeroMVs, ApplyNormalConeTest,
				     SubdivTol, Depth + 1);

	for (i = 0; i < NumOfMVs; i++) {
	    if (MVAR_IS_BEZIER_MV(MVs[i])) {
		IritFree(MVs1[i] -> PAux);
		IritFree(MVs2[i] -> PAux);
	    }

	    MvarMVFree(MVs1[i]);
	    MvarMVFree(MVs2[i]);
	}
	IritFree(MVs1);
	IritFree(MVs2);

	return (MvarPtStruct *) CagdListAppend(PtList1, PtList2);
    }

#   ifdef DEBUG
    {
        IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugPrintZeroPts, FALSE) {
	    CagdRType Min, Max;

	    fprintf(stderr, "Zero by subdivision level at");
	    for (i = 0; i < MVs[0] -> Dim; i++) {
	        MvarGetSubdivParamDomains(MVs[0], &Min, &Max, i);
		fprintf(stderr, " %f [%f %f]", (Min + Max) * 0.5, Min, Max);
	    }
	    fprintf(stderr, "\n");
	}
    }
#   endif /* DEBUG */

    /* If we got here then in all parametric directions, the domain is       */
    /* smaller than the SubdivTol.  Return central location of this patch.   */
    return MvarGenPtMidMvar(MVs[0], FALSE);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Apply a numerical improvement stage, as a first order minimization       *
* procedure of gradient computation and marching.                            *
*                                                                            *
* PARAMETERS:                                                                *
*   ZeroSet:	  Set of approximated solution, derived from a subdivision   *
*		  process.						     *
*   MVs:          Vector of multivariate constraints.                        *
*   Constraints:  Either an equality or an inequality type of constraint.    *
*   NumOfMVs:     Size of the MVs and Constraints vector.                    *
*   Neighborhood: Neighborhood of the exact solution.		             *
*   NumericTol:   Tolerance of the numerical process.  Tolerance is measured *
*		  in the deviation of the scalar multivariates from their    *
*		  equality. Inequalities are ignored here.  If NumericTol is *
*		  negative, points that fail to improve numerically are      *
*		  purged away.						     *
*                                                                            *
* RETURN VALUE:                                                              *
*   MvarPtStruct *:   List of points on the solution set.  Dimension of the  *
*		      points will be the same as the dimensions of all MVs.  *
*		      Points that failed to improve are purged away.	     *
*****************************************************************************/
static MvarPtStruct *MvarMVsNumericZeros(MvarPtStruct **ZeroSet,
					 MvarMVStruct **MVs,
					 MvarConstraintType *Constraints,
					 int NumOfMVs,
					 CagdRType Neighborhood,
					 CagdRType NumericTol)
{
    MvarPtStruct
	*ImprovedSet = NULL;
    CagdBType
	PurgeFailure = NumericTol < 0.0;
    int i,
	Dim = MVs[0] -> Dim;
    CagdRType *A, *x, *b;

    NumericTol = FABS(NumericTol);

    /* Make sure this is not an over-constrained system of equations. */
    if (Dim < NumOfMVs || NumOfMVs == 0)
	return *ZeroSet;

    /* Prepare gradient functions for faster evaluation. */
    for (i = 0; i < NumOfMVs; i++)
        AttrSetPtrAttrib(&MVs[i] -> Attr, "Gradient",
			 MvarMVPrepGradient(MVs[i]));

    A = (CagdRType *) IritMalloc(sizeof(CagdRType) * Dim * NumOfMVs);
    x = (CagdRType *) IritMalloc(sizeof(CagdRType) * Dim);
    b = (CagdRType *) IritMalloc(sizeof(CagdRType) * Dim);

    while (*ZeroSet != NULL) {
	int Count = 0;
	MvarPtStruct
	    *Pt = *ZeroSet;
	CagdRType NewErr,
	    StartErr = MvarEvalErrorL1(MVs, Pt -> Pt, NumOfMVs),
	    Err = StartErr;

#	ifdef DEBUG
	{
	    IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugMvarZeroErr, FALSE) {
	        fprintf(stderr, "Start Error = %17.14f ", StartErr);
		for (i = 0; i < MVs[0] -> Dim; i++)
		    fprintf(stderr, " %f", Pt -> Pt[i]);
		fprintf(stderr, "\n");
	    }
	}
#	endif /* DEBUG */

	*ZeroSet = (*ZeroSet) -> Pnext;
	Pt -> Pnext = NULL;

	do {
	    CagdRType r, MinR, Min, Max;

	    /* We now seek a point that is on each tangent hyperplanes.      */
	    /* Hence, we have Dim linear equations and Dim dofs.             */
	    for (i = 0; i < NumOfMVs; i++) {
		/* Derive tangent hyperplane for multivariate at current pos.*/
		MvarPlaneStruct
		    *Pln = MvarMVEvalTanPlane(MVs[i], Pt -> Pt);

		/* Copy the constraint into the matrix form. */
		CAGD_GEN_COPY(&A[i * Dim], Pln -> Pln,
			      sizeof(CagdRType) * Dim);
		b[i] = -Pln -> Pln[Dim + 1];

		/* Free the computed plane. */
		MvarPlaneFree(Pln);
	    }
	    
	    /* Solve, possibly under constrainted, system of lin. equations. */
	    if (IritQRUnderdetermined(A, NULL, NULL, NumOfMVs, Dim))
		break;

	    /* Add current position to b vector, seeking relative solution. */
	    for (i = 0; i < NumOfMVs; i++) {
		int j;
		CagdRType
		    t = 0.0;

		for (j = 0; j < Dim; j++)
		    t += A[i * Dim + j] * Pt -> Pt[j];

		b[i] -= t;
	    }

	    /* Solve for the new multivariate(s) coefficients. */
	    IritQRUnderdetermined(NULL, x, b, NumOfMVs, Dim);

	    /* And add the current location back, getting absolute result. */
	    for (i = 0; i < Dim; i++)
		x[i] += Pt -> Pt[i];

#	    ifdef DEBUG
	    {
	        IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugMvarZeroErr, FALSE) {
		    fprintf(stderr, "\tOriginal solution ");
		    for (i = 0; i < MVs[0] -> Dim; i++)
		        fprintf(stderr, " %f", x[i]);
		    fprintf(stderr, "\n");
		}
	    }
#	    endif /* DEBUG */

#	    ifdef DEBUG
	    {
	        /* See if x is in the neighborhood and is in the domain. */
	        IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugExamineX, FALSE) {
		    for (i = 0; i < Dim; i++) {
		        MvarMVDomain(MVs[0], &Min, &Max, i);

			if (FABS(x[i] - Pt -> Pt[i]) > Neighborhood) {
			    if (x[i] > Pt -> Pt[i])
			        x[i] = Pt -> Pt[i] + Neighborhood;
			    else
			        x[i] = Pt -> Pt[i] - Neighborhood;
			}
			x[i] = BOUND(x[i], Min, Max);
		    }
		}
	    }
#	    endif /* DEBUG */

	    /* If we bounce our limits, compute exact distance to the        */
	    /* boundary and move that much exactly.			     */
	    MinR = 1.0;
	    for (i = 0; i < Dim; i++) {
		MvarMVDomain(MVs[0], &Min, &Max, i);
		if (x[i] != BOUND(x[i], Min, Max)) {
		    if (x[i] < Min && !APX_EQ(Pt -> Pt[i], Min))
			r = (Min - x[i]) / (Pt -> Pt[i] - x[i]);
		    else if (x[i] > Max && !APX_EQ(Pt -> Pt[i], Max))
			r = (x[i] - Max) / (x[i] - Pt -> Pt[i]);
		    else
		        r = IRIT_INFNTY;

		    if (MinR > r)
			MinR = r;
		}
	    }
	    if (MinR < 1.0) {
		/* Walk in direction MinR distance, getting to boundary. */
		for (i = 0; i < Dim; i++)
		    x[i] = x[i] * (1.0 - MinR) + Pt -> Pt[i] * MinR;
	    }

	    for (i = 0; i < Dim; i++) {
		MvarMVDomain(MVs[0], &Min, &Max, i);
		x[i] = BOUND(x[i], Min, Max);   /* Prevent round off errors. */
	    }

	    /* Now the million dollar question - is x a better solution!? */
	    NewErr = MvarEvalErrorL1(MVs, x, NumOfMVs);

#	    ifdef DEBUG
	    {
	        IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugMvarZeroErr, FALSE) {
		    fprintf(stderr, "\tError = %17.14f ", NewErr);
		    for (i = 0; i < MVs[0] -> Dim; i++)
		        fprintf(stderr, " %f", x[i]);
		    fprintf(stderr, "\n");
		}
	    }
#	    endif /* DEBUG */

	    if (NewErr < Err) {
		Err = NewErr;
		CAGD_GEN_COPY(Pt -> Pt, x, sizeof(CagdRType) * Dim);
	    }
	    else {
		/* Try half that step size, as a last resort. */
		for (i = 0; i < Dim; i++)
		    x[i] = (x[i] + Pt -> Pt[i]) * 0.5;
		NewErr = MvarEvalErrorL1(MVs, x, NumOfMVs);

#	        ifdef DEBUG
		{
		    IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugMvarZeroErr, FALSE) {
		        fprintf(stderr, "\tError Half = %17.14f ", NewErr);
			for (i = 0; i < MVs[0] -> Dim; i++)
			    fprintf(stderr, " %f", x[i]);
			fprintf(stderr, "\n");
		    }
		}
#	        endif /* DEBUG */

		if (NewErr < Err) {
		    Err = NewErr;
		    CAGD_GEN_COPY(Pt -> Pt, x, sizeof(CagdRType) * Dim);
		}
		else
		    break;
	    }
	}
	while (Err > NumericTol && Count++ <  MVAR_NUMER_ZERO_NUM_STEPS);

	if (PurgeFailure && Err > NumericTol) {
	    MvarPtFree(Pt);
	}
	else {
	    AttrSetRealAttrib(&Pt -> Attr, "Error", Err);
	    LIST_PUSH(Pt, ImprovedSet);
	}
    }

    IritFree(A);
    IritFree(x);
    IritFree(b);

    /* Prepare gradient functions for faster evaluation. */
    for (i = 0; i < NumOfMVs; i++) {
	MvarMVFreeGradient(AttrGetPtrAttrib(MVs[i] -> Attr, "Gradient"));
	AttrFreeOneAttribute(&MVs[i] -> Attr, "Gradient");
    }

    *ZeroSet = ImprovedSet;
    return ImprovedSet;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Computes the L1 error of the current position.                           *
*                                                                            *
* PARAMETERS:                                                                *
*   MVs:        The multivariates to evaluate error for.                     *
*   Params:     The location where the error is to be evaluated.             *
*   NumOfMVs:   Number of multivariates we have.                             *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdRType:     Computed error.                                           *
*****************************************************************************/
static CagdRType MvarEvalErrorL1(MvarMVStruct **MVs,
				 CagdRType *Params,
				 int NumOfMVs)
{
    int i;
    CagdRType
	Err = 0.0;

    for (i = 0; i < NumOfMVs; i++) {
	CagdRType
	    *R = MvarMVEval(MVs[i], Params);

	if (MVAR_IS_RATIONAL_MV(MVs[i]))
	    Err += FABS(R[1] / R[0]);
	else
	    Err += FABS(R[1]);
    }

    return Err;
}