/* See ia.h */
/* Last edited on 2009-11-23 01:12 by ceaz */

#include <affirm.h>
#include <math.h>
#include <stdlib.h>
#include <ia.h>

#define ia_FULL  (ia_Full)
#define ia_ISFULL(x)  (((x).lo <= MinusInfinity) || ((x).hi >= PlusInfinity))
#define ia_NORMFULL(x)  if (ia_ISFULL(x)) (x) = ia_FULL

static Interval ia_Full;

void ia_init(void)
  { flt_init();
    ia_Full.lo = MinusInfinity; ia_Full.hi = PlusInfinity;
  }
  
Interval ia_full(void)
  { return (ia_FULL); }

int ia_is_full(Interval *x)
  { return (ia_ISFULL(*x)); }
    
Interval ia_add (Interval x, Interval y)
  { Interval z;
    //if (ia_ISFULL(x) || ia_ISFULL(y)) { return ia_FULL; }
    z.lo = x.lo + y.lo;
    z.hi = x.hi + y.hi;
    return z;
  }

Interval ia_sub(Interval x, Interval y)
  { Interval z;
    //if (ia_ISFULL(x) || ia_ISFULL(y)) { return ia_FULL; }
    z.lo = x.lo - y.hi;
    z.hi = x.hi - y.lo;
    return z;
  }

Interval ia_neg(Interval x)
  { Interval z;
    if ( ia_ISFULL(x) ) { return ia_FULL; }
    z.lo = -x.hi;
    z.hi = -x.lo;
    return z;
  }

Interval ia_scale (Interval x, Float alpha, Float zeta)
  { if ( ia_ISFULL(x) ) { return ia_FULL; }
    if (alpha == zeta)
      { return(x); }
    else if (alpha == -zeta)
      { return(ia_neg(x)); }
    else if ((zeta == Zero) || (alpha == PlusInfinity))
      { return(ia_FULL); }
    else if ((alpha == Zero) || (zeta == PlusInfinity))
      { return((Interval){Zero, Zero}); }
    else
      {	Interval z;
        double tmp;
        if ((alpha < Zero) == (zeta < Zero))
          { 
	    tmp = (double)alpha * (double)x.lo;
            z.lo = tmp/zeta ;
	    tmp = (double)alpha * (double)x.hi;
	    z.hi = tmp/zeta;
          }
        else
          { 
	    tmp = (double)alpha * (double)x.hi;
            z.lo = tmp/zeta ;
	    tmp = (double)alpha * (double)x.lo;
	    z.hi = tmp/zeta;
          }
	return z;
      }
  }
  
Interval ia_shift (Interval x, Float gamma)
  { Interval z;
    if (ia_ISFULL(x)) { return ia_FULL; }
    z.lo = x.lo + gamma;
    z.hi = x.hi + gamma;
    return z;
  }

Interval ia_mul(Interval x, Interval y)
  { Interval z;
    if (ia_ISFULL(x) || ia_ISFULL(y)) { return ia_FULL; }
    if (x.lo >= Zero)
      {	if (y.lo >= Zero)
	  {
	    z.lo = x.lo * y.lo;
	    z.hi = x.hi * y.hi;
	  }
	else if (y.hi <= Zero)
	  { 
	    z.lo = x.hi * y.lo;
	    z.hi = x.lo * y.hi;
	  }
	else
	  { 
	    z.lo = x.hi * y.lo;
	    z.hi = x.hi * y.hi;
	  }
      }
    else if (x.hi <= Zero)
      {	if (y.lo >= Zero)
	  { 
	    z.lo = x.lo * y.hi;
	    z.hi = x.hi * y.lo;
	  }
	else if (y.hi <= Zero)
	  { 
	    z.lo = x.hi * y.hi;
	    z.hi = x.lo * y.lo;
	  }
	else
	  { 
	    z.lo = x.lo * y.hi;
	    z.hi = x.lo * y.lo;
	  }
      }
    else
      {	if (y.lo >= Zero)
	  { 
	    z.lo = x.lo * y.hi;
	    z.hi = x.hi * y.hi;
	  }
	else if (y.hi <= Zero)
	  { 
	    z.lo = x.hi * y.lo;
	    z.hi = x.lo * y.lo;
	  }
	else
	  { 
	    Float z1 = x.lo * y.hi;
            Float z2 = x.hi * y.lo;
            z.lo = FMIN(z1, z2);
	    Float z3 = x.lo * y.lo;
            Float z4 = x.hi * y.hi;
            z.hi = FMAX(z3, z4);
	  }
      }
    return z;
  }

Interval ia_div (Interval x, Interval y)
  { Interval z = ia_FULL;
    if (ia_ISFULL(x) || ia_ISFULL(y)) { return z; }
    if ((y.lo <= Zero) && (y.hi >= Zero)) { return z; }
    if (x.lo >= Zero)
      {	if (y.lo >= Zero)
	  { 
	    z.lo = x.lo / y.hi;
	    z.hi = x.hi / y.lo;
	  }
	else if (y.hi <= Zero)
	  { 
	    z.lo = x.hi / y.hi;
	    z.hi = x.lo / y.lo;
	  }
      }
    else if (x.hi <= Zero)
      {	if (y.lo >= Zero)
	  { 
	    z.lo = x.lo / y.lo;
	    z.hi = x.hi / y.hi;
	  }
	else if (y.hi <= Zero)
	  {
	    z.lo = x.hi / y.lo;
	    z.hi = x.lo / y.hi;
	  }
      }
    else
      {	if (y.lo >= Zero)
	  { 
	    z.lo = x.lo / y.lo;
	    z.hi = x.hi / y.lo;
	  }
	else if (y.hi <= Zero)
	  { 
	    z.lo = x.hi / y.hi;
	    z.hi = x.lo / y.hi;
	  }
      }
    return z;
  }

Interval ia_sqrt(Interval x)
  { Interval z;
    if (ia_ISFULL(x)) { return ia_FULL; }
    if (x.hi < Zero)
      {	fatalerror("ia_sqrt: negative interval"); return ia_FULL; }
    else if (x.hi == Zero)
      {	z.hi = Zero; }
    else
      {	
	z.hi = sqrt(x.hi);
      }
    if (x.lo > Zero)
      {
	z.lo = sqrt(x.lo);
      }
    else
      {	z.lo = Zero; }
    return z;
  }


Interval ia_abs   (Interval x)
  { if (x.lo >= Zero)
      { return (x); }
    else if (x.hi <= Zero) 
      { return ((Interval){-x.hi, -x.lo}); }
    else 
      { return ((Interval){Zero, FMAX((-x.lo), x.hi)}); }
  }
  
Interval ia_max   (Interval x, Interval y)
  { return ((Interval){FMAX(x.lo, y.lo), FMAX(x.hi, y.hi)}); }
  
Interval ia_min   (Interval x, Interval y)
  { return ((Interval){FMIN(x.lo, y.lo), FMIN(x.hi, y.hi)}); }

Float ia_mid (Interval x)
  { if (ia_ISFULL(x))
      { return 0; }
    else if (x.lo == x.hi)
      { return x.lo; }
    else
      { Float m;
        m = (x.lo * Half) + (x.hi * Half);
        affirm((m >= x.lo) && (m <= x.hi), "mean failed");
        return m;
      }
  }

Float ia_mean (Interval x)
{
  return (Float)((x.lo + x.hi)/2);
}

Interval ia_pow (Interval x,double e)
{
  Interval z;
  if (ia_ISFULL(x))
      { return ; }
  else if (x.lo >= 0)
  {
    z.lo = pow(x.lo,e);
    z.hi = pow(x.hi,e);
  }
  else if (x.hi < 0)
  {
    z.lo = pow(x.hi,e);
    z.hi = pow(x.lo,e);
  }
  else
  {
    z.lo = 0;
    z.hi = pow(FMAX(-x.lo,x.hi),e);
  }
  return z;
}

Interval ia_big (Interval x, Interval y)
{
  return ((Interval){FMIN(x.lo, y.lo), FMAX(x.hi, y.hi)});
}