equations.c

Go to the documentation of this file.

#include "equations.h"

#include "defines.h"
#include "error.h"
#include "geom.h"

#include <math.h>

void SetEquationInfo(Tequation *e,TequationInfo *ei);

Tequation *GetOriginalEquation(TequationInfo *ei);

void CopyEquationInfo(TequationInfo *ei_dst,TequationInfo *ei_src);

void GetFirstOrderApproximationToEquation(Tbox *b,double *c,
                                          TLinearConstraint *lc,
                                          TequationInfo *ei);
void ErrorDueToVariable(unsigned int m,double *c,Tinterval *is,double *ev,TequationInfo *ei);

boolean  LinearizeSaddleEquation(double epsilon,
                                 Tbox *b,
                                 TSimplex *lp,TequationInfo *ei);


boolean LinearizeParabolaEquation(double epsilon,
                                  Tbox *b,
                                  TSimplex *lp,TequationInfo *ei);

boolean LinearizeCircleEquation(double epsilon,
                                Tbox *b,
                                TSimplex *lp,TequationInfo *ei);

boolean LinearizeSphereEquation(double epsilon,
                                Tbox *b,
                                TSimplex *lp,TequationInfo *ei);

boolean LinearizeGeneralEquationInt(unsigned int cmp,
                                    double epsilon,
                                    Tbox *b,
                                    double *c,
                                    TSimplex *lp,
                                    TequationInfo *ei);
boolean LinearizeGeneralEquation(double epsilon,
                                 Tbox *b,
                                 TSimplex *lp,TequationInfo *ei);

void DeleteEquationInfo(TequationInfo *ei);


/*******************************************************************************/

void SetEquationInfo(Tequation *eq,TequationInfo *ei)
{
  unsigned int i,j,nf;
  Tvariable_set *vars,*varsf;
  Tequation dd;
  Tmonomial *f;
  unsigned int s;
  Tinterval bounds;
  double v;


  if (LinearEquation(eq))
    ei->EqType=LINEAR_EQUATION;
  else
    {
      if (ParabolaEquation(eq))
        ei->EqType=PARABOLA_EQUATION;
      else
        {
          if (SaddleEquation(eq))
            ei->EqType=SADDLE_EQUATION;
          else
            {
              if (CircleEquation(eq))
                ei->EqType=CIRCLE_EQUATION;
              else
                {
                  if (SphereEquation(eq))
                    ei->EqType=SPHERE_EQUATION;
                  else
                    ei->EqType=GENERAL_EQUATION;
                }
            }
        }
    }

  if (GetNumMonomials(eq)==0)
    Error("Only linear, quadratic, and bilinear equations are considered");
  else
    {
      NEW(ei->equation,1,Tequation);
      CopyEquation(ei->equation,eq);

      vars=GetEquationVariables(eq);
      ei->n=VariableSetSize(vars);

      if (ei->EqType==LINEAR_EQUATION)
        {
          /* Define a LinearConstraint from the Linear equation.
             The linear constraint is the structure to be directly added to the simplex */

          NEW(ei->lc,1,TLinearConstraint);
          InitLinearConstraint(ei->lc);
          v=GetEquationValue(eq);
          NewInterval(v,v,&bounds);
          SimplexExpandBounds(GetEquationCmp(eq),&bounds);
          SetLinearConstraintError(&bounds,ei->lc);

          nf=GetNumMonomials(eq);
          for(i=0;i<nf;i++)
            {
              f=GetMonomial(i,eq);
              varsf=GetMonomialVariables(f);
              s=VariableSetSize(varsf);
              if (s!=1) /* The ct terms of the equation are all in the right-hand part */
                Error("Non-linear equation disguised as linear (SetEquationInfo)");

              AddTerm2LinearConstraint(GetVariableN(0,varsf),GetMonomialCt(f),ei->lc);
            }

          ei->Jacobian=NULL;
          ei->Hessian=NULL;
        }
      else
        {
          /*Compute the Jacobain/Hessian to be used in the linear approximation*/
              
          NEW(ei->Jacobian,ei->n,Tequation *);
          for(i=0;i<ei->n;i++)
            {
              NEW(ei->Jacobian[i],1,Tequation);
              DeriveEquation(GetVariableN(i,vars),ei->Jacobian[i],eq);
            }
              
          NEW(ei->Hessian,ei->n,double *);
          for(i=0;i<ei->n;i++)
            NEW(ei->Hessian[i],ei->n,double);
              
          for(i=0;i<ei->n;i++)
            {
              for(j=i;j<ei->n;j++)
                {
                  DeriveEquation(GetVariableN(j,vars),&dd,ei->Jacobian[i]);
                      
                  nf=GetNumMonomials(&dd);
                      
                  if (nf==0)
                    {
                      ei->Hessian[i][j]=ei->Hessian[j][i]=-GetEquationValue(&dd);
                    }
                  else
                    Error("Non-quadratic equation in SetEquationInfo");
                      
                  DeleteEquation(&dd);
                }
            }
              
          ei->lc=NULL;  
        }
    }
}

Tequation *GetOriginalEquation(TequationInfo *ei)
{
  return(ei->equation);
}

void CopyEquationInfo(TequationInfo *ei_dst,TequationInfo *ei_src)
{
  unsigned int i,j;

  NEW(ei_dst->equation,1,Tequation);
  CopyEquation(ei_dst->equation,ei_src->equation);

  ei_dst->n=ei_src->n;

  ei_dst->EqType=ei_src->EqType;

  if (ei_dst->EqType==LINEAR_EQUATION)
    {
      NEW(ei_dst->lc,1,TLinearConstraint);
      CopyLinearConstraint(ei_dst->lc,ei_src->lc);

      ei_dst->Jacobian=NULL;
      ei_dst->Hessian=NULL;
    }
  else
    {
      NEW(ei_dst->Jacobian,ei_dst->n,Tequation *);
      for(i=0;i<ei_dst->n;i++)
        {
          NEW(ei_dst->Jacobian[i],1,Tequation);
          CopyEquation(ei_dst->Jacobian[i],ei_src->Jacobian[i]);
        }
      
      NEW(ei_dst->Hessian,ei_dst->n,double *);
      for(i=0;i<ei_dst->n;i++)
        {
          NEW(ei_dst->Hessian[i],ei_dst->n,double);
          
          for(j=0;j<ei_dst->n;j++)
            ei_dst->Hessian[i][j]=ei_src->Hessian[i][j];
        }
      
      ei_dst->lc=NULL; 
    }
}


void GetFirstOrderApproximationToEquation(Tbox *b,double *c,
                                          TLinearConstraint *lc,
                                          TequationInfo *ei)
{
  double d;
  Tinterval a;
  unsigned int i,j,k;
  Tvariable_set *vars;
  Tinterval *is;
  unsigned int m;
  Tinterval *cInt,error,dfc,o;
  double v;

  m=GetBoxNIntervals(b);
  is=GetBoxIntervals(b);
 
  /* f(x) = f(c) + df(c) (x-c) + E2 */
  /*     with E2 a quadratic error term */
  /* f(x) = df(c) x + f(c) - df(c) c + E2 */
  /*     f(c) and df(c) are evaluated UP and DOWN
         to avoid missing solutions */
  /* f(x) = [df(c)_l df(c)_u] x + [f(c)_l f(c)_u] - [df(c)_l df(c)_u] c + E2 */
  /* f(x) = df(c)_l x + [0 df(c)_u-df(c)_l] x + [f(c)_l f(c)_u] - [df(c)_l df(c)_u] c + E2  */
  /*     We accumulate all the intervals into one error term */
  /*     Error = [0 df(c)_u-df(c)_l] x + [f(c)_l f(c)_u] - [df(c)_l df(c)_u] c + E2  */
  /*     Error = [0 df(c)_u-df(c)_l] x + [f(c)_l f(c)_u] + [-df(c)_u -df(c)_l] c + E2  */
  /* f(x) = df(c)_l x +  Error */
  

  /* To take into account floating point errors, we perform interval evaluation of the
     equations but on punctual intervals (defined from the given points).  */
  
  vars=GetEquationVariables(ei->equation);

  NEW(cInt,m,Tinterval);
  for(i=0;i<ei->n;i++)
    {
      k=GetVariableN(i,vars);
      NewInterval(c[k],c[k],&(cInt[k]));
    }

  EvaluateEquationInt(cInt,&error,ei->equation); 
  v=GetEquationValue(ei->equation);
  
  /* The coefficients and value of the equations can be affected by some noise
     due to floating point roundings when operating them. We add a small
     range (1e-11) to compensate for those possible errors. */
  NewInterval(-v-ZERO,-v+ZERO,&o);
  
  IntervalAdd(&error,&o,&error);      

  InitLinearConstraint(lc);

  for(i=0;i<ei->n;i++)
    {
      k=GetVariableN(i,vars);

      EvaluateEquationInt(cInt,&dfc,ei->Jacobian[i]); 
      v=GetEquationValue(ei->Jacobian[i]);
      /* The coefficients and value of the equations can be affected by some noise
         due to floating point roundings when operating them. We add a small
         range (1e-11) to compensate for those possible errors. */
      NewInterval(-v-ZERO,-v+ZERO,&o);

      IntervalAdd(&dfc,&o,&dfc);

      d=LowerLimit(&dfc);

      IntervalOffset(&dfc,-d,&a);
      IntervalProduct(&a,&(is[k]),&a);
      IntervalAdd(&error,&a,&error);

      IntervalInvert(&dfc,&a);
      IntervalScale(&a,c[k],&a);
      IntervalAdd(&error,&a,&error);

      AddTerm2LinearConstraint(k,d,lc);
    }
  
  /* Now we bound the error using interval arithmetics */
  /*   Error= 0.5 * \sum_{i,j} r[i] H_{i,j} r[j] */

  /* The interval-based bound is exact as far as the Hessian is ct (i.e., we only
     have quadratic functions) and if there are not repeated variables in the
     error expression ( 0.5 * \sum_{i,j} r[i] H_{i,j} r[j]). This is the case
     of the circle, vector product and scalar product that appear in our problems.
     In other cases we could over-estimate the error (this delays the convergence,
     but does not avoid it)
  */

  /* We re-use the cInt vector previously allocated. This is not a problem since cInt
     is not used any more and the size of cInt (num var in the problem) is larger than
     ei->n (number of variables in the current equation) */
  for(i=0;i<ei->n;i++)
    {
      k=GetVariableN(i,vars);
      IntervalOffset(&(is[k]),-c[k],&(cInt[i])); 
    }

  for(i=0;i<ei->n;i++)
    {
      for(j=0;j<ei->n;j++)
        {
          if (ei->Hessian[i][j]!=0.0)
            {
              if (i==j)
                IntervalPow(&(cInt[i]),2,&a);
              else
                IntervalProduct(&(cInt[i]),&(cInt[j]),&a);
              
              /* The coefficients and value of the equations can be affected by some noise
                 due to floating point roundings when operating them. We add a small
                 range (1e-11) to compensate for those possible errors. */
              NewInterval((0.5*ei->Hessian[i][j])-ZERO,(0.5*ei->Hessian[i][j])+ZERO,&o);
              IntervalProduct(&a,&o,&a);

              IntervalAdd(&error,&a,&error);
            }
        }
    }
  
  free(cInt);
  
  /*  If one of the input ranges infinite, we override the computed error
      (probably have overflow....) */
  if (GetBoxMaxSizeVarSet(vars,b)<INF)
    IntervalInvert(&error,&error);
  else
    NewInterval(-INF,INF,&error);

  SetLinearConstraintError(&error,lc);
}

void ErrorDueToVariable(unsigned int m,double *c,Tinterval *is,
                        double *ev,TequationInfo *ei)
{

  if (ei->EqType!=LINEAR_EQUATION)
    {
      Tvariable_set *vars;
      unsigned int i,k;
      Tinterval error;
      Tinterval *r,a;
      unsigned int j;
      /*
        Error= 1/2 * \sum_{i,j} r[i] * H_{i,j} * r[j] -> 
        Error= \sum{i} Error(i) 
        with
        Error(i)= 1/2 * r[i] * \sum_{j}  H_{i,j}* r[j]

        The output of this function is ev[i]=Error(i) (defined as above)
      */

      vars=GetEquationVariables(ei->equation);
  
      NEW(r,ei->n,Tinterval);

      for(i=0;i<ei->n;i++)
        {
          k=GetVariableN(i,vars);
          IntervalOffset(&(is[k]),-c[k],&(r[i])); 
        }

      for(i=0;i<ei->n;i++)
        {
          NewInterval(0,0,&error);
          for(j=0;j<ei->n;j++)
            {
              if (ei->Hessian[i][j]!=0.0)
                {
                  if (i==j)
                    IntervalPow(&(r[i]),2,&a);
                  else
                    IntervalProduct(&(r[i]),&(r[j]),&a);
                  
                  IntervalScale(&a,ei->Hessian[i][j],&a);
                  IntervalAdd(&error,&a,&error);
                }
            }
          
          ev[i]=IntervalSize(&error)/2.0;
        }

      free(r);  
    }

}

void DeleteEquationInfo(TequationInfo *ei)
{
  unsigned int i;

  DeleteEquation(ei->equation);
  free(ei->equation);

  if (ei->EqType==LINEAR_EQUATION)
    {
      DeleteLinearConstraint(ei->lc);
      free(ei->lc);
    }
  else
    {
      for(i=0;i<ei->n;i++)
        {
          DeleteEquation(ei->Jacobian[i]);
          free(ei->Jacobian[i]);
        }
      free(ei->Jacobian);
      
      for(i=0;i<ei->n;i++)
        free(ei->Hessian[i]);
      
      free(ei->Hessian);
    }
}

void InitEquations(Tequations *eqs)
{
  unsigned int i;

  eqs->n=0;
  eqs->s=0;
  eqs->c=0;
  eqs->d=0;
  eqs->e=0;
  eqs->max=INIT_NUM_EQUATIONS;
  NEW(eqs->equation,eqs->max,TequationInfo *);
  
  for(i=0;i<eqs->max;i++)
    eqs->equation[i]=NULL;
}

void CopyEquations(Tequations *eqs_dst,Tequations *eqs_src)
{
  unsigned int i;

  eqs_dst->n=eqs_src->n;
  eqs_dst->s=eqs_src->s;
  eqs_dst->c=eqs_src->c;
  eqs_dst->d=eqs_src->d;
  eqs_dst->e=eqs_src->e;

  eqs_dst->max=eqs_src->max;

  NEW(eqs_dst->equation,eqs_dst->max,TequationInfo *);

  for(i=0;i<eqs_src->max;i++)
    {
      if (eqs_src->equation[i]!=NULL)
        {
          NEW(eqs_dst->equation[i],1,TequationInfo);
          CopyEquationInfo(eqs_dst->equation[i],eqs_src->equation[i]);
        }
      else
        eqs_dst->equation[i]=NULL;
    }
}

boolean UsedVarInNonDummyEquations(unsigned int nv,Tequations *eqs)
{
  Tequation *eq;
  boolean found;
  unsigned int i;
  
  i=0;
  found=FALSE;
  while((i<eqs->n)&&(!found))
    {
      eq=GetEquation(i,eqs);
      if (GetEquationType(eq)!=DUMMY_EQ)
        found=VarIncluded(nv,GetEquationVariables(eq));

      if (!found) i++;
    }

  return(found);
}

boolean UsedVarInEquations(unsigned int nv,Tequations *eqs)
{
  boolean found;
  unsigned int i;
  
  i=0;
  found=FALSE;
  while((i<eqs->n)&&(!found))
    {
      found=VarIncluded(nv,GetEquationVariables(GetEquation(i,eqs)));
      if (!found) i++;
    }

  return(found);
}

void RemoveEquationsWithVar(double epsilon,unsigned int nv,Tequations *eqs)
{
  Tequation *eq;
  Tequations newEqs;
  unsigned int j;

  InitEquations(&newEqs);
  for(j=0;j<eqs->n;j++)
    {
      
      eq=GetEquation(j,eqs);
      if (!VarIncluded(nv,GetEquationVariables(eq)))
        {
          /*The equations with the given variables are not included in the new set*/
          /*The given variables is to be removed from the problem (it is not used
            any more!). Removing a variable causes a shift in all the variable
            indexes above the removed one. This effect can be easily simulated
            by fixing the variable-to-be-removed for the equations to be kept in 
            the new set. */
          FixVariableInEquation(epsilon,nv,1,eq);
          AddEquation(eq,&newEqs);
        }
    }

  DeleteEquations(eqs);
  CopyEquations(eqs,&newEqs);
  DeleteEquations(&newEqs);
}

boolean ReplaceVariableInEquations(double epsilon,
                                   unsigned int nv,TLinearConstraint *lc,
                                   Tequations *eqs)
{
  boolean hold;
  unsigned int j,r;
  Tequations newEqs;
  Tequation *eqn;
  
  hold=TRUE;
  InitEquations(&newEqs);
  for(j=0;((hold)&&(j<eqs->n));j++)
    {
      eqn=GetEquation(j,eqs);

      r=ReplaceVariableInEquation(epsilon,nv,lc,eqn);

      if (r==2)
        hold=FALSE; /*the equation became ct and does not hold -> the whole system fails*/
      else
        {
          if (r==0)
            AddEquation(eqn,&newEqs);

          /*if r==1, a ct equation that holds -> no need to add it to the new set*/
        }
    }

  DeleteEquations(eqs);
  CopyEquations(eqs,&newEqs);
  DeleteEquations(&newEqs);

  return(hold);
}

boolean GaussianElimination(Tequations *eqs)
{
  unsigned int i,j,k,neq,neqn,meq;
  Tequation *eqi,*eqj,eqMin,eqNew;
  Tequations newEqs;
  unsigned int nfi,nfj,r;
  double ct;
  boolean changes;
  Tmonomial *fi,*fj;

  /*We try to get the simplest possible set of equations
    Simple equation = as less number of monomials as possible*/

  changes=FALSE;

  neq=NEquations(eqs);
  InitEquations(&newEqs);
  for(i=0;i<neq;i++)
    {
      eqi=GetEquation(i,eqs);

      if ((GetEquationCmp(eqi)!=EQU)||(GetEquationType(eqi)==DUMMY_EQ))
        AddEquation(eqi,&newEqs);
      else
        {
          nfi=GetNumMonomials(eqi);

          CopyEquation(&eqMin,eqi);
          
          /* We try to simplify eqi linearly operating it with the
             already simplified equations and the equations still to be
             simplified. */
          neqn=NEquations(&newEqs);
          meq=neqn+(neq-i-1);
          for(j=0;j<meq;j++)
            {
              if (j<neqn)
                eqj=GetEquation(j,&newEqs);
              else
                eqj=GetEquation(i+1+(j-neqn),eqs);
                
              if ((GetEquationCmp(eqj)==EQU)&&(GetEquationType(eqj)!=DUMMY_EQ))
                {
                  nfj=GetNumMonomials(eqj);

                  for(k=0;k<nfi;k++)
                    {
                      fi=GetMonomial(k,eqi);
                      r=FindMonomial(fi,eqj);

                      if (r!=NO_UINT) 
                        {
                          /*If the monomial of eqi is also in eqj*/
                          fj=GetMonomial(r,eqj);
                          ct=-GetMonomialCt(fi)/GetMonomialCt(fj);

                          CopyEquation(&eqNew,eqi);
                          AccumulateEquations(eqj,ct,&eqNew);
                          
                          /*If the equation after operation has less monomials than
                            the minimal equation we have up to now -> we have a new
                            minimal */
                          /*
                          if ((GetNumMonomials(&eqNew)<GetNumMonomials(&eqMin))||
                              ((GetNumMonomials(&eqNew)==GetNumMonomials(&eqMin))&&
                              (fabs(GetEquationValue(&eqNew))<fabs(GetEquationValue(&eqMin)))))*/
                          if ((GetNumMonomials(&eqNew)<GetNumMonomials(&eqMin))||
                              ((GetNumMonomials(&eqNew)==GetNumMonomials(&eqMin))&&
                               (GetEquationNumVariables(&eqNew)<GetEquationNumVariables(&eqMin))))
                            {
                              DeleteEquation(&eqMin);
                              CopyEquation(&eqMin,&eqNew);
                              changes=TRUE;
                            }
                          DeleteEquation(&eqNew);
                        }
                    }
                }
            }
            
          AddEquation(&eqMin,&newEqs); 
          DeleteEquation(&eqMin); 
        }
    }

  DeleteEquations(eqs);
  CopyEquations(eqs,&newEqs);
  DeleteEquations(&newEqs);

  return(changes);
}

unsigned int NEquations(Tequations *eqs)
{
  return(eqs->n);
}

unsigned int NSystemEquations(Tequations *eqs)
{
  return(eqs->s);
}

unsigned int NCoordEquations(Tequations *eqs)
{
  return(eqs->c);
}

unsigned int NDummyEquations(Tequations *eqs)
{
  return(eqs->d);
}

unsigned int NEqualityEquations(Tequations *eqs)
{
  return(eqs->e);
}

boolean IsSystemEquation(unsigned int i,Tequations *eqs)
{
  return((i<eqs->n)&&
         (GetEquationType(GetOriginalEquation(eqs->equation[i]))==SYSTEM_EQ));
}

boolean IsCoordEquation(unsigned int i,Tequations *eqs)
{
  return((i<eqs->n)&&
         (GetEquationType(GetOriginalEquation(eqs->equation[i]))==COORD_EQ));
}

boolean IsDummyEquation(unsigned int i,Tequations *eqs)
{
  return((i<eqs->n)&&
         (GetEquationType(GetOriginalEquation(eqs->equation[i]))==DUMMY_EQ));
}

unsigned int GetEquationTypeN(unsigned int i,Tequations *eqs)
{
  if (eqs->equation[i]==NULL)
    return(NOTYPE_EQ);
  else
    return(GetEquationType(GetOriginalEquation(eqs->equation[i])));
}

boolean HasEquation(Tequation *equation,Tequations *eqs)
{
  boolean found;
  unsigned int i;
  
  i=0;
  found=FALSE;
  while((i<eqs->n)&&(!found))
    {
      found=(CmpEquations(equation,GetOriginalEquation(eqs->equation[i]))==0);
      if (!found) i++;
    }
  return(found);
}

unsigned int CropEquation(unsigned int ne,
                          double epsilon,double rho,
                          Tbox *b,
                          Tequations *eqs)
{
  unsigned int status;
  TequationInfo *ei;
  unsigned int cmp;
  Tequation *eq;
  Tvariable_set *vars;
  
  double val,alpha;
  Tvariable_set *varsf;;
  Tinterval x_new,y_new,z_new;
  unsigned int nx,ny,nz; 
  unsigned int m,i,k;
  Tinterval *is,bounds;

  ei=eqs->equation[ne];

  eq=GetOriginalEquation(ei);

  status=NOT_REDUCED_BOX;

  cmp=GetEquationCmp(eq);

  #if (_DEBUG>6)
    printf("\n       Croping Equation: %u\n",ne);
  #endif

  if (cmp==EQU)
    {
      m=GetBoxNIntervals(b);
      is=GetBoxIntervals(b);

      vars=GetEquationVariables(eq);

      val=GetEquationValue(eq);

      #if (_DEBUG>6)
        printf("         Eq: ");
        PrintEquation(stdout,NULL,eq);
        printf("         Input ranges: ");
        for(i=0;i<ei->n;i++)
          {
            k=GetVariableN(i,vars);
            printf(" v[%u]=",k);PrintInterval(stdout,&(is[k]));printf("  ");
          }
        printf("\n");
      #endif

            
      /* Special equations cropped with special functions that get tighter bounds
         (they use non-linear functions) */
      switch(ei->EqType)
        {
        case LINEAR_EQUATION:
          status=CropLinearConstraint(epsilon,b,ei->lc);
          break;

        case PARABOLA_EQUATION:
          /* x^2 - \alpha y = 0 */
          nx=GetVariableN(0,GetMonomialVariables(GetMonomial(0,eq)));
          ny=GetVariableN(0,GetMonomialVariables(GetMonomial(1,eq)));
          alpha=GetMonomialCt(GetMonomial(1,eq));

          if (((IntervalSize(&(is[nx]))>=epsilon)||
               (IntervalSize(&(is[ny]))>=epsilon))&&
              (IntervalSize(&(is[nx]))<INF))
            {
              if (RectangleParabolaClipping(&(is[nx]),-alpha,&(is[ny]),&x_new,&y_new))
                { 
                  status=REDUCED_BOX;
                  CopyInterval(&(is[nx]),&x_new);
                  CopyInterval(&(is[ny]),&y_new);
                }
              else
                status=EMPTY_BOX;
            }
          break;

        case SADDLE_EQUATION:
          /* x*y - \alpha z = 0 */
          varsf=GetMonomialVariables(GetMonomial(0,eq));
          alpha=GetMonomialCt(GetMonomial(1,eq));

          nx=GetVariableN(0,varsf);
          ny=GetVariableN(1,varsf);
          nz=GetVariableN(0,GetMonomialVariables(GetMonomial(1,eq)));
      
          if (((IntervalSize(&(is[nx]))>=epsilon)||
               (IntervalSize(&(is[ny]))>=epsilon)||
               (IntervalSize(&(is[nz]))>=epsilon))&&
              (IntervalSize(&(is[nx]))<INF)&&
              (IntervalSize(&(is[ny]))<INF))
            {
              if (BoxSaddleClipping(&(is[nx]),&(is[ny]),-alpha,&(is[nz]),
                                    &x_new,&y_new,&z_new))
                { 
                  status=REDUCED_BOX;
                  CopyInterval(&(is[nx]),&x_new);
                  CopyInterval(&(is[ny]),&y_new);
                  CopyInterval(&(is[nz]),&z_new);
                }
              else
                status=EMPTY_BOX;
            }
          break;

        case CIRCLE_EQUATION:
          /* x^2 + y^2 = val */
          nx=GetVariableN(0,vars);
          ny=GetVariableN(1,vars);

          if ((IntervalSize(&(is[nx]))>=epsilon)||
              (IntervalSize(&(is[ny]))>=epsilon))
            {
              if (RectangleCircleClipping(val,
                                          &(is[nx]),&(is[ny]),
                                          &x_new,&y_new))
                { 
                  status=REDUCED_BOX;
                  CopyInterval(&(is[nx]),&x_new);
                  CopyInterval(&(is[ny]),&y_new);
                }
              else
                status=EMPTY_BOX;
            }
          break;

        case SPHERE_EQUATION:
          /* x^2 + y^2 + z^2 = val */
          nx=GetVariableN(0,vars);
          ny=GetVariableN(1,vars);
          nz=GetVariableN(2,vars);
      
          if ((IntervalSize(&(is[nx]))>=epsilon)||
              (IntervalSize(&(is[ny]))>=epsilon)||
              (IntervalSize(&(is[nz]))>=epsilon))
            {
              if (BoxSphereClipping(val,
                                    &(is[nx]),&(is[ny]),&(is[nz]),
                                    &x_new,&y_new,&z_new))
                { 
                  status=REDUCED_BOX;
                  CopyInterval(&(is[nx]),&x_new);
                  CopyInterval(&(is[ny]),&y_new);
                  CopyInterval(&(is[nz]),&z_new);
                }
              else
                status=EMPTY_BOX;
            }
          break;

        case GENERAL_EQUATION:
          /* When simplifying systems it is quite common to introduce equations
             of type x^2=ct or x*y=ct. We treat them as special cases and taking into
             accound the non-linealities -> stronger crop */
          if ((GetNumMonomials(eq)==1)&&
              (GetMonomialCt(GetMonomial(0,eq))==1))
            {
              Tinterval r;

              if (VariableSetSize(vars)>2)
                Error("Three variables in a single monomial?");
                

              if (VariableSetSize(vars)==1) /*only one variable in the equation*/
                {
                  /* Only one variable included in the equation-> x^2=ct
                     (the linal case x=ct is treated above */
                  /* ax^2=ct is transformed to x^2+y^2=ct with y=[0,0] */
                  nx=GetVariableN(0,vars);

                  NewInterval(-ZERO,+ZERO,&r);

                  if (IntervalSize(&(is[nx]))>=epsilon)
                    {
                      if (RectangleCircleClipping(val,
                                                  &(is[nx]),&r,
                                                  &x_new,&y_new))
                        { 
                          status=REDUCED_BOX;
                          CopyInterval(&(is[nx]),&x_new);
                        }
                      else
                        status=EMPTY_BOX;
                    }
                }
              else
                {
                  /* xy=ct is transformed to xy-z=0 with z=ct*/
                  nx=GetVariableN(0,vars);
                  ny=GetVariableN(1,vars);
                  NewInterval(val-ZERO,val+ZERO,&r);

                  if ((IntervalSize(&(is[nx]))>=epsilon)||
                      (IntervalSize(&(is[ny]))>=epsilon))
                    {
                      if (BoxSaddleClipping(&(is[nx]),&(is[ny]),1,&r,
                                            &x_new,&y_new,&z_new))
                        { 
                          status=REDUCED_BOX;
                          CopyInterval(&(is[nx]),&x_new);
                          CopyInterval(&(is[ny]),&y_new);
                        }
                      else
                        status=EMPTY_BOX;
                    }
                }
            }
          else
            {
              double *c; /*central point*/
              TLinearConstraint lc;

              NEW(c,m,double); 

              do {
                /*possible improvement -> use a random point inside the ranges
                  this will cut increase the crop */
                for(i=0;i<ei->n;i++)
                  {
                    k=GetVariableN(i,vars);
                    c[k]=IntervalCenter(&(is[k]));
                  }

                GetFirstOrderApproximationToEquation(b,c,&lc,ei);

                GetLinearConstraintError(&bounds,&lc);
                SimplexExpandBounds(cmp,&bounds);
                SetLinearConstraintError(&bounds,&lc);          

                status=CropLinearConstraint(epsilon,b,&lc);
                DeleteLinearConstraint(&lc);
             
              } while(status==REDUCED_BOX);

              free(c);
            }
          break;
          
        default:
          Error("Unknown equation type in CropEquation");
        }
    }
  #if (_DEBUG>6)
  else
    printf("           Inequalities are not cropped\n");
    
  printf("         Output ranges: ");
  for(i=0;i<ei->n;i++)
    {
      k=GetVariableN(i,vars);
      printf(" v[%u]=",k);PrintInterval(stdout,&(is[k]));printf("  ");
    }
  printf("\n         ");
  #endif     
 
  return(status);
}

void AddEquation(Tequation *equation,Tequations *eqs)
{
  if ((GetEquationType(equation)==NOTYPE_EQ)||(GetEquationCmp(equation)==NOCMP))
    Error("Adding an equation with unknown type of cmp");

  NormalizeEquation(equation);

  if ((GetNumMonomials(equation)>0)&&(!HasEquation(equation,eqs)))
    {
      unsigned int t,j;
      TequationInfo *s;

      if (eqs->n==eqs->max)
        {
          unsigned int i,k;

          k=eqs->max;
          MEM_DUP(eqs->equation,eqs->max,TequationInfo *);

          for(i=k;i<eqs->max;i++)
            eqs->equation[i]=NULL;
        }

      NEW(eqs->equation[eqs->n],1,TequationInfo);
      SetEquationInfo(equation,eqs->equation[eqs->n]);

      t=GetEquationType(equation);
      if (t==SYSTEM_EQ) eqs->s++;
      if (t==COORD_EQ) eqs->c++;
      if (t==DUMMY_EQ) eqs->d++;
      if (GetEquationCmp(equation)==EQU) eqs->e++;

      eqs->n++;

      j=eqs->n-1;
      while((j>0)&&(CmpEquations(GetOriginalEquation(eqs->equation[j-1]),
                                 GetOriginalEquation(eqs->equation[j]))==1))
        {
          s=eqs->equation[j-1];
          eqs->equation[j-1]=eqs->equation[j];
          eqs->equation[j]=s;

          j--;
        }
    }
}

Tequation *GetEquation(unsigned int n,Tequations *eqs)
{
  if (n<eqs->n)
    return(GetOriginalEquation(eqs->equation[n]));
  else
    return(NULL);
}

boolean LinearizeSaddleEquation(double epsilon,Tbox *b,
                                TSimplex *lp,TequationInfo *ei)
{ 
  Tvariable_set *varsf;
  unsigned int nx,ny,nz;
  double x_min,x_max;
  double y_min,y_max;
  double *c;
  boolean full;
  Tequation *eq;
  Tinterval *is;
  unsigned int m;

  m=GetBoxNIntervals(b);
  is=GetBoxIntervals(b);
  eq=GetOriginalEquation(ei);

  /* x*y - z =0  */
  varsf=GetMonomialVariables(GetMonomial(0,eq));
  
  nx=GetVariableN(0,varsf);
  ny=GetVariableN(1,varsf);
  nz=GetVariableN(0,GetMonomialVariables(GetMonomial(1,eq)));

  x_min=LowerLimit(&(is[nx]));
  x_max=UpperLimit(&(is[nx]));

  y_min=LowerLimit(&(is[ny]));
  y_max=UpperLimit(&(is[ny]));

  NEW(c,m,double); 
  
  {
    double px[4]={x_min,x_min,x_max,x_max};
    double py[4]={y_min,y_max,y_min,y_max};
    unsigned int cmp[4]={LEQ,GEQ,GEQ,LEQ};
    
    unsigned int i; 
    
    full=TRUE;
    for(i=0;((full)&&(i<4));i++)
      {
        c[nx]=px[i];
        c[ny]=py[i];
        c[nz]=0;
        
        full=LinearizeGeneralEquationInt(cmp[i],epsilon,b,c,lp,ei);
      }
  }

  free(c);
  
  return(full);
}

boolean LinearizeParabolaEquation(double epsilon,
                                  Tbox *b,
                                  TSimplex *lp,TequationInfo *ei)
{
  unsigned int nx,ny;
  double x_min,x_max,x_c;
  double y_min,y_max;
  double *c;
  boolean full;
  Tequation *eq;
  Tinterval *is;
  unsigned int m;

  m=GetBoxNIntervals(b);
  is=GetBoxIntervals(b);
  eq=GetOriginalEquation(ei);

  /* x^2 - y =0  */
  nx=GetVariableN(0,GetMonomialVariables(GetMonomial(0,eq)));
  ny=GetVariableN(0,GetMonomialVariables(GetMonomial(1,eq)));

  x_min=LowerLimit(&(is[nx]));
  x_max=UpperLimit(&(is[nx]));
  x_c=IntervalCenter(&(is[nx]));

  y_min=LowerLimit(&(is[ny]));
  y_max=UpperLimit(&(is[ny]));

  NEW(c,m,double); 
  
  {
    double px[3]={x_c,x_min,x_max};
    unsigned int cmp[3]={EQU,LEQ,LEQ};
    unsigned int i;

    full=TRUE;
    for(i=0;((full)&&(i<3));i++)
      {
        c[nx]=px[i];
        c[ny]=0;  /*The value of 'y' is not used in the linealization*/
        
        full=LinearizeGeneralEquationInt(cmp[i],epsilon,b,c,lp,ei);
      }
  }

  free(c);

  return(full);
}

boolean LinearizeCircleEquation(double epsilon,
                                Tbox *b,
                                TSimplex *lp,TequationInfo *ei)
{
  Tequation *eq;
  Tvariable_set *vars;
  unsigned int nx,ny;
  double x_min,x_max,x_c;
  double y_min,y_max,y_c;
  double *c;
  double r2; 
  /*At most 5 linearization points are selected (one
    for each box corner and one from the center of the
    interval x,y)*/
  double cx[5];
  double cy[5];
  unsigned int cmp[5];
  unsigned int i,nc; 
  boolean full;
  Tinterval *is;
  unsigned int m;

  m=GetBoxNIntervals(b);
  is=GetBoxIntervals(b);
  eq=GetOriginalEquation(ei);

  vars=GetEquationVariables(eq);

  r2=GetEquationValue(eq);

  nx=GetVariableN(0,vars);
  ny=GetVariableN(1,vars);

  x_min=LowerLimit(&(is[nx]));
  x_max=UpperLimit(&(is[nx]));
  x_c=IntervalCenter(&(is[nx]));

  y_min=LowerLimit(&(is[ny]));
  y_max=UpperLimit(&(is[ny]));
  y_c=IntervalCenter(&(is[ny]));

  NEW(c,m,double); /* space for the linealization point */
  
  /* for large input ranges (both in x and y) we add linealization constraints
     from the corners of the box*/
  nc=0;
  {
    /* Check which of the 4 box corners are out of the circle.
       For the external points we add a linealization point
       selected on the circle: intersection of the circle and a 
       line from the origin to the corner (a normalization of
       point (x,y) to norm sqrt(r2))
    */
    double px[4]={x_min,x_min,x_max,x_max};
    double py[4]={y_min,y_max,y_min,y_max};
    double norm;
    
    for(i=0;i<4;i++)
      {
        ROUNDDOWN;
        norm=px[i]*px[i]+py[i]*py[i];
        ROUNDNEAR;
        if (norm>(r2+epsilon)) /*corner (clearly) out of the circle*/
          {
            norm=sqrt(r2/norm);
            cx[nc]=px[i]*norm;
            cy[nc]=py[i]*norm; 
            cmp[nc]=LEQ;
            nc++;
          }
      }
  }
  
  /* A point selected from the center of x range is added always
     'y' is selected so that it is on the circle, if possible.
     If not, y is just eh center of the y range
  */
  cx[nc]=x_c;
  cy[nc]=y_c;
  //cmp[nc]=EQU;
  cmp[nc]=GetEquationCmp(eq);
  nc++;

  full=TRUE;
  for(i=0;((full)&&(i<nc));i++)
    {
      c[nx]=cx[i];
      c[ny]=cy[i];

      full=LinearizeGeneralEquationInt(cmp[i],epsilon,b,c,lp,ei);
    }

  free(c);

  return(full);
}

boolean LinearizeSphereEquation(double epsilon,
                                Tbox *b,
                                TSimplex *lp,TequationInfo *ei)
{
  Tequation *eq;
  Tvariable_set *vars;
  unsigned int nx,ny,nz;
  double x_min,x_max,x_c;
  double y_min,y_max,y_c;
  double z_min,z_max,z_c;
  double *c;
  double r2; 
  /*At most 9 linearization points are selected (one
    for each box corner and one from the center of the
    interval x,y,z)*/
  double cx[9];
  double cy[9];
  double cz[9];
  unsigned int cmp[9];
  unsigned int i,nc; 
  boolean full;
  Tinterval *is;
  unsigned int m;

  m=GetBoxNIntervals(b);
  is=GetBoxIntervals(b);
  eq=GetOriginalEquation(ei);

  vars=GetEquationVariables(eq);

  r2=GetEquationValue(eq);

  nx=GetVariableN(0,vars);
  ny=GetVariableN(1,vars);
  nz=GetVariableN(2,vars);

  x_min=LowerLimit(&(is[nx]));
  x_max=UpperLimit(&(is[nx]));
  x_c=IntervalCenter(&(is[nx]));

  y_min=LowerLimit(&(is[ny]));
  y_max=UpperLimit(&(is[ny]));
  y_c=IntervalCenter(&(is[ny]));

  z_min=LowerLimit(&(is[nz]));
  z_max=UpperLimit(&(is[nz]));
  z_c=IntervalCenter(&(is[nz]));

  NEW(c,m,double); /* space for the linealization point */
  
  /* for large input ranges (both in x and y) we add linealization constraints
     from the corners of the box*/
  nc=0;

  {
    /* Check which of the 8 box corners are out of the sphere.
       For the external points we add a linealization point
       selected on the sphere: intersection of the sphere and a 
       line from the origin to the corner (a normalization of
       point (x,y,z) to norm sqrt(r2))
    */
    double px[8]={x_min,x_min,x_min,x_min,x_max,x_max,x_max,x_max};
    double py[8]={y_min,y_min,y_max,y_max,y_min,y_min,y_max,y_max};
    double pz[8]={z_min,z_max,z_min,z_max,z_min,z_max,z_min,z_max};
    double norm;
    
    for(i=0;i<8;i++)
      {
        ROUNDDOWN;
        norm=px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i];
        ROUNDNEAR;
        if (norm>(r2+epsilon)) /*corner (clearly) out of the sphere*/
          {
            norm=sqrt(r2/norm);
            cx[nc]=px[i]*norm;
            cy[nc]=py[i]*norm;
            cz[nc]=pz[i]*norm; 
            cmp[nc]=LEQ;
            nc++;
          }
      }
  }
  
  /* A point selected from the center of x range is added always
     'y' is selected so that it is on the circle, if possible.
     If not, y is just eh center of the y range
  */
  cx[nc]=x_c;
  cy[nc]=y_c;
  cz[nc]=z_c; 
  cmp[nc]=EQU;
  nc++;

  full=TRUE;
  for(i=0;((full)&&(i<nc));i++)
    {
      c[nx]=cx[i];
      c[ny]=cy[i];
      c[nz]=cz[i];

      full=LinearizeGeneralEquationInt(cmp[i],epsilon,b,c,lp,ei);
    }

  free(c);

  return(full);
}

boolean LinearizeGeneralEquationInt(unsigned int cmp,
                                    double epsilon,
                                    Tbox *b,double *c,
                                    TSimplex *lp,TequationInfo *ei)
{  
  Tvariable_set *vars;
  TLinearConstraint lc;
  boolean full=TRUE;

  vars=GetEquationVariables(GetOriginalEquation(ei));

  if (GetBoxMaxSizeVarSet(vars,b)<INF)
    {
      GetFirstOrderApproximationToEquation(b,c,&lc,ei);

      full=SimplexAddNewConstraint(epsilon,&lc,cmp,GetBoxIntervals(b),lp);

      DeleteLinearConstraint(&lc);
    }

  return(full);
}

boolean LinearizeGeneralEquation(double epsilon,
                                 Tbox *b,
                                 TSimplex *lp,TequationInfo *ei)
{
  Tequation *eq;
  unsigned int i,k;
  Tvariable_set *vars;
  double *c;
  boolean full;
  Tinterval *is;
  unsigned int m;
  double l,u;

  m=GetBoxNIntervals(b);
  is=GetBoxIntervals(b);

  eq=GetOriginalEquation(ei);

  vars=GetEquationVariables(eq);

  NEW(c,m,double); 
  for(i=0;i<ei->n;i++)
    {
      k=GetVariableN(i,vars);
      l=LowerLimit(&(is[k]));
      u=UpperLimit(&(is[k]));
      if ((l==-INF)||(u==INF))
        c[k]=0.0;
      else
        c[k]=IntervalCenter(&(is[k]));
    }

  full=LinearizeGeneralEquationInt(GetEquationCmp(eq),epsilon,b,c,lp,ei);

  free(c);

  return(full);
}

boolean AddEquation2Simplex(unsigned int ne, /*eq identifier*/
                            double lr2tm_q,double lr2tm_s,
                            double epsilon,
                            double rho,
                            Tbox *b,
                            Tvariables *vs,
                            TSimplex *lp,
                            Tequations *eqs)
{
  boolean full;

  full=TRUE;

  if (ne<eqs->n)/* No equation -> nothing done */
    {
      TequationInfo *ei;
      Tequation *eq;
      Tinterval r;
      double val;
      double s;
      boolean ctEq;

      ei=eqs->equation[ne];
      eq=GetOriginalEquation(ei);

      #if (_DEBUG>4)
        printf("     Adding equation %u of type %u to the simplex \n",ne,ei->EqType);
        printf("       ");
        PrintEquation(stdout,NULL,eq);
      #endif

      if (CropEquation(ne,epsilon,rho,b,eqs)==EMPTY_BOX)
        {
          #if (((!_USE_MPI)&&(_DEBUG==1))||(_DEBUG>4))
            fprintf(stderr,"C");
          #endif
          full=FALSE;
        }
      else
        {
          /*There is no need to add constant equations to the simplex.
            Thus, we evaluate the equations and if the result is almost
            constant we just check if the equation holds or not.
            If not, the system has no solution.
            The same applies to inequalities that already hold
          */

          EvaluateEquationInt(GetBoxIntervals(b),&r,eq);
          val=GetEquationValue(eq);
          IntervalOffset(&r,-val,&r);
      
          #if (_DEBUG>4)
            printf("         EvalIntEq:v");
            PrintInterval(stdout,&r);
            if (IntervalSize(&r)<INF)
              printf(" -> %g\n",IntervalSize(&r));
            else
              printf(" -> +inf\n");
          #endif

          ctEq=FALSE;

          switch(GetEquationCmp(eq))
            {
            case EQU:   
              if (IntervalSize(&r)<=(10*epsilon))
                {
                  /* We have an (almost) constant equation -> just check whether
                     it holds or not but do not add it to the simplex */
                  /* We allow for some numerical inestabilities: numerical values
                     below epsilon are converted to zero when defining the linear
                     constraints and, consequently, we have to admit epsilon-like
                     errors
                  */
                  ctEq=TRUE;
                  full=((IsInside(0.0,&r))||
                        (fabs(LowerLimit(&r))<10*epsilon)||
                        (fabs(UpperLimit(&r))<10*epsilon));

                }
              break;

            case LEQ:
              if (UpperLimit(&r)<0.0)
                {
                  ctEq=TRUE;
                  full=TRUE;
                }
              else
                {
                  if (LowerLimit(&r)>0.0)
                    {
                      ctEq=TRUE;
                      full=FALSE;
                    }
                }
              break;
            case GEQ:
              if (UpperLimit(&r)<0.0)
                {
                  ctEq=TRUE;
                  full=FALSE;
                }
              else
                {
                  if (LowerLimit(&r)>0.0)
                    {
                      ctEq=TRUE;
                      full=TRUE;
                    }
                }
              break;
            }
          
          #if (_DEBUG>4)
            if (ctEq)
              {
                printf("     Constant Equation\n");
                printf("       Eq eval:[%g %g] -> ",LowerLimit(&r),UpperLimit(&r));
                if (full) 
                  printf("ct full equation\n\n");
                else
                  printf("ct empty equation\n\n");
              }
          #endif

          if (!ctEq)
            {
              /* The equation is not constant, we proceed to add it into
                 the simplex tableau, linearizing if needed.*/

              s=GetBoxMinSizeVarSet(GetEquationVariables(eq),b);

              /* Use a particular linear relaxation depending on the
                 equation type */
              switch (ei->EqType)
                {
                case LINEAR_EQUATION:
                  full=SimplexAddNewConstraint(epsilon,
                                               ei->lc,
                                               GetEquationCmp(eq),
                                               GetBoxIntervals(b),lp);
                  break;
                  
                case SADDLE_EQUATION:
                  if (s<lr2tm_s)
                    full=LinearizeGeneralEquation(epsilon,b,lp,eqs->equation[ne]);
                  else
                    full=LinearizeSaddleEquation(epsilon,b,lp,eqs->equation[ne]);
                  break;
                  
                case PARABOLA_EQUATION:
                  if (s<lr2tm_q)
                    full=LinearizeGeneralEquation(epsilon,b,lp,eqs->equation[ne]);
                  else
                    full=LinearizeParabolaEquation(epsilon,b,lp,eqs->equation[ne]);
                  break;
                  
                case CIRCLE_EQUATION:
                  if (s<lr2tm_q)
                    full=LinearizeGeneralEquation(epsilon,b,lp,eqs->equation[ne]);
                  else
                    full=LinearizeCircleEquation(epsilon,b,lp,eqs->equation[ne]);
                  break;
                      
                case SPHERE_EQUATION:
                  if (s<lr2tm_q)
                    full=LinearizeGeneralEquation(epsilon,b,lp,eqs->equation[ne]);
                  else
                    full=LinearizeSphereEquation(epsilon,b,lp,eqs->equation[ne]);
                  break;
                  
                case GENERAL_EQUATION:
                  full=LinearizeGeneralEquation(epsilon,b,lp,eqs->equation[ne]);
                  break;
                  
                default:
                  Error("Unknown equation type in AddEquation2Simplex");
                }

            } 
          #if ((!_USE_MPI)&&(_DEBUG==1))
            if (!full)
              fprintf(stderr,"A");
          #endif     
        }
    }

  return(full);
}

void UpdateSplitWeight(unsigned int ne,double *splitDim,
                       Tbox *b,Tequations *eqs)
{
  unsigned int i,k;
  TequationInfo *ei;
  Tvariable_set *vars;
  double *c,*ev;
  Tinterval *is;
  unsigned int m;

  if (ne<eqs->n)
    {
      m=GetBoxNIntervals(b);
      is=GetBoxIntervals(b);

      ei=eqs->equation[ne];
 
      if (ei->EqType!=LINEAR_EQUATION)
        {
          vars=GetEquationVariables(ei->equation);

          NEW(c,m,double);
          NEW(ev,m,double); 

          for(i=0;i<ei->n;i++)
            {
              k=GetVariableN(i,vars);
              c[k]=IntervalCenter(&(is[k]));
            }

          ErrorDueToVariable(m,c,is,ev,ei);

          for(i=0;i<ei->n;i++)
            {
              k=GetVariableN(i,vars);
              splitDim[k]+=ev[i];
            }

          free(c);
          free(ev);
        }
    }
}

void PrintEquations(FILE *f,char **varNames,Tequations *eqs)
{
  unsigned int i;
  
  if (eqs->s>0)
    {
      fprintf(f,"\n[SYSTEM EQS]\n\n");
      for(i=0;i<eqs->n;i++)
        {
          if (IsSystemEquation(i,eqs))
            PrintEquation(f,varNames,GetOriginalEquation(eqs->equation[i]));
        }
    }

  if (eqs->c>0)
    {
      fprintf(f,"\n[COORD EQS]\n\n");
      for(i=0;i<eqs->n;i++)
        {
          if (IsCoordEquation(i,eqs))
            PrintEquation(f,varNames,GetOriginalEquation(eqs->equation[i]));
        }
    }

  if (eqs->d>0)
    {
      fprintf(f,"\n[DUMMY EQS]\n\n");
      for(i=0;i<eqs->n;i++)
        {
          if (IsDummyEquation(i,eqs))
            PrintEquation(f,varNames,GetOriginalEquation(eqs->equation[i]));
        }
    }
}

void DeleteEquations(Tequations *eqs)
{
  unsigned int i;

  for(i=0;i<eqs->n;i++)
    {    
      if (eqs->equation[i]!=NULL)
        {
          DeleteEquationInfo(eqs->equation[i]);
          free(eqs->equation[i]);
          eqs->equation[i]=NULL;
        }
    }
  eqs->n=0;
  eqs->s=0; 
  eqs->c=0; 
  eqs->d=0; 
  eqs->e=0; 

  free(eqs->equation);
  eqs->equation=NULL;
}