Islope.h

Go to the documentation of this file.

#ifndef __ISLOPE_H__
#define __ISLOPE_H__

#include "interval.h"
#include <stdlib.h>
#include <iostream>

using namespace coco;

typedef interval I;


/* This class is used to automatically compute intervalslopes of functions. */

class Islope
{
        private:
                I *data;
        
        public:
                int order;
                double *minmaxval;
                inline Islope(void);
                inline Islope(int n);
                inline Islope(double a,int n);
                inline Islope(double a, double b, int n);
                inline Islope(const Islope &a);
                virtual ~Islope();
                inline Islope& operator=(const Islope &a);
                /* Operator [] for the slope class. Returns data[pos]. */
                inline I& operator[] (int pos) const;
                inline static Islope getVar(const double &x,int n);
                inline static Islope getVar(const double &x, const double &y, int n);
                inline static Islope getVar(const I &x, int n);
                inline static Islope getVar(const double &x, const double &y, const double &z, int n);
                inline void setDeg(int n);
};

//OPERATORS
inline Islope operator +(const Islope &a, const Islope &b);
inline Islope operator +(const double &a, const Islope &b);
inline Islope operator +(const Islope &a, const double &b);
inline Islope operator -(const Islope &a, const Islope &b);
inline Islope operator -(const double &a, const Islope &b);
inline Islope operator -(const Islope &a, const double &b);
inline Islope operator -(const Islope &a);
inline Islope operator *(const Islope &a, const Islope &b);
inline Islope operator *(const double &a, const Islope &b);
inline Islope operator *(const Islope &a, const double &b);
inline Islope operator /(const Islope &a, const Islope &b);
inline Islope operator /(const double &a, const Islope &b);
inline Islope operator /(const Islope &a, const double &b);

//OUTPUT
std::ostream &operator << (std::ostream &s, const Islope &a);

//FUNCTIONS
inline Islope exp(const Islope &a);
inline Islope sin(const Islope &a);
inline Islope cos(const Islope &a);
inline Islope sqrt(const Islope &a);
inline Islope power(const Islope &a, int n);
inline Islope log(const Islope &a);
inline Islope atan(const Islope &a);
inline Islope atan2(const Islope &a, const Islope &b);


/****************/
/*INLINECODE*/
/****************/

/*Islope*/

//CONSTRUCTORS
inline Islope::Islope(void)
{                                       //all values are set to 0 in this constructor
        order=1;                        //slopeorder is 1
        data= new I[3];                 //
        data[0] = I(0.);                //=F_x
        data[1] = I(0.);                //=F_{x_c}=F_z
        data[2] = I(0.);                //=F_s=F_{s,x}
        minmaxval= new double[2];       //
        minmaxval[0] = 0.;              //=f(x_min)
        minmaxval[1] = 0.;              //=f(x_max)

}

inline Islope::Islope(int n)
{
        if ((n!=2)&&(n!=1)){
                std::cerr<<"Islope(double,int) is for initializing slope of order 1 or 2 ";
                abort();
        }
        else if(n==1){                  //all values are set to 0 in this constructor
                order=1;                        //slopeorder is 1
                data= new I[3];                 //
                data[0] = I(0.);                //=F_x
                data[1] = I(0.);                //=F_{x_c}=F_z
                data[2] = I(0.);                //=F_s=F_{s,x}
                minmaxval= new double[2];       //
                minmaxval[0] = 0.;              //=f(x_min)
                minmaxval[1] = 0.;              //=f(x_max)
        }
        else if(n==2){                          //all values are set to 0 in this constructor
                order=2;                //slopeorder is 2
                data= new I[7];         
                data[0] = I(0.);        //=F_x
                data[1] = I(0.);        //=F_{x_c}=F_z
                data[2] = I(0.);        //=F_s=F_{s,x}
                data[3] = I(0.);        //=F_y
                data[4] = I(0.);        //=f'(z)
                data[5] = I(0.);        //=F2_s
                data[6] = I(0.);        //=x
                minmaxval= new double[2];
                minmaxval[0] = 0.;      //=f(x_min)
                minmaxval[1] = 0.;      //=f(x_max)
        }
}

inline Islope::Islope(double a, int n)
{
        if ((n!=2)&&(n!=1)){
                std::cerr<<"Islope(double,int) is for initializing slope of order 1 or 2 ";
                abort();
        }
        else if(n==1){                  //the function range and the middlepoint are set to a, the slope is 0
                order=1;                //slopeorder is 1
                data= new I[3];         //
                data[0] = I(a);         //=F_x
                data[1] = I(a);         //=F_{x_c}=F_z
                data[2] = I(0.);        //=F_s=F_{s,x}
                minmaxval= new double[2];       //
                minmaxval[0] = a;       //=f(x_min)
                minmaxval[1] = a;       //=f(x_max)
        }
        else if(n==2){                  //the function range and the middlepoints are set to a, the slopes of order 1 and 2 andthe derivative at the middlepoint are 0
                order=2;                //slopeorder is 2
                data= new I[7];         
                data[0] = I(a);         //=F_x
                data[1] = I(a);         //=F_{x_c}=F_z
                data[2] = I(0.);        //=F_s=F_{s,x}
                data[3] = I(a);         //=F_y
                data[4] = I(0.);        //=f'(z)
                data[5] = I(0.);        //=F2_s
                data[6] = I(0.);        //=x
                minmaxval= new double[2];
                minmaxval[0] = a;       //=f(x_min)
                minmaxval[1] = a;       //=f(x_max)
        }
}


inline Islope::Islope(double a, double b, int n)
{
        if ((n!=2)&&(n!=1)){
                std::cerr<<"Islope(double,int) is for initializing slope of order 1 or 2 ";
                abort();
        }
        else if(n==1){                  //the function range is set to [a,b] and the middlepoint is set to (a+b)/2, the slope is 0
                order=1;                //slopeorder is 1
                data= new I[3];         //
                data[0] = I(a,b);       //=F_x
                data[1] = I((a+b)/2);   //=F_{x_c}=F_z
                data[2] = I(0.);        //=F_s=F_{s,x}
                minmaxval= new double[2];       //
                minmaxval[0] = a;       //=f(x_min)
                minmaxval[1] = b;       //=f(x_max)
        }
        else if(n==2){                  //the function range is set to [a,b] and the middlepoints are set to (a+b)/2, the slopes of order 1 and 2 andthe derivative at the middlepoint are 0
                order=2;                //slopeorder is 2
                data= new I[7];         
                data[0] = I(a,b);       //=F_x
                data[1] = I((a+b)/2);   //=F_{x_c}=F_z
                data[2] = I(0.);        //=F_s=F_{s,x}
                data[3] = I((a+b)/2);   //=F_y
                data[4] = I(0.);        //=f'(z)
                data[5] = I(0.);        //=F2_s
                data[6] = I(0.);        //=x
                minmaxval= new double[2];
                minmaxval[0] = a;       //=f(x_min)
                minmaxval[1] = b;       //=f(x_max)
        }
}

inline Islope::Islope(const Islope &d)
{
        order = d.order;
        int _tmp=3;
        if (order==2) _tmp=7;
        data = new I[_tmp];
        minmaxval= new double[2];
        memcpy(data, d.data, (_tmp)*sizeof(I));
        memcpy(minmaxval, d.minmaxval, 2*sizeof(double));
}

inline Islope& Islope::operator=(const Islope &a)
{
        if (this != &a){
                if (data) delete [] data;
                if (a.order==1){
                        data = new I[3];
                        for (int i=0; i<3; i++)
                                data[i]=a[i];
                }
                if (a.order==2){
                        data = new I[7];
                        for (int i=0; i<7; i++){
                                data[i]=a[i];
                        }
                }
                minmaxval[0]=a.minmaxval[0];
                minmaxval[1]=a.minmaxval[1];
        }
        return (*this);
}


inline I &Islope::operator[] (int pos) const
{
        return data[pos];
}



//VARIABLES

inline Islope Islope::getVar(const double &x, int n)
{
        if ((n!=2)&&(n!=1)){
                std::cerr<<"Islope(double,int) is for initializing slope of order 1 or 2 ";
                abort();
        }
        Islope d;
        if(n==1){
                Islope c;                       //set slope for function f(x)=x at specific point x
                c[0]=I(x);                      //=F_x
                c[1]=I(x);                      //=F_z
                c[2]=I(1.);                     //=F_s=F_ {s,x}
                c.minmaxval[0] = x;             //=f(x_min)
                c.minmaxval[1] = x;             //=f(x_max)
                d=c;
        }
        if(n==2){
                d.order=2;
                Islope c(2);                    //set slope of order 2 for function f(x)=x at specific point x.
                c[0] = I(x);                    //=F_x
                c[1] = I(x);                    //=F_{x_c}=F_z
                c[2] = I(1.);                   //=F_s=F_{s,x}
                c[3] = I(x);                    //=F_y
                c[4] = I(1.);                   //=f'(z)
                c[5] = I(0.);                   //=F2_s
                c[6] = I(x);                    //=x
                c.minmaxval[0] = x;             //=f(x_min)
                c.minmaxval[1] = x;             //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope Islope::getVar(const double &x, const double &y, int n)
{
        if ((n!=2)&&(n!=1)){
                std::cerr<<"Islope(double,int) is for initializing slope of order 1 or 2 ";
                abort();
        }
        Islope d;
        if(n==1){
                Islope c;                       //set slope for function f(x)=x in specific interval [x,y]. Middlepoint is (x+y)/2
                c[0]=I(x,y);                    //=F_x
                c[1]=I((x+y)/2);                //=F_z
                c[2]=I(1.);                     //=F_s=F_{s,x}
                c.minmaxval[0] = x;             //=f(x_min)
                c.minmaxval[1] = y;             //=f(x_max)
                d=c;
        }
        if (n==2){
                d.order=2;
                Islope c(2);                    //set slope of order 2 for function f(x)=x in specific interval [x,y]. Middlepoints are (x+y)/2
                c[0] = I(x,y);                  //=F_x
                c[1] = I((x+y)/2);              //=F_{x_c}=F_z
                c[2] = I(1.);                   //=F_s=F_{s,x}
                c[3] = I((x+y)/2);              //=F_y
                c[4] = I(1.);                   //=f'(z)
                c[5] = I(0.);                   //=F2_s
                c[6] = I(x,y);                  //=x
                c.minmaxval[0] = x;             //=f(x_min)
                c.minmaxval[1] = y;             //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope Islope::getVar(const I &x, int n)
{
        if ((n!=2)&&(n!=1)){
                std::cerr<<"Islope(double,int) is for initializing slope of order 1 or 2 ";
                abort();
        }
        Islope d;
        if(n==1){
                Islope c;                       //set slope for function f(x)=x in specific interval [x,y]. Middlepoint is (x+y)/2
                c[0]=x;                         //=F_x
                c[1]=I((x.inf()+x.sup())/2);    //=F_z
                c[2]=I(1.);                     //=F_s=F_{s,x}
                c.minmaxval[0] = x.inf();       //=f(x_min)
                c.minmaxval[1] = x.sup();       //=f(x_max)
                d=c;
        }
        if (n==2){
                d.order=2;
                Islope c(2);                    //set slope of order 2 for function f(x)=x in specific interval [x,y]. Middlepoints are (x+y)/2
                c[0] = x;                       //=F_x
                c[1] = I((x.inf()+x.sup())/2);  //=F_{x_c}=F_z
                c[2] = I(1.);                   //=F_s=F_{s,x}
                c[3] = I((x.inf()+x.sup())/2);  //=F_y
                c[4] = I(1.);                   //=f'(z)
                c[5] = I(0.);                   //=F2_s
                c[6] = x;                       //=x
                c.minmaxval[0] = x.inf();       //=f(x_min)
                c.minmaxval[1] = x.sup();       //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope Islope::getVar(const double &x, const double &y, const double &z, int n)
{
        if ((n!=2)&&(n!=1)){
                std::cerr<<"Islope(double,int) is for initializing slope of order 1 or 2 ";
                abort();
        }
        Islope d;
        if(n==1){
                Islope c;                       //set slope for function f(x)=x in specific interval [x,y]. Middlepoint is (x+y)/2
                c[0]=I(x,y);                    //=F_x
                c[1]=I(z);                      //=F_z
                c[2]=I(1.);                     //=F_s=F_{s,x}
                c.minmaxval[0] = x;             //=f(x_min)
                c.minmaxval[1] = y;             //=f(x_max)
                d=c;
        }
        if (n==2){
                d.order=2;
                Islope c(2);                    //set slope of order 2 for function f(x)=x in specific interval [x,y]. Middlepoints are z
                c[0] = I(x,y);                  //=F_x
                c[1] = I(z);                    //=F_{x_c}=F_z
                c[2] = I(1.);                   //=F_s=F_{s,x}
                c[3] = I(z);                    //=F_y
                c[4] = I(1.);                   //=f'(z)
                c[5] = I(0.);                   //=F2_s
                c[6] = I(x,y);                  //=x
                c.minmaxval[0] = x;             //=f(x_min)
                c.minmaxval[1] = y;             //=f(x_max)
                d=c;
        }
        return d;
}

inline void Islope::setDeg(int n)
{
        I tmp[7];
        if ((n==1)&&(order==2)){        //
                for (int i=0; i<3; i++) //
                        tmp[i]=data[i]; //save data[0]-data[2]
                delete [] data;         //delete all data
                data=new I[3];          //initialize data new
                for (int i=0; i<3; i++) //
                        data[i]=tmp[i]; //copy saved data
        }
        if ((n==2)&&(order==1)) {
                for (int i=0; i<3; i++) //
                        tmp[i]=data[i]; //save data[0]-data[2]
                delete [] data;         //delete all data
                data=new I[7];          //initialize data new
                for (int i=0; i<3; i++)
                        data[i]=tmp[i]; //copy saved data
                for (int j=3; j<7; j++)
                        data[j]=0.;     //initialize data[3]-data[6]
        }
        order=n;
}


//OPERATORS
inline Islope operator +(const Islope &a, const Islope &b)
{
        if (a.order!=b.order){
                std::cerr<<"+ slope order should be the same";
                abort();
        }
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=a[0]+b[0];                                 //=F_x
                c[1]=a[1]+b[1];                                 //=F_z
                c[2]=a[2]+b[2];                                 //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=a[0]+b[0];                                 //=F_x
                c[1]=a[1]+b[1];                                 //=F_z
                c[2]=a[2]+b[2];                                 //=F_s
                c[3]=a[3]+b[3];                                 //=F_y
                c[4]=a[4]+b[4];                                 //=f'(z)
                c[5]=a[5]+b[5];                                 //=F2_s
                if ((a[6].inf()==0) && (a[6].sup()==0)){
                        c[6]=b[6];                              //=x
                }
                else {
                        c[6]=a[6];                              //=x
                }
                c.minmaxval[0]=a.minmaxval[0]+b.minmaxval[0];   //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]+b.minmaxval[1];   //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator +(const double &a, const Islope &b)
{
        Islope d;
        if (b.order==1){
                Islope c;
                c[0]=I(a)+b[0];                                 //=F_x
                c[1]=I(a)+b[1];                                 //=F_z
                c[2]=b[2];                                      //=F_s
                d=c;
        }
        if (b.order==2){
                Islope c(2);
                d.order=2;
                c[0]=I(a)+b[0];                                 //=F_x
                c[1]=I(a)+b[1];                                 //=F_z
                c[2]=b[2];                                      //=F_s
                c[3]=I(a)+b[3];                                 //=F_y
                c[4]=b[4];                                      //=f'(z)
                c[5]=b[5];                                      //=F2_s
                c[6]=b[6];                                      //=x
                c.minmaxval[0]=a+b.minmaxval[0];                //=f(x_min)
                c.minmaxval[1]=a+b.minmaxval[1];                //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator +(const Islope &a, const double &b)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=a[0]+I(b);                                 //=F_x
                c[1]=a[1]+I(b);                                 //=F_z
                c[2]=a[2];                                      //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=a[0]+I(b);                                 //=F_x
                c[1]=a[1]+I(b);                                 //=F_z
                c[2]=a[2];                                      //=F_s
                c[3]=a[3]+I(b);                                 //=F_y
                c[4]=a[4];                                      //=f'(z)
                c[5]=a[5];                                      //=F2_s
                c[6]=a[6];                                      //=x
                c.minmaxval[0]=a.minmaxval[0]+b;                //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]+b;                //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator -(const Islope &a, const Islope &b)
{
        if (a.order!=b.order){
                std::cerr<<"+ slope order should be the same";
                abort();
        }
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=a[0]-b[0];                                 //=F_x
                c[1]=a[1]-b[1];                                 //=F_z
                c[2]=a[2]-b[2];                                 //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=a[0]-b[0];                                 //=F_x
                c[1]=a[1]-b[1];                                 //=F_z
                c[2]=a[2]-b[2];                                 //=F_s
                c[3]=a[3]-b[3];                                 //=F_y
                c[4]=a[4]-b[4];                                 //=f'(z)
                c[5]=a[5]-b[5];                                 //=F2_s
                if ((a[6].inf()==0) && (a[6].sup()==0)){
                        c[6]=b[6];                              //=x
                }
                else {
                        c[6]=a[6];                              //=x
                }
                c.minmaxval[0]=a.minmaxval[0]-b.minmaxval[0];   //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]-b.minmaxval[1];   //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator -(const double &a, const Islope &b)
{
        Islope d;
        if (b.order==1){
                Islope c;
                c[0]=I(a)-b[0];                                 //=F_x
                c[1]=I(a)-b[1];                                 //=F_z
                c[2]=-b[2];                                     //=F_s
                d=c;
        }
        if (b.order==2){
                Islope c(2);
                d.order=2;
                c[0]=I(a)-b[0];                                 //=F_x
                c[1]=I(a)-b[1];                                 //=F_z
                c[2]=-b[2];                                     //=F_s
                c[3]=I(a)-b[3];                                 //=F_y
                c[4]=-b[4];                                     //=f'(z)
                c[5]=-b[5];                                     //=F2_s
                c[6]=b[6];                                      //=x
                c.minmaxval[0]=a-b.minmaxval[0];                //=f(x_min)
                c.minmaxval[1]=a-b.minmaxval[1];                //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator -(const Islope &a, const double &b)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=a[0]-I(b);                                 //=F_x
                c[1]=a[1]-I(b);                                 //=F_z
                c[2]=a[2];                                      //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=a[0]-I(b);                                 //=F_x
                c[1]=a[1]-I(b);                                 //=F_z
                c[2]=a[2];                                      //=F_s
                c[3]=a[3]-I(b);                                 //=F_y
                c[4]=a[4];                                      //=f'(z)
                c[5]=a[5];                                      //=F2_s
                c[6]=a[6];                                      //=x
                c.minmaxval[0]=a.minmaxval[0]+b;                //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]+b;                //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator -(const Islope &a)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=-a[0];                                     //=F_x
                c[1]=-a[1];                                     //=F_z
                c[2]=-a[2];                                     //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=-a[0];                                     //=F_x
                c[1]=-a[1];                                     //=F_z
                c[2]=-a[2];                                     //=F_s
                c[3]=-a[3];                                     //=F_y
                c[4]=-a[4];                                     //=f'(z)
                c[5]=-a[5];                                     //=F2_s
                c[6]=a[6];                                      //=x
                c.minmaxval[0]=-a.minmaxval[0];                 //=f(x_min)
                c.minmaxval[1]=-a.minmaxval[1];                 //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator *(const Islope &a, const Islope &b)
{
        if (a.order!=b.order){
                std::cerr<<"* slope order should be the same";
                abort();
        }
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=a[0]*b[0];                                 //=F_x
                c[1]=a[1]*b[1];                                 //=F_z
                c[2]=a[2]*b[1]+a[0]*b[2];                       //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=a[0]*b[0];                                 //=F_x
                c[1]=a[1]*b[1];                                 //=F_z
                c[2]=a[2]*b[1]+a[0]*b[2];                       //=F_s
                c[3]=a[3]*b[3];                                 //=F_y
                c[4]=a[4]*b[1]+a[1]*b[4];                       //=f'(z)
                c[5]=a[5]*b[1]+a[2]*b[2]+a[3]*b[5];             //=F2_s
                if ((a[6].inf()==0) && (a[6].sup()==0)){        
                        c[6]=b[6];                              //=x
                }
                else {
                        c[6]=a[6];                              //=x
                }
                c.minmaxval[0]=a.minmaxval[0]*b.minmaxval[0];   //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]*b.minmaxval[1];   //=f(x_max)
                d=c;
        }       
        return d;
}

inline Islope operator *(const double &a, const Islope &b)
{
        Islope d;
        if (b.order==1){
                Islope c;
                c[0]=I(a)*b[0];                                 //=F_x
                c[1]=I(a)*b[1];                                 //=F_z
                c[2]=I(a)*b[2];                                 //=F_s
                d=c;
        }
        if (b.order==2){
                Islope c(2);
                d.order=2;
                c[0]=I(a)*b[0];                                 //=F_x
                c[1]=I(a)*b[1];                                 //=F_z
                c[2]=I(a)*b[2];                                 //=F_s
                c[3]=I(a)*b[3];                                 //=F_y
                c[4]=I(a)*b[4];                                 //=f'(z)
                c[5]=I(a)*b[5];                                 //=F2_s
                c[6]=b[6];                                      //=x
                c.minmaxval[0]=a*b.minmaxval[0];                //=f(x_min)
                c.minmaxval[1]=a*b.minmaxval[1];                //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator *(const Islope &a, const double &b)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=I(b)*a[0];                                 //=F_x
                c[1]=I(b)*a[1];                                 //=F_z
                c[2]=I(b)*a[2];                                 //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=I(b)*a[0];                                 //=F_x
                c[1]=I(b)*a[1];                                 //=F_z
                c[2]=I(b)*a[2];                                 //=F_s
                c[3]=I(b)*a[3];                                 //=F_y
                c[4]=I(b)*a[4];                                 //=f'(z)
                c[5]=I(b)*a[5];                                 //=F2_s
                c[6]=a[6];                                      //=x
                c.minmaxval[0]=a.minmaxval[0]*b;                //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]*b;                //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope operator /(const double &a, const Islope &b)
{
        Islope d;
        if (b.order==1){
                Islope c(a,1);                  //make a slope of order 1 out of the double a
                d=c/b;                          //use the slope division operator
        }
        if (b.order==2){
                Islope c(a,2);                  //make a slope of order 2 out of the double a
                d.order=2;                      
                d=c/b;                          //use the slope division operator
        }
        return d;
}

inline Islope operator /(const Islope &a, const Islope &b)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=a[0]/b[0];                 //=F_x
                c[1]=a[1]/b[1];                 //=F_z
                c[2]=(a[2]-c[1]*b[2])/b[0];     //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=a[0]/b[0];                                 //=F_x
                c[1]=a[1]/b[1];                                 //=F_z
                c[2]=(a[2]-c[1]*b[2])/b[0];                     //=F_s
                c[3]=a[3]/b[3];                                 //=F_y
                c[4]=(a[4]*b[1]-a[1]*b[4])/power(b[1],2);       //=f'(z)
                c[5]=a[5]/b[0]-c[1]*b[5]/b[0]+(c[1]*b[2]*(b[2]-a[2]))/(b[0]*b[3]);      //=F2_s
                if ((a[6].inf()==0) && (a[6].sup()==0)){        
                        c[6]=b[6];                              //=x
                }
                else {
                        c[6]=a[6];                              //=x
                }
                c.minmaxval[0]=a.minmaxval[0]/b.minmaxval[0];   //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]/b.minmaxval[1];   //=f(x_max)
                d=c;
                }
        return d;
}

inline Islope operator /(const Islope &a, const double &b)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=a[0]/I(b);                         //=F_x
                c[1]=a[1]/I(b);                         //=F_z
                c[2]=a[2]/I(b);                         //=F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=a[0]/I(b);                         //=F_x
                c[1]=a[1]/I(b);                         //=F_z
                c[2]=a[2]/I(b);                         //=F_s
                c[3]=a[3]/I(b);                         //=F_y
                c[4]=a[4]/I(b);                         //=f'(z)
                c[5]=a[5]/I(b);                         //=F2_s
                c[6]=a[6];                              //=x
                c.minmaxval[0]=a.minmaxval[0]/b;        //=f(x_min)
                c.minmaxval[1]=a.minmaxval[1]/b;        //=f(x_max)     
                d=c;
        }
        return d;
}

inline Islope exp(const Islope &a)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=exp(a[0]);                                                 //=F_x
                c[1]=exp(a[1]);                                                 //=F_z
                double lowdiff=a[0].inf()-a[1].inf();                           //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                           //=x_max - z
                if ((lowdiff==0)&&(supdiff==0)){                                //if lowdiff and supdiff is 0 the derivative has to be used
                        c[2]=c[0]*a[2];                                         //because division by lowdiff or supdiff is not possible
                }
                if ((lowdiff!=0)&&(supdiff!=0)){
                        double low=(exp(a[0].inf())-exp(a[1].inf()))/lowdiff;   //use formula from kolev (exp is convex everywhere)
                        double sup=(exp(a[0].sup())-exp(a[1].sup()))/supdiff;   //
                        c[2]=I(low,sup)*a[2];                                   //
                }
                if ((lowdiff==0)&&(supdiff!=0)){
                        double low=exp(a[0].inf());                             //use derivative (lowdiff is 0)
                        double sup=(exp(a[0].sup())-exp(a[1].sup()))/supdiff;   //use formula from kolev (exp is convex everywhere)
                        c[2]=I(low,sup)*a[2];                                   //
                }
                if ((lowdiff!=0)&&(supdiff==0)){
                        double low=(exp(a[0].inf())-exp(a[1].inf()))/lowdiff;   //use formula from kolev (exp is convex everywhere)
                        double sup=exp(a[0].sup());                             //use derivative (supdiff is 0)
                        c[2]=I(low,sup)*a[2];                                   //
                }
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=exp(a[0]);                                                 //=F_x
                c[1]=exp(a[1]);                                                 //=F_z
                c[3]=exp(a[3]);                                                 //=F_y
                c[6]=a[6];                                                      //=x
                double lowdiff=a[0].inf()-a[1].inf();                           //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                           //=x_max - z
                if ((lowdiff==0)&&(supdiff==0)){                                //if lowdiff or supdiff is 0 the derivative has to be used because division by lowdiff or supdiff is not possible
                        c[2]=c[0]*a[2];                                         //=F_s
                        c[4]=c[1]*a[4];                                         //=f'(z)
                        c[5]=c[0]*(0.5*power(a[2],2)+a[5]);                     //=F2_s
                }
                if ((lowdiff!=0)&&(supdiff!=0)){
                        double low=(exp(a[0].inf())-exp(a[1].inf()))/lowdiff;   //use formula from kolev (exp is convex everywhere)
                        double sup=(exp(a[0].sup())-exp(a[1].sup()))/supdiff;   //
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[4]=c[1]*a[4];                                         //=f'(z)
                        double low2=(low-c[1].inf())/lowdiff;                   //use formula from kolev for second order slope (exp' is convex everywhere)
                        double sup2=(sup-c[1].sup())/supdiff;                   //
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                if ((lowdiff==0)&&(supdiff!=0)){
                        double low=exp(a[0].inf());                             //use derivative
                        double sup=(exp(a[0].sup())-exp(a[1].sup()))/supdiff;   //use formula from kolev (exp is convex everywhere)
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[4]=c[1]*a[4];                                         //=f'(z)
                        double low2=0.5*exp(a[0].inf());                        //use derivative
                        double sup2=(sup-c[1].sup())/supdiff;                   //use formula from kolev for second order slope (exp' is convex everywhere)
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                if ((lowdiff!=0)&&(supdiff==0)){
                        double low=(exp(a[0].inf())-exp(a[1].inf()))/lowdiff;   //use formula from kolev (exp is convex everywhere)
                        double sup=exp(a[0].sup());                             //use derivative
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[4]=c[1]*a[4];                                         //=f'(z)
                        double low2=(low-c[1].inf())/lowdiff;                   //use formula from kolev for second order slope (exp' is convex everywhere)
                        double sup2=0.5*exp(a[0].sup());                        //use derivative
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                if (c[5].empty()){                                              //if roundingerrors are to big, it can happen that low2 > sup2, so c[5] is empty. Then the derivative is used instead of the kolev formula.
                        c[2]=c[0]*a[2];                                         //=F_s
                        c[5]=c[0]*(0.5*power(a[2],2)+a[5]);                     //=F2_s
                }
                c.minmaxval[0]=exp(a.minmaxval[0]);                             //=f(x_min)
                c.minmaxval[1]=exp(a.minmaxval[1]);                             //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope sin(const Islope &a)
{
        Islope d;
        double low;
        double sup;
        double low2;
        double sup2;
        double lowdiff;
        double supdiff;
        
        if (a.order==1){
                Islope c;
                c[0]=sin(a[0]);                                                 //=F_x
                c[0]=sin(a[0]);                                                 //muss da stehen sonst ist c[0]=[empty]
                c[1]=sin(a[1]);                                                 //=F_z
                if (width(a[0])>=U_PI){                                         //if distance between points is greater than pi, then sin in this interval is neither convex nor concave and the derivative has to be used.
                        c[2]=cos(a[0])*a[2];                                    //=F_s
                }
                
                else if ((c[0].inf()<=0) && (c[0].sup()<=0)){                   //if both functionvalues are negative sin is convex in the interval between the points
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                low=(sin(a[0].inf())-sin(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                low=cos(a[0].inf());                            //if lowdiff is 0 the derivative has to be used
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                sup=(sin(a[0].sup())-sin(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                sup=cos(a[0].sup());                            //if supdiff is 0 the derivative has to be used
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                }
                
                else if ((c[0].inf()>=0) && (c[0].sup()>=0)){                   //if both functionvalues are positive sin is concave in the interval between the points
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                sup=(sin(a[0].inf())-sin(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                sup=cos(a[0].inf());                            //if lowdiff is 0 the derivative has to be used
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                low=(sin(a[0].sup())-sin(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                low=cos(a[0].sup());                            //if supdiff is 0 the derivative has to be used
                        }
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                }
                
                else {                                                          //if the functionvalues are not both positive or negative sin is neither convex nor concav and the derivative has to be used
                        c[2]=cos(a[0])*a[2];                                    //=F_s
                }
                
                d=c;
        }
        
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=sin(a[0]);                                                 //=F_x
                c[0]=sin(a[0]);                                                 //muss da stehen sonst ist c[0]=[empty]
                c[1]=sin(a[1]);                                                 //=F_z
                c[3]=sin(a[3]);                                                 //=F_y
                c[4]=cos(a[1])*a[4];                                            //=f'(z)
                c[6]=a[6];                                                      //=x
                I tmp=cos(a[0]);                                                //outer derivative of sin(a[0])
                if (width(a[0])>=U_PI/2.){                                      //if the distance between the points is greater than pi/2, sin and sin' are not convex or concave and the derivative has to be used
                        c[2]=cos(a[0])*a[2];                                    //=F_s
                        c[5]=0.5*(-sin(a[0])*power(a[2],2))+cos(a[0])*a[5];     //=F2_s
                }
                
                else if ((c[0].inf()>=0) && (c[0].sup()>=0) && (tmp.inf()>=0) && (tmp.sup()>=0)){       //if all functionvalues are positiv sin and sin' are concave in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                sup=(sin(a[0].inf())-sin(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                sup2=(sup-cos(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                sup=cos(a[0].inf());
                                sup2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                low=(sin(a[0].sup())-sin(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                low2=(low-cos(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                low=cos(a[0].sup());
                                low2=-0.5*c[0].sup();
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                
                else if ((c[0].inf()>=0) && (c[0].sup()>=0) && (tmp.inf()<=0) && (tmp.sup()<=0)){       //if functionvalues for sin are positiv and for cos negative sin is concave and sin' is convex in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z 
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                sup=(sin(a[0].inf())-sin(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                low2=(sup-cos(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                sup=cos(a[0].inf());
                                low2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                low=(sin(a[0].sup())-sin(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                sup2=(low-cos(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                low=cos(a[0].sup());
                                sup2=-0.5*c[0].sup();
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                        
                        //c[2]=c[2].intersect(cos(a[0])*a[2]);
                        //c[5]=c[5].intersect(0.5*(-sin(a[0])*power(a[2],2)+cos(a[0])*a[5]));
                }
                
                else if ((c[0].inf()<=0) && (c[0].sup()<=0) && (tmp.inf()<=0) && (tmp.sup()<=0)){       //if functionvalues for sin are negative and for cos negative sin is convex and sin' is convex in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                low=(sin(a[0].inf())-sin(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                low2=(low-cos(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                low=cos(a[0].inf());
                                low2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                sup=(sin(a[0].sup())-sin(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                sup2=(sup-cos(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                sup=cos(a[0].sup());
                                sup2=-0.5*c[0].sup();
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                        
                        //c[2]=c[2].intersect(cos(a[0])*a[2]);
                        //c[5]=c[5].intersect(0.5*(-sin(a[0])*power(a[2],2)+cos(a[0])*a[5]));
                }
                
                else if ((c[0].inf()<=0) && (c[0].sup()<=0) && (tmp.inf()>=0) && (tmp.sup()>=0)){       //if functionvalues for sin are negative and for cos positive sin is convex and sin' is concave in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                low=(sin(a[0].inf())-sin(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                sup2=(low-cos(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                low=cos(a[0].inf());
                                sup2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                sup=(sin(a[0].sup())-sin(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                low2=(sup-cos(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                sup=cos(a[0].sup());                            
                                low2=-0.5*c[0].sup();                           
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                        
                        //c[2]=c[2].intersect(cos(a[0])*a[2]);
                        //c[5]=c[5].intersect(0.5*(-sin(a[0])*power(a[2],2)+cos(a[0])*a[5]));
                }
                
                else {                                                          //if sin and sin' are not concave or convex the derivative has to be used instead of the kolev formula
                        c[2]=cos(a[0])*a[2];                                    //=F_s
                        c[5]=0.5*(-sin(a[0])*power(a[2],2))+cos(a[0])*a[5];     //=F2_s
                }
                
                if (c[5].empty()) {                                             //if roundingerrors are to big, it can happen that low2 > sup2, so c[5] is empty. Then the derivative is used instead of the kolev formula.
                        c[2]=cos(a[0])*a[2];                                    //=F_s
                        c[5]=0.5*(-sin(a[0])*power(a[2],2))+cos(a[0])*a[5];     //=F2_s
                }
                
                c.minmaxval[0]=sin(a.minmaxval[0]);                             //=f(x_min)
                c.minmaxval[1]=sin(a.minmaxval[1]);                             //=f(x_max)
                
                d=c;
        }       
        return d;
}

inline Islope cos(const Islope &a)
{
        Islope d;
        double low;
        double sup;
        double low2;
        double sup2;
        double lowdiff;
        double supdiff;
        
        if (a.order==1){
                Islope c;
                c[0]=cos(a[0]);                                                 //=F_x
                c[0]=cos(a[0]);                                                 //muss da stehen sonst ist c[0]=[empty]
                c[1]=cos(a[1]);                                                 //=F_z
                if (width(a[0])>=U_PI){                                         //if distance between points is greater than pi, then cos in this interval is neither convex nor concave and the derivative has to be used.
                        c[2]=-sin(a[0])*a[2];                                   //=F_s
                }
                
                else if ((c[0].inf()<=0) && (c[0].sup()<=0)){                   //if both functionvalues are negative cos is convex in the interval between the points
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                low=(cos(a[0].inf())-cos(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                low=-sin(a[0].inf());                           //if lowdiff is 0 the derivative has to be used
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                sup=(cos(a[0].sup())-cos(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                sup=-sin(a[0].sup());                           //if supdiff is 0 the derivative has to be used
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                }
                
                else if ((c[0].inf()>=0) && (c[0].sup()>=0)){                   //if both functionvalues are positive cos is concave in the interval between the points
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                sup=(cos(a[0].inf())-cos(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                sup=-sin(a[0].inf());                           //if lowdiff is 0 the derivative has to be used
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                low=(cos(a[0].sup())-cos(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                        }
                        else {
                                low=-sin(a[0].sup());                           //if supdiff is 0 the derivative has to be used
                        }
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                }
                
                else {                                                          //if the functionvalues are not both positive or negative cos is neither convex nor concav and the derivative has to be used
                        c[2]=-sin(a[0])*a[2];                                   //=F_s
                }
                
                d=c;
        }
        
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=cos(a[0]);                                                 //=F_x
                c[0]=cos(a[0]);                                                 //muss da stehen sonst ist c[0]=[empty]
                c[1]=cos(a[1]);                                                 //=F_z
                c[3]=cos(a[3]);                                                 //=F_y
                c[4]=-sin(a[1])*a[4];                                           //=f'(z)
                c[6]=a[6];                                                      //=x
                I tmp=sin(a[0]);                                                //outer derivative of cos(a[0])
                if (width(a[0])>=U_PI/2.){                                      //if the distance between the points is greater than pi/2, cos and cos' are not convex or concave and the derivative has to be used
                        c[2]=cos(a[0])*a[2];                                    //=F_s
                        c[5]=0.5*(-cos(a[0])*power(a[2],2))-sin(a[0])*a[5];     //=F2_s
                }
                
                else if ((c[0].inf()>=0) && (c[0].sup()>=0) && (tmp.inf()>=0) && (tmp.sup()>=0)){       //if all functionvalues are positiv cos and cos' are concave in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                sup=(cos(a[0].inf())-cos(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                sup2=(sup+sin(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                sup=-sin(a[0].inf());
                                sup2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                low=(cos(a[0].sup())-cos(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                low2=(low+sin(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                low=-sin(a[0].sup());
                                low2=-0.5*c[0].sup();
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                
                else if ((c[0].inf()>=0) && (c[0].sup()>=0) && (tmp.inf()<=0) && (tmp.sup()<=0)){       //if functionvalues for cos are positiv and for sin negative cos is concave and cos' is convex in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z 
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                sup=(cos(a[0].inf())-cos(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                low2=(sup+sin(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                sup=-sin(a[0].inf());
                                low2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                low=(cos(a[0].sup())-cos(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                sup2=(low+sin(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                low=-sin(a[0].sup());
                                sup2=-0.5*c[0].sup();
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                        
                        //c[2]=c[2].intersect(cos(a[0])*a[2]);
                        //c[5]=c[5].intersect(0.5*(-sin(a[0])*power(a[2],2)+cos(a[0])*a[5]));
                }
                
                else if ((c[0].inf()<=0) && (c[0].sup()<=0) && (tmp.inf()<=0) && (tmp.sup()<=0)){       //if functionvalues for cos are negative and for sin negative cos is convex and cos' is convex in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                low=(cos(a[0].inf())-cos(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                low2=(low+sin(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                low=-sin(a[0].inf());
                                low2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                sup=(cos(a[0].sup())-cos(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                sup2=(sup+sin(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                sup=-sin(a[0].sup());
                                sup2=-0.5*c[0].sup();
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                        
                        //c[2]=c[2].intersect(cos(a[0])*a[2]);
                        //c[5]=c[5].intersect(0.5*(-sin(a[0])*power(a[2],2)+cos(a[0])*a[5]));
                }
                
                else if ((c[0].inf()<=0) && (c[0].sup()<=0) && (tmp.inf()>=0) && (tmp.sup()>=0)){       //if functionvalues for cos are negative and for sin positive cos is convex and cos' is concave in the given intervall
                        lowdiff=a[0].inf()-a[1].inf();                          //=x_min - z
                        supdiff=a[0].sup()-a[1].sup();                          //=x_max - z
                        
                        if (lowdiff!=0){                                        //if lowdiff is not 0 kolev formula can be used
                                low=(cos(a[0].inf())-cos(a[1].inf()))/lowdiff;  //use formula from kolev for first order slopes
                                sup2=(low+sin(a[1].inf()))/lowdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if lowdiff is 0 the derivative has to be used
                                low=-sin(a[0].inf());
                                sup2=-0.5*c[0].inf();
                        }
                        
                        if (supdiff!=0){                                        //if supdiff is not 0 kolev formula can be used
                                sup=(cos(a[0].sup())-cos(a[1].sup()))/supdiff;  //use formula from kolev for first order slopes
                                low2=(sup+sin(a[1].sup()))/supdiff;             //use formula from kolev for second order slopes
                        }
                        else {                                                  //if supdiff is 0 the derivative has to be used
                                sup=-sin(a[0].sup());                           
                                low2=-0.5*c[0].sup();                           
                        }
                        
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                        
                        //c[2]=c[2].intersect(cos(a[0])*a[2]);
                        //c[5]=c[5].intersect(0.5*(-sin(a[0])*power(a[2],2)+cos(a[0])*a[5]));
                }
                
                else {                                                          //if cos and cos' are not concave or convex the derivative has to be used instead of the kolev formula
                        c[2]=-sin(a[0])*a[2];                                   //=F_s
                        c[5]=0.5*(-cos(a[0])*power(a[2],2))-sin(a[0])*a[5];     //=F2_s
                }
                
                if (c[5].empty()) {                                             //if roundingerrors are to big, it can happen that low2 > sup2, so c[5] is empty. Then the derivative is used instead of the kolev formula.
                        c[2]=-sin(a[0])*a[2];                                   //=F_s
                        c[5]=0.5*(-cos(a[0])*power(a[2],2))-sin(a[0])*a[5];     //=F2_s
                }
                
                c.minmaxval[0]=cos(a.minmaxval[0]);                             //=f(x_min)
                c.minmaxval[1]=cos(a.minmaxval[1]);                             //=f(x_max)
                
                d=c;
        }       
        return d;
}

inline Islope sqrt(const Islope &a)
{
        if (a[0].inf()<0){
                std::cerr<<"sqrt: interval contains negative numbers";
                abort();
        }
        Islope d;
        
        if (a.order==1){
                Islope c;
                c[0]=sqrt(a[0]);                                                //=F_x
                c[1]=sqrt(a[1]);                                                //=F_z
                c[2]=a[2]/(c[0]+c[1]);                                          //=F_s
                d=c;
        }
                
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=sqrt(a[0]);                                                //=F_x
                c[1]=sqrt(a[1]);                                                //=F_z
                c[3]=sqrt(a[3]);                                                //=F_y
                c[6]=a[6];                                                      //=x
                /*double lowdiff=a[0].inf()-a[1].inf();
                double supdiff=a[0].sup()-a[1].sup();
                if ((lowdiff==0)||(supdiff==0)){
                        c[2]=a[2]/(c[0]+c[1]);
                        c[4]=a[4]/(2.*(c[1]));
                //      c[5]=;
                }*/
                
                c[2]=a[2]/(c[0]+c[1]);                                          //=F_s
                c[4]=a[4]/(2.*c[1]);                                            //=f'(z)
                c[5]=a[5]/(c[0]+c[1])-a[4]*a[2]/(power(c[0]+c[1],2)*4.*c[1]);   //=F2_s
                
                c.minmaxval[0]=sqrt(a.minmaxval[0]);                            //=f(x_min)
                c.minmaxval[1]=sqrt(a.minmaxval[1]);                            //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope power(const Islope &a, int n)
{
        if (n<2){
                std::cerr<<"power degree not supported: "<< n;
                abort();
        }
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=power(a[0], n);                                                                    //=F_x
                c[1]=power(a[1], n);                                                                    //=F_z
                double lowdiff=a[0].inf()-a[1].inf();                                                   //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                                                   //=x_max - z
                switch (n % 2) {                                                                        //if n even power(x,n) is convex
                        case 0: if (n == 2) {                                                           //if n is 2 a direct formula can be used
                                                c[2]=a[2] * (a[0] + a[1]);                              //=F_s
                                }
                                else if ((lowdiff != 0) && (supdiff != 0)) {                            //if lowidff and supdiff are both not 0 the kolev formula can be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else if ((lowdiff == 0) && (supdiff != 0)) {                            //if lowdiff is 0 the derivative at x_min has to be used
                                        double low=n*pow(a[0].inf(),n-1);                               //use derivative
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else if ((lowdiff != 0) && (supdiff == 0)) {                            //if supdiff is 0 the derivative at x_max has to be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=n*pow(a[0].sup(),n-1);                               //use derivative
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else c[2]=n*power(a[0],n-1)*a[2];                                       //else use intervalderivative =F_s
                                break;
                        case 1: if ((a[0].inf() >= 0) && (lowdiff != 0) && (supdiff != 0)) {            //if n is odd power(x,n) is convex for positiv numbers and concave for negative ones. if lowiff and supdiff are both not 0 and x_min is positive, power(x,n) is convex and the kolev formula can be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else if ((a[0].inf() >= 0) && (lowdiff == 0) && (supdiff != 0)) {       //if lowdiff is 0 the derivative at x_min has to be used
                                        double low=n*pow(a[0].inf(),n-1);                               //use derivative
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else if ((a[0].inf() >= 0) && (lowdiff != 0) && (supdiff == 0)) {       //if supdiff is 0 the derivative at x_max has to be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=n*pow(a[0].sup(),n-1);                               //use derivative
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff != 0) && (supdiff != 0)) {       //if lowdiff and supdiff are both not 0 and x_max is negative, power(x,n) is concave and the kolev formula can be used.
                                        double sup=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double low=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff == 0) && (supdiff != 0)) {       //if lowdiff is 0 the derivative at x_min has to be used
                                        double sup=n*pow(a[0].inf(),n-1);                               //use derivative
                                        double low=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff != 0) && (supdiff == 0)) {       //if supdiff is 0 the derivative at x_max has to be used
                                        double sup=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double low=n*pow(a[0].sup(),n-1);                               //use derivative
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                }
                                else c[2]=n*power(a[0],n-1)*a[2];                                       //else use intervalderivative. =F_s
                                break;
                }
        d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=power(a[0], n);                                                                    //=F_x
                c[1]=power(a[1], n);                                                                    //=F_z
                c[3]=power(a[3], n);                                                                    //=F_y
                c[6]=a[6];                                                                              //=x
                double lowdiff=a[0].inf()-a[1].inf();                                                   //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                                                   //=x_max - z
                switch (n % 2) {                                                                        //if n even power(x,n) is convex
                        case 0: if (n == 2) {                                                           //if n is 2 a direct formula can be used
                                                c[2]=a[2] * (a[0] + a[1]);                              //=F_s
                                                c[4]=2.*a[1]*a[4];                                      //=f'(z)
                                                c[5]=(a[0]+a[1])*a[5]+a[4]*a[2];                        //=F2_s
                                }
                                else if ((a[0].inf() >= 0) && (lowdiff != 0) && (supdiff != 0)) {       //if lowidff and supdiff are both not 0 and x_min is positive, power(x,n) is convex and power(x,n)'is convex and the kolev formula can be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=(low-n*pow(a[1].inf(),n-1))/lowdiff;                //use formula from kolev for second order slopes
                                        double sup2=(sup-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                //      cout<<"power"<<low<<","<<sup<<","<<low2<<","<<sup2<<","<<n*pow(a[1].inf(),n-1)<<","<<c[5].empty()<<"\n";
                                }
                                else if ((a[0].inf() >= 0) && (lowdiff == 0) && (supdiff != 0)) {       //if lowdiff is 0 the derivative at x_min has to be used
                                        double low=n*pow(a[0].inf(),n-1);                               //use derivative
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=n*(n-1)*pow(a[0].inf(),n-2);                        //use derivative
                                        double sup2=(sup-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                
                                }
                                else if ((a[0].inf() >= 0) && (lowdiff != 0) && (supdiff == 0)) {       //if supdiff is 0 the derivative at x_max has to be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=n*pow(a[0].sup(),n-1);                               //use derivative
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=(low-c[1].inf())/lowdiff;                           //use formula from kolev for second order slopes
                                        double sup2=n*(n-1)*pow(a[0].sup(),n-2);                        //use derivative
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff != 0) && (supdiff != 0)) {       //if lowidff and supdiff are both not 0 and x_max is negative, power(x,n) is convex and power(x,n)'is concave and the kolev formula can be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double sup2=(sup-n*pow(a[1].inf(),n-1))/lowdiff;                //use formula from kolev for second order slopes
                                        double low2=(low-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff == 0) && (supdiff != 0)) {       //if lowdiff is 0 the derivative at x_min has to be used
                                        double low=n*pow(a[0].inf(),n-1);                               //use derivative
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double sup2=n*(n-1)*pow(a[0].inf(),n-2);                        //use derivative
                                        double low2=(low-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff != 0) && (supdiff == 0)) {       //if supdiff is 0 the derivative at x_max has to be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=n*pow(a[0].sup(),n-1);                               //use derivative
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double sup2=(sup-n*pow(a[1].inf(),n-1))/lowdiff;                //use formula from kolev for second order slopes
                                        double low2=n*(n-1)*pow(a[0].sup(),n-2);                        //use derivative
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else {                                                                  //if power(x,n) is neither convex nor concave use intervalderivative
                                        c[2]=n*power(a[0],n-1)*a[2];                                    //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        c[5]=0.5*(n*(n-1)*power(a[0],n-2)*power(a[2],2))+n*power(a[0],n-1)*a[5];        //=F2_s
                                }
                                
                                if (c[5].empty()){                                                      //if roundingerrors are to big, it can happen that low2 > sup2, so c[5] is empty. Then the derivative is used instead of the kolev formula.
                                        c[2]=n*power(a[0],n-1)*a[2];                                    //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        c[5]=0.5*(n*(n-1)*power(a[0],n-2)*power(a[2],2))+n*power(a[0],n-1)*a[5];        //=F2_s
                        //              cout<<"else"<<c[5];
                                }
                                break;
                        case 1: if ((a[0].inf() >= 0) && (lowdiff != 0) && (supdiff != 0)) {            //if n is odd power(x,n) is convex for positiv numbers and concave for negative ones. if lowidff and supdiff are both not 0 and x_min is positive, power(x,n) is convex and power(x,n)'is convex and the kolev formula can be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=(low-n*pow(a[1].inf(),n-1))/lowdiff;                //use formula from kolev for second order slopes
                                        double sup2=(sup-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else if ((a[0].inf() >= 0) && (lowdiff == 0) && (supdiff != 0)) {       //if lowdiff is 0 the derivative at x_min has to be used
                                        double low=n*pow(a[0].inf(),n-1);                               //use derivative
                                        double sup=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=n*(n-1)*pow(a[0].inf(),n-2);                        //use derivative
                                        double sup2=(sup-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //F2_s
                                }
                                else if ((a[0].inf() >= 0) && (lowdiff != 0) && (supdiff == 0)) {       //if supdiff is 0 the derivative at x_max has to be used
                                        double low=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double sup=n*pow(a[0].sup(),n-1);                               //use derivative
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=(low-c[1].inf())/lowdiff;                           //use formula from kolev for second order slopes
                                        double sup2=n*(n-1)*pow(a[0].sup(),n-2);                        //use derivative
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff != 0) && (supdiff != 0)) {       //if n is odd power(x,n) is convex for positiv numbers and concave for negative ones. if lowidff and supdiff are both not 0 and x_max is negative, power(x,n) is concave and power(x,n)'is convex and the kolev formula can be used
                                        double sup=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double low=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=(sup-n*pow(a[1].inf(),n-1))/lowdiff;                //use formula from kolev for second order slopes
                                        double sup2=(low-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff == 0) && (supdiff != 0)) {       //if lowdiff is 0 the derivative at x_min has to be used
                                        double sup=n*pow(a[0].inf(),n-1);                               //use derivative
                                        double low=(pow(a[0].sup(),n)-pow(a[1].sup(),n))/supdiff;       //use formula from kolev for first order slopes
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=n*(n-1)*pow(a[0].inf(),n-2);                        //use derivative
                                        double sup2=(low-n*pow(a[1].sup(),n-1))/supdiff;                //use formula from kolev for second order slopes
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else if ((a[0].sup() <= 0) && (lowdiff != 0) && (supdiff == 0)) {       //if supdiff is 0 the derivative at x_max has to be used
                                        double sup=(pow(a[0].inf(),n)-pow(a[1].inf(),n))/lowdiff;       //use formula from kolev for first order slopes
                                        double low=n*pow(a[0].sup(),n-1);                               //use derivative
                                        c[2]=I(low,sup)*a[2];                                           //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        double low2=(sup-n*pow(a[1].inf(),n-1))/lowdiff;                //use formula from kolev for second order slopes
                                        double sup2=n*(n-1)*pow(a[0].sup(),n-2);                        //use derivative
                                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                                }
                                else {                                                                  //if power(x,n) is neither convex nor concave use intervalderivative
                                        c[2]=n*power(a[0],n-1)*a[2];                                    //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        c[5]=0.5*(n*(n-1)*power(a[0],n-2)*power(a[2],2))+n*power(a[0],n-1)*a[5];        //=F2_s
                                }
                                
                                if (c[5].empty()) {                                                     //if roundingerrors are to big, it can happen that low2 > sup2, so c[5] is empty. Then the derivative is used instead of the kolev formula.
                                        c[2]=n*power(a[0],n-1)*a[2];                                    //=F_s
                                        c[4]=n*power(a[1],n-1)*a[4];                                    //=f'(z)
                                        c[5]=0.5*(n*(n-1)*power(a[0],n-2)*power(a[2],2))+n*power(a[0],n-1)*a[5];        //=F2_s
                                }
                                break;
                }
                
        c.minmaxval[0]=pow(a.minmaxval[0],n);                                                           //=f(x_min)
        c.minmaxval[1]=pow(a.minmaxval[1],n);                                                           //=f(x_max)
        d=c;
        }
        return d;
}

inline Islope log(const Islope &a)
{
        if (a[0].inf()<0){
                std::cerr<<"log: interval contains negative numbers";
                abort();
        }
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=log(a[0]);                                                 //=F_x
                c[1]=log(a[1]);                                                 //=F_z
                double lowdiff=a[0].inf()-a[1].inf();                           //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                           //=x_max - z
                if ((lowdiff==0)&&(supdiff==0)){                                //if lowdiff and supdiff are 0 use intervalderivative
                        c[2]=a[2]/a[0];                                         //=F_s
                }
                if ((lowdiff!=0)&&(supdiff!=0)) {                               //if lowdiff and supdiff are both not 0, log is concave and the kolev formula can be used
                        double sup=(log(a[0].inf())-log(a[1].inf()))/lowdiff;   //use formula from kolev for first order slopes
                        double low=(log(a[0].sup())-log(a[1].sup()))/supdiff;   //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                }
                if ((lowdiff!=0)&&(supdiff==0)) {                               //if supdiff is 0 the derivative at x_max has to be used
                        double sup=(log(a[0].inf())-log(a[1].inf()))/lowdiff;   //use formula from kolev for first order slopes
                        double low=1./a[0].sup();                               //use derivative
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                }
                if ((lowdiff==0)&&(supdiff!=0)) {                               //if lowdiff is 0 the derivative at x_min has to be used
                        double sup=1./a[0].inf();                               //use derivative
                        double low=(log(a[0].sup())-log(a[1].sup()))/supdiff;   //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                }
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=log(a[0]);                                                 //=F_x
                c[1]=log(a[1]);                                                 //=F_z
                c[3]=log(a[3]);                                                 //=F_y
                c[6]=a[6];                                                      //=x
                double lowdiff=a[0].inf()-a[1].inf();                           //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                           //=x_max - z
                if ((lowdiff==0)&&(supdiff==0)){                                //if lowdiff and supdiff is 0 use intervalderivative
                        c[2]=a[2]/a[0];                                         //=F_s
                        c[4]=a[4]/a[1];                                         //=f'(z)
                        c[5]=a[5]/a[0]-0.5*power(a[2],2)/power(a[0],2);         //=F2_s
                }
                if ((lowdiff!=0)&&(supdiff!=0)){                                //if lowdiff and supdiff are both not 0, log is concave and the kolev formula can be used
                        double sup=(log(a[0].inf())-log(a[1].inf()))/lowdiff;   //use formula from kolev for first order slopes
                        double low=(log(a[0].sup())-log(a[1].sup()))/supdiff;   //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[4]=a[4]/a[1];                                         //=f'(z)
                        double sup2=(low-1./a[1].sup())/supdiff;                //use formula from kolev for second order slopes
                        double low2=(sup-1./a[1].inf())/lowdiff;                //use formula from kolev for second order slopes
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                if ((lowdiff!=0)&&(supdiff==0)){                                //if supdiff is 0 the derivative at x_max has to be used
                        double sup=(log(a[0].inf())-log(a[1].inf()))/lowdiff;   //use formula from kolev for first order slopes
                        double low=1./a[0].sup();                               //use derivative
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[4]=a[4]/a[1];                                         //=f'(z)
                        double sup2=-0.5/pow(a[0].sup(),2);                     //use derivative
                        double low2=(sup-1./a[1].inf())/lowdiff;                //use formula from kolev for second order slopes
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                if ((lowdiff==0)&&(supdiff!=0)){                                //if lowdiff is 0 the derivative at x_min has to be used
                        double sup=1./a[0].inf();                               //use derivative
                        double low=(log(a[0].sup())-log(a[1].sup()))/supdiff;   //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                   //=F_s
                        c[4]=a[4]/a[1];                                         //=f'(z)
                        double sup2=(low-1./a[1].sup())/supdiff;                //use formula from kolev for second order slopes
                        double low2=-0.5/pow(a[0].inf(),2);                     //use derivative
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];            //=F2_s
                }
                if (c[5].empty()){                                              //if roundingerrors are to big, it can happen that low2 > sup2, so c[5] is empty. Then the derivative is used instead of the kolev formula.
                        c[2]=a[2]/a[0];                                         //=F_s
                        c[5]=a[5]/a[0]-0.5*power(a[2],2)/power(a[0],2);         //=F2_s
                }
                c.minmaxval[0]=log(a.minmaxval[0]);                             //=f(x_min)
                c.minmaxval[1]=log(a.minmaxval[1]);                             //=f(x_max)
                d=c;
        }
        return d;
}

inline Islope atan(const Islope &a)
{
        Islope d;
        if (a.order==1){
                Islope c;
                c[0]=atan(a[0]);                                                        //=F_x
                c[1]=atan(a[1]);                                                        //=F_z
                double lowdiff=a[0].inf()-a[1].inf();                                   //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                                   //=x_max - z
        
                if ((a[0].inf() >= 0) && (lowdiff != 0) && (supdiff != 0)) {            //if lowdiff and supdiff are both not 0 and x_min is positive atan(x) is concave and the kolev formula can be used
                        double sup=(atan(a[0].inf())-atan(a[1].inf()))/lowdiff;         //use formula from kolev for first order slopes
                        double low=(atan(a[0].sup())-atan(a[1].sup()))/supdiff;         //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                           //=F_s
                }
                else if ((a[0].sup() <= 0) && (lowdiff != 0) && (supdiff != 0)) {       //if lowdiff and supdiff are both not 0 and x_max is negative atan(x) is convex and the kolev formula can be used
                        double low=(atan(a[0].inf())-atan(a[1].inf()))/lowdiff;         //use formula from kolev for first order slopes
                        double sup=(atan(a[0].sup())-atan(a[1].sup()))/supdiff;         //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                           //=F_s
                }
                else c[2]=a[2]/(1.+power(a[0],2));                                      //else use intervalderivative =F_s
                d=c;
        }
        if (a.order==2){
                Islope c(2);
                d.order=2;
                c[0]=atan(a[0]);                                                        //=F_x
                c[1]=atan(a[1]);                                                        //=F_z
                c[3]=atan(a[3]);                                                        //=F_y
                c[6]=a[6];                                                              //=x
                double lowdiff=a[0].inf()-a[1].inf();                                   //=x_min - z
                double supdiff=a[0].sup()-a[1].sup();                                   //=x_max - z
        
                if ((a[0].inf() >= sqrt(3.)/3) && (lowdiff != 0) && (supdiff != 0)) {   //if supdiff and lowdiff are both not 0 and x_min is bigger than sqrt3/3 then atan(x) is concave and atan'(x) is concave
                        double sup=(atan(a[0].inf())-atan(a[1].inf()))/lowdiff;         //use formula from kolev for first order slopes
                        double low=(atan(a[0].sup())-atan(a[1].sup()))/supdiff;         //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                           //=F_s
                        c[4]=a[4]/(1.+power(a[1],2));                                   //=f'(z)
                        double sup2=(low-1./(1.+power(a[1],2).sup()))/supdiff;          //use formula from kolev for second order slopes
                        double low2=(sup-1./(1.+power(a[1],2).inf()))/lowdiff;          //use formula from kolev for second order slopes
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                        std::cout<<"daaaaaaaaaa"<<c[5];
                }
                else if ((a[0].inf() >= 0) && (a[0].sup() <= sqrt(3.)/3) && (lowdiff != 0) && (supdiff != 0)) { //if supdiff and lowdiff are both not 0 and x is in [0,sqrt3/3] then atan(x) is concave and atan'(x) is convex
                        double sup=(atan(a[0].inf())-atan(a[1].inf()))/lowdiff;         //use formula from kolev for first order slopes
                        double low=(atan(a[0].sup())-atan(a[1].sup()))/supdiff;         //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                           //=F_s
                        c[4]=a[4]/(1.+power(a[1],2));                                   //=f'(z)
                        double low2=(low-1./(1.+power(a[1],2).sup()))/supdiff;          //use formula from kolev for second order slopes
                        double sup2=(sup-1./(1.+power(a[1],2).inf()))/lowdiff;          //use formula from kolev for second order slopes
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                }
                else if ((a[0].inf() >= -sqrt(3.)/3) && (a[0].sup() <= 0) && (lowdiff != 0) && (supdiff != 0)) {        //if supdiff and lowdiff are both not 0 and x is in [-sqrt3/3,0] then atan(x) is convex and atan'(x) is convex
                        double low=(atan(a[0].inf())-atan(a[1].inf()))/lowdiff;         //use formula from kolev for first order slopes
                        double sup=(atan(a[0].sup())-atan(a[1].sup()))/supdiff;         //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                           //=F_s
                        c[4]=a[4]/(1.+power(a[1],2));                                   //=f'(z)
                        double low2=(sup-1./(1.+power(a[1],2).sup()))/supdiff;          //use formula from kolev for second order slopes
                        double sup2=(low-1./(1.+power(a[1],2).inf()))/lowdiff;          //use formula from kolev for second order slopes
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                }
                else if ((a[0].sup() <= -sqrt(3.)/3) && (lowdiff != 0) && (supdiff != 0)) {     //if supdiff and lowdiff are both not 0 and x_max is smaller than -sqrt3/3 then atan(x) is convex and atan'(x) is concave
                        double low=(atan(a[0].inf())-atan(a[1].inf()))/lowdiff;         //use formula from kolev for first order slopes
                        double sup=(atan(a[0].sup())-atan(a[1].sup()))/supdiff;         //use formula from kolev for first order slopes
                        c[2]=I(low,sup)*a[2];                                           //=F_s
                        c[4]=a[4]/(1.+power(a[1],2));                                   //=f'(z)
                        double sup2=(sup-1./(1.+power(a[1],2).sup()))/supdiff;          //use formula from kolev for second order slopes
                        double low2=(low-1./(1.+power(a[1],2).inf()))/lowdiff;          //use formula from kolev for second order slopes
                        c[5]=I(low,sup)*a[5]+I(low2,sup2)*a[4]*a[2];                    //=F2_s
                }
                else {                                                                  //else use intervalderivative
                        c[2]=a[2]/(1.+power(a[0],2));                                   //=F_s
                        c[4]=a[4]/(1.+power(a[1],2));                                   //=f'(z)
                        c[5]=-a[0]*power(a[2],2)/(power(1.+power(a[0],2),2))+a[5]/(1.+power(a[0],2));   //=F2_s
                        std::cout<<"hier";
                        }
                
                if (c[5].empty()){                                                      //if roundingerrors are to big, it can happen that low2 > sup2, so c[5] is empty. Then the derivative is used instead of the kolev formula.
                        c[2]=a[2]/(1.+power(a[0],2));                                   //=F_s
                        c[5]=-a[0]*power(a[2],2)/(power(1.+power(a[0],2),2))+a[5]/(1.+power(a[0],2));   //=F2_s
                }
                c.minmaxval[0]=atan(a.minmaxval[0]);                                    //=f(x_min)
                c.minmaxval[1]=atan(a.minmaxval[1]);                                    //=f(x_max)
                d=c;
        }
                
        return d;
}

inline Islope atan2(const Islope &a, const Islope &b)
{
        Islope c;                       //atan2(a,b) = atan(x/y) so this formula is used
        c=b/a;
        c=atan(c);
        return c;
}

/*INTERVALL SLOPE*/
inline Islope::~Islope()
{
        if (data) delete data;
}

inline std::ostream &operator << (std::ostream &s, const Islope &a)
{
        if (a.order==1){
                s<<"< Range=";
                s<<a[0];
                s<<", Fixed Point=";
                s<<a[1];
                s<<", Slope=";
                s<<a[2];
                s<<" >\n";
        }
        if (a.order==2){
                s.precision(6);
                s<<std::showpoint;
                double slopemin_fmin=a[5].inf()*pow((a[6].inf()-a[6].sup()),2)/4.+(a[6].inf()-a[6].sup())*a[4].inf()/2.+a[1].inf();
                double slopemax_fmin=a[5].sup()*pow((a[6].inf()-a[6].sup()),2)/4.+(a[6].inf()-a[6].sup())*a[4].sup()/2.+a[1].sup();
                double slopemin_fmax=a[5].inf()*pow((a[6].inf()-a[6].sup()),2)/4.-(a[6].inf()-a[6].sup())*a[4].inf()/2.+a[1].inf();
                double slopemax_fmax=a[5].sup()*pow((a[6].inf()-a[6].sup()),2)/4.-(a[6].inf()-a[6].sup())*a[4].sup()/2.+a[1].sup();
                double slopemin1_fmin=-a[2].inf()*(a[6].sup()-a[6].inf())/2.+a[1].inf();
                double slopemax1_fmin=-a[2].sup()*(a[6].sup()-a[6].inf())/2.+a[1].sup();
                double slopemin1_fmax=a[2].inf()*(a[6].sup()-a[6].inf())/2.+a[1].inf();
                double slopemax1_fmax=a[2].sup()*(a[6].sup()-a[6].inf())/2.+a[1].sup();
                s<<"Range_x=";
                s<<a[0];
                s<<", Fixed Point_z=";
                s<<a[1];
                s<<"\n";
                s<<"Slope_x=";
                s<<a[2];
                s<<", Fixed Point_y=";
                s<<a[3];
                s<<"\n";
                s<<"Slope_2=";
                s<<a[5];
                s<<", Diff_z=       ";
                s<<a[4];
                s<<"\n";
                s<<"+---------------------------------------------------------+\n";
                s<<"|             slope_min      real       slope_max         |\n";    
                s<<"|               f(min)=     f(min)=     f(min)=           |\n";
                s<<"| order 1  ";
                s.width(12);
                s<<slopemin1_fmin;
                s.width(12);
                s<<a.minmaxval[0];
                s.width(12);
                s<<slopemax1_fmin;
                s<<"           |\n";
                s<<"| order 2  ";
                s.width(12);
                s<<slopemin_fmin;
                s.width(12);
                s<<a.minmaxval[0];
                s.width(12);
                s<<slopemax_fmin;
                s<<"           |\n";
                s<<"|               f(max)=     f(max)=     f(max)=           |\n";
                s<<"| order 1  ";
                s.width(12);
                s<<slopemin1_fmax;
                s.width(12);
                s<<a.minmaxval[1];
                s.width(12);
                s<<slopemax1_fmax;
                s<<"           |\n";
                s<<"| order 2  ";
                s.width(12);
                s<<slopemin_fmax;
                s.width(12);
                s<<a.minmaxval[1];
                s.width(12);
                s<<slopemax_fmax;
                s<<"           |\n";
                s<<"+---------------------------------------------------------+\n";
                if (((slopemin_fmin-a.minmaxval[0]<0) && (slopemax_fmin-a.minmaxval[0]<0)) || ((slopemin_fmin-a.minmaxval[0]>0) && (slopemax_fmin-a.minmaxval[0]>0))){
                        s<<"  !!slopefehler min!! "<<slopemin_fmin-a.minmaxval[0]<<" "<<slopemax_fmin-a.minmaxval[0]<<"\n";
                }
                if (((slopemin_fmax-a.minmaxval[1]<0) && (slopemax_fmax-a.minmaxval[1]<0)) || ((slopemin_fmax-a.minmaxval[1]>0) && (slopemax_fmax-a.minmaxval[1]>0))){
                        s<<"  !!slopefehler max!! "<<slopemin_fmax-a.minmaxval[1]<<" "<<slopemax_fmax-a.minmaxval[1]<<"\n";
                }
                s<<"\n";
        }
        return s;
}

#endif

---