interval.c

Go to the documentation of this file.

#include "interval.h"

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

#include <math.h>

double NormalizeAngle(double a);

/*
 * moves the angle 'a' but in the interval [0,2pi]
 */
double NormalizeAngle(double a)
{
  double new_a;

  new_a=a;
  if ((new_a>M_PI)||(new_a<-M_PI))
    new_a=atan2(sin(new_a),cos(new_a));

  if (new_a<0.0) new_a+=M_2PI;

  return(new_a);
}

/*
 * Defines a new interval with limits: lower and upper
 */
void NewInterval(double lower,      /*lower limit of the new interval*/
                 double upper,      /*upper limit of the new interval*/
                 Tinterval *i       /*new interval*/
                 )
{
  i->lower_limit=lower;
  i->upper_limit=upper;
}

/*
 * Copies i_org into i_dst.
 */
void CopyInterval(Tinterval *i_dst,Tinterval *i_org)
{
  i_dst->lower_limit=i_org->lower_limit;
  i_dst->upper_limit=i_org->upper_limit;
}



/*
 * Returns true if i1 and i2 are exactly the same.
 * This only works for intervals that are exactly the same (copies, the same
 * initialization,....)
 */
boolean CmpIntervals(Tinterval *i1,Tinterval *i2)
{
  return((i1->lower_limit==i2->lower_limit)&&
         (i1->upper_limit==i2->upper_limit));
}

/*
 * Returns the lower limit of interval 'i'.
 */
double LowerLimit(Tinterval *i)
{
  return(i->lower_limit);
}

/*
 * Returns the upper limit of interval 'i'.
 */
double UpperLimit(Tinterval *i)
{
  return(i->upper_limit);
}

/*
 * Returns the size of interval 'i'.
 * NOTE: For empty intervals it can return a negative size
 */
double IntervalSize(Tinterval *i)
{
  double s;

  ROUNDUP;
  s=i->upper_limit-i->lower_limit;
  ROUNDNEAR;
  
  return(s);
}

/*
 * Returns the center of the interval
 */
double IntervalCenter(Tinterval *i)
{
  double c;

  if (i->upper_limit==i->lower_limit)
    c=i->upper_limit;
  else
    { 
      if ((i->lower_limit==-INF)&&(i->upper_limit==INF))
        c=0.0;
      else
        {
          c=(i->upper_limit+i->lower_limit)/2.0;
      
          /* Conditions below can occur because of roundings */
          if (c<i->lower_limit) c=i->lower_limit;
          if (c>i->upper_limit) c=i->upper_limit;
        }
    }
  return(c);
}

void SetLowerLimit(double v,Tinterval *i)
{
  i->lower_limit=v;
}

void SetUpperLimit(double v,Tinterval *i)
{
  i->upper_limit=v;
}

/*
 * Returns TRUE if i1 is fully included in i2
 */
boolean IntervalInclusion(Tinterval *i1,Tinterval *i2)
{
  return((i2->lower_limit<i1->lower_limit)&&(i1->upper_limit<i2->upper_limit));
}

/*
 * Returns TRUE if intervals 'i1' and 'i2' actually intersect.
 * NOTE: This rutine is only valid for intervals where lower_limit<upper_limit
 * and thus it fails in the intersection between certain rotational intervals.
 */
boolean Intersect(Tinterval *i1,Tinterval *i2)
{
  Tinterval i3;

  return(Intersection(i1,i2,&i3));
}

/*
 * Returns in 'i_out' the intersection of intervals 'i1' and 'i2'.
 * NOTE: This rutine is only valid for intervals where lower_limit<upper_limit
 * and thus it can fail for intersection between certain rotational intervals.
 * RETURNS true if the intersection is NOT EMPTY
 */
boolean Intersection(Tinterval *i1,Tinterval *i2,Tinterval *i_out)
{
  i_out->lower_limit=(i1->lower_limit>i2->lower_limit?i1->lower_limit:i2->lower_limit);
  i_out->upper_limit=(i1->upper_limit<i2->upper_limit?i1->upper_limit:i2->upper_limit);
  
  return(i_out->lower_limit<=i_out->upper_limit);
}

/*
 * Returns in i_out the union of i1 and i2
 * Returns false if the the input intervals are empty
*/
boolean Union(Tinterval *i1,Tinterval *i2,Tinterval *i_out)
{ 
  boolean full;

  full=TRUE;

  if (i1->lower_limit>i1->upper_limit) /*i1 empty*/
    {
      if (i2->lower_limit>i2->upper_limit)/*i2 empty*/
        full=FALSE;
      else
        {
          i_out->lower_limit=i2->lower_limit;
          i_out->upper_limit=i2->upper_limit;
        }
    }
  else
    {
      if (i2->lower_limit>i2->upper_limit) /*i2 empty*/
        {
          i_out->lower_limit=i1->lower_limit;
          i_out->upper_limit=i1->upper_limit;
        }
      else
        {  
          i_out->lower_limit=(i1->lower_limit<i2->lower_limit?i1->lower_limit:i2->lower_limit);
          i_out->upper_limit=(i1->upper_limit>i2->upper_limit?i1->upper_limit:i2->upper_limit);
        }
    }
  return(full);
}


/*
 * Returns TRUE if interval 'i' is empty
 * NOTE: This rutine is only valid for intervals where lower_limit<upper_limit
 * and thus it can return wrong results for certain rotational intervals.
 */
boolean EmptyInterval(Tinterval *i)
{
  return(i->lower_limit>i->upper_limit);
}

/*
 * returns true if p is inside the interval i
 */ 
boolean IsInside(double p,Tinterval *i)
{
  return((p>=i->lower_limit)&&(p<=i->upper_limit));
}

/*
 * Products and interval and a scalar
 * Be carefull i_out can be i1 or i2!!!!!
 */
void IntervalScale(Tinterval *i1,double e,Tinterval *i_out)
{
  double e1,e2;
      
  if (e>0)
    {
      /*The sign is preserved: extremes are keept */
      ROUNDDOWN;
      e1=i1->lower_limit*e;
      
      ROUNDUP; 
      e2=i1->upper_limit*e;
    }
  else
    {
      /*The sign is changed: extremes are swap */
      ROUNDDOWN;
      e1=i1->upper_limit*e;
      
      ROUNDUP; 
      e2=i1->lower_limit*e;
    }
  ROUNDNEAR;

  i_out->lower_limit=e1;
  i_out->upper_limit=e2;
}

/*
 * Product of two intervals
 * Be carefull i_out can be i1 or i2!!!!!
 */
void IntervalProduct(Tinterval *i1,Tinterval *i2,Tinterval *i_out)
{
  double e[3];
  unsigned int i;
  double a,b;

  ROUNDDOWN;
  e[0]=i1->lower_limit*i2->upper_limit;
  e[1]=i1->upper_limit*i2->lower_limit;
  e[2]=i1->upper_limit*i2->upper_limit;

  a=i1->lower_limit*i2->lower_limit;
  for(i=0;i<3;i++)
    if (e[i]<a) a=e[i];

  ROUNDUP; 
  e[0]=i1->lower_limit*i2->upper_limit;
  e[1]=i1->upper_limit*i2->lower_limit;
  e[2]=i1->upper_limit*i2->upper_limit;

  b=i1->lower_limit*i2->lower_limit;
  for(i=0;i<3;i++)
    if (e[i]>b) b=e[i];
   
  ROUNDNEAR;

  i_out->lower_limit=a;
  i_out->upper_limit=b;
}

/*
 * Additon of two intervals
 * Be carefull i_out can be i1 or i2!!!!!
 */
void IntervalAdd(Tinterval *i1,Tinterval *i2,Tinterval *i_out)
{
  double a,b;

  ROUNDDOWN;
  a=i1->lower_limit+i2->lower_limit;
         
  ROUNDUP;
  b=i1->upper_limit+i2->upper_limit;

  ROUNDNEAR;

  i_out->lower_limit=a;
  i_out->upper_limit=b;
}

/*
 * i_out=i1-i2;
 */
void IntervalSubstract(Tinterval *i1,Tinterval *i2,Tinterval *i_out)
{ 
  double a,b;

  ROUNDDOWN;
  a=i1->lower_limit-i2->upper_limit;
  
  ROUNDUP;
  b=i1->upper_limit-i2->lower_limit;

  ROUNDNEAR;

  i_out->lower_limit=a;
  i_out->upper_limit=b;
}

/*
 * Changes the sign of a given interval
 */
void IntervalInvert(Tinterval *i,Tinterval *i_out)
{  
  double a,b;
  
  a=-i->upper_limit;
  b=-i->lower_limit;

  i_out->lower_limit=a;
  i_out->upper_limit=b;
}

void IntervalPow(Tinterval *i,unsigned int p,Tinterval *i_out)
{
  double a,b,e1,e2;

  ROUNDDOWN;
  e1=pow(i->lower_limit,(double)p);
  e2=pow(i->upper_limit,(double)p);

  if (e1<e2)
    a=e1;
  else
    a=e2;
 
  ROUNDUP;
  e1=pow(i->lower_limit,(double)p);
  e2=pow(i->upper_limit,(double)p);

  if (e1>e2)
    b=e1;
  else
    b=e2;

  ROUNDNEAR;

  if (((p%2)==0)&&(IsInside(0,i)))
    a=0;  /*x^p with p=2*n always have a minimum at 0  (if included in the input interval)*/

  i_out->lower_limit=a;
  i_out->upper_limit=b;
}

/*
 * Square root 
 */
boolean IntervalSqrt(Tinterval *i,Tinterval *i_out)
{
  double a,b;

  if (i->upper_limit<0.0)
    return(FALSE);

  ROUNDDOWN;
  if (i->lower_limit<=0.0)
    a=0.0;
  else
    a=sqrt(i->lower_limit);

  ROUNDUP;
  b=sqrt(i->upper_limit);

  ROUNDNEAR;
  
  i_out->lower_limit=a;
  i_out->upper_limit=b;

  return(TRUE);
}

/*
*/
void IntervalDivision(Tinterval *i1,Tinterval *i2,Tinterval *i_out)
{
  double e[3];
  unsigned int i;
  double a,b;
  
  if ((i2->lower_limit<=0.0)&&(i2->upper_limit>=0.0)) 
    Error("Interval division by Zero");

  ROUNDDOWN;
  e[0]=i1->lower_limit/i2->upper_limit;
  e[1]=i1->upper_limit/i2->lower_limit;
  e[2]=i1->upper_limit/i2->upper_limit;

  a=i1->lower_limit/i2->lower_limit;
  for(i=0;i<3;i++)
    if (e[i]<a) a=e[i];

  ROUNDUP;
  e[0]=i1->lower_limit/i2->upper_limit;
  e[1]=i1->upper_limit/i2->lower_limit;
  e[2]=i1->upper_limit/i2->upper_limit;

  b=i1->lower_limit/i2->lower_limit;
  for(i=0;i<3;i++)
    if (e[i]>b) b=e[i];

  ROUNDNEAR;

  i_out->lower_limit=a;
  i_out->upper_limit=b;
}

/*
 * Shifts a interval
*/
void IntervalOffset(Tinterval *i,double offset,Tinterval *i_out)
{
  double a,b;

  ROUNDDOWN;
  a=i->lower_limit+offset;

  ROUNDUP;
  b=i->upper_limit+offset;

  ROUNDNEAR;

  i_out->lower_limit=a;
  i_out->upper_limit=b;
}


/*
 * Computes the sinus of an interval
 */
void IntervalSinus(Tinterval *i,Tinterval *i_out)
{
  double l,u;

  if (IntervalSize(i)>=M_2PI)
    NewInterval(-1,1,i_out);
  else
    {
      l=NormalizeAngle(i->lower_limit);
      u=NormalizeAngle(i->upper_limit);
  
      if (u<l)
        {
          Tinterval i_in;
          Tinterval i_out1,i_out2;

          i_in.lower_limit=l;
          i_in.upper_limit=M_2PI;

          IntervalSinus(&i_in,&i_out1);

          i_in.lower_limit=0.0;
          i_in.upper_limit=u;

          IntervalSinus(&i_in,&i_out2);

          i_out->lower_limit=(i_out1.lower_limit<i_out2.lower_limit?i_out1.lower_limit:i_out2.lower_limit);
          i_out->upper_limit=(i_out1.upper_limit>i_out2.upper_limit?i_out1.upper_limit:i_out2.upper_limit);
        }
      else
        {
          /*this is only valid for*/
          double a,b;

          ROUNDUP;
          a=sin(l);
          b=sin(u);
          ROUNDNEAR;

          if ((l<=M_PI_2)&&(M_PI_2<=u))
            i_out->upper_limit=1.0;
          else
            i_out->upper_limit=(a>b?a:b);

          ROUNDDOWN;
          a=sin(l);
          b=sin(u);
          ROUNDNEAR;

          if ((l<=(3.0*M_PI_2))&&((3.0*M_PI_2)<=u))
            i_out->lower_limit=-1.0;
          else
            i_out->lower_limit=(a<b?a:b);
        }
    }
}

/*
 * Computes the cosinus of an interval
 */
void IntervalCosinus(Tinterval *i,Tinterval *i_out)
{
  double l,u;

  if (IntervalSize(i)>=M_2PI)
    NewInterval(-1,1,i_out);
  else
    {
      l=NormalizeAngle(i->lower_limit);
      u=NormalizeAngle(i->upper_limit);
  
      if (u<l)
        {
          Tinterval i_in;
          Tinterval i_out1,i_out2;

          i_in.lower_limit=l;
          i_in.upper_limit=M_2PI;

          IntervalCosinus(&i_in,&i_out1);

          i_in.lower_limit=0.0;
          i_in.upper_limit=u;

          IntervalCosinus(&i_in,&i_out2);

          i_out->lower_limit=(i_out1.lower_limit<i_out2.lower_limit?i_out1.lower_limit:i_out2.lower_limit);
          i_out->upper_limit=(i_out1.upper_limit>i_out2.upper_limit?i_out1.upper_limit:i_out2.upper_limit);
        }
      else
        {
          double a,b;

          a=cos(l);
          b=cos(u);

          if ((l<=0.0)&&(0.0<=u))
            i_out->upper_limit=1.0;
          else
            i_out->upper_limit=(a>b?a:b);

          if ((l<=M_PI)&&(M_PI<=u))
            i_out->lower_limit=-1.0;
          else
            i_out->lower_limit=(a<b?a:b);
        }
    }
}


void IntervalAtan2(Tinterval *is,Tinterval *ic,Tinterval *i_out)
{
  double a[4];
  unsigned int i;  
  double t1,t2;

  if ((IntervalSize(is)>0.5)||(IntervalSize(ic)>0.5)) 
    Error ("Interval atan2 only works for small intervals");

  a[0]=atan2(is->lower_limit,ic->lower_limit);
  a[1]=atan2(is->lower_limit,ic->upper_limit);
  a[2]=atan2(is->upper_limit,ic->lower_limit);
  a[3]=atan2(is->upper_limit,ic->upper_limit);
    
  t1=t2=a[0];
  for(i=1;i<4;i++)
    {
      if (a[i]<t1) t1=a[i];
      else {if (a[i]>t2) t2=a[i];}
    }
  
  i_out->lower_limit=t1;
  i_out->upper_limit=t2;

  if ((i_out->upper_limit-i_out->lower_limit)>M_PI) 
    {
      double d;
      
      d=i_out->lower_limit;
      i_out->lower_limit=i_out->upper_limit;
      i_out->upper_limit=d+M_2PI;
    }
}


/*
 * Prints the contents of interval 'i' on file 'f'.
 */
void PrintInterval(FILE *f,Tinterval *i)
{
  if (EmptyInterval(i))
    fprintf(f,"(Empty Interval)");

  fprintf(f,"[");

  if (i->lower_limit==-INF) fprintf(f,"-inf");
  else fprintf(f,"%.12g",i->lower_limit);

  fprintf(f,",");

  if (i->upper_limit==+INF) fprintf(f,"+inf");
  else fprintf(f,"%.12g",i->upper_limit);

  fprintf(f,"]");
}

/*
 * Deletes the internal structures of interval 'i'.
 */
void DeleteInterval(Tinterval *i)
{
}