// ** Duesseldorf, den 20.12.1994
// ******* 
// ** interval.C
// ******

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "interval.h"
#include "fpu_ctrl.h"  // Include-Datei fuer Einstellung der Rundungsrichtung.
#include "tools.h"     // Definition von min, max.

/*----------------------------------------------------------------------------
Funktion : Show
Ausgabe des Intervalls in den String c.
----------------------------------------------------------------------------*/
char *interval::Show(char *c)
{
 sprintf(c,"[%1.16e, %1.16e]",inf,sup);
 return c;
}

/*----------------------------------------------------------------------------
Funktion : operator+=, -=, *=, /=, +, -, *, /
Definition der Operatoren auf Intervallen.
----------------------------------------------------------------------------*/
interval &interval::operator+=(interval iv)
{
 RmDn();
 inf+= iv.inf;
 RmUp();
 sup+= iv.sup;
 return *this;
}

interval &interval::operator-=(interval iv)
{
 RmDn();
 inf-= iv.sup;
 RmUp();
 sup-= iv.inf;
 return *this;
}

interval &interval::operator*=(interval iv)
{
 float _inf,a,b;

 RmDn();
 a= min(iv.inf*inf,iv.inf*sup);
 b= min(iv.sup*inf,iv.sup*sup);
 _inf= min(a,b);
 RmUp();
 a= max(iv.inf*inf,iv.inf*sup);
 b= max(iv.sup*inf,iv.sup*sup);
 sup= max(a,b);
 inf= _inf;
 return *this;
}

interval &interval::operator/=(interval iv)
{
 interval tmp(*this);

 if (!(0<<iv)) // not 0 in iv
 {
  RmDn();
  tmp.inf= 1/iv.sup;
  RmUp();
  tmp.sup= 1/iv.inf;
  tmp*= *this;
  *this= tmp;
 }
 else
 {
  printf("interval::operator/=: Division by zero.\n");
  exit(1);
 }
 return *this;
}

interval operator+(interval a, interval b)
{
 interval tmp(a);
 return tmp+=b; 
}

interval operator-(interval a, interval b)
{
 interval tmp(a);
 return tmp-=b; 
}

interval operator*(interval a, interval b)
{
 interval tmp(a);
 return tmp*=b; 
}

interval operator/(interval a, interval b)
{
 interval tmp(a);
 return tmp/=b; 
}

/*----------------------------------------------------------------------------
Funktion :Mid
Berechnet die Mitte von *this.
----------------------------------------------------------------------------*/
float interval::Mid(void)
{
 float d;
 
 d= (inf+sup)/2;
 if (!(d<<(*this))) d= inf;
 return d;
}

/*----------------------------------------------------------------------------
Funktion :operator<<
Test, ob das Intervall d in dem Intervall iv enthalten ist.
----------------------------------------------------------------------------*/
int operator<<(interval d, interval iv)
{
 return d.inf>=iv.inf && d.sup<=iv.sup;
}

/*----------------------------------------------------------------------------
Funktion : operator&&
Test, ob der Schnitt zweier Intervalle nicht leer ist.
----------------------------------------------------------------------------*/
int interval::operator&&(interval v)
{
 if (v<<*this) return 1;
 if (*this<<v) return 1;
 if (v.Inf()<<*this) return 1;
 if (v.Sup()<<*this) return 1;
 return 0;
}

/*----------------------------------------------------------------------------
Funktion : operator&
Berechnung des Schnitts zweier Intervalle.
----------------------------------------------------------------------------*/
interval interval::operator&(interval v)
{
 if (!((*this)&&v))
 {
  printf("interval::operator&: Cut is empty\n");
  exit(1);
 }
 return interval(max(inf,v.inf),min(sup,v.sup));
}