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

#define EPS 2*DBL_EPSILON 

typedef struct {double d; double e;} double_error;

double_error convert_int(int i){
	double_error de;
	de.d = (double) i;
	de.e = EPS*abs(i);
	return de;
}

double_error to_exptype(double d){
	double_error de;
	de.d = d;
	de.e = max(EPS*fabs(d),DBL_MIN);
	return de;
}

double_error sum(double_error de1, double_error de2){
	double_error de;
	de.d = de1.d + de2.d;
	de.e = de1.e + de2.e + ( fabs(de1.d) + fabs(de2.d) + de1.e + de2.e )*EPS + DBL_MIN;
	return de;
}

double_error difference(double_error de1, double_error de2){
	double_error de2m;
	de2m.d = -de2.d;
	de2m.e = de2.e;
	return sum(de1,de2m);
}

bool is_greater(double_error de1, double_error de2){ return (de1.d > de2.d); }

double_error maximum(double_error de1, double_error de2){
	double_error de;
	double max1 = de1.d+de1.e;
	double max2 = de2.d+de2.e;
	if( max1 >= max2 ){
		de.d = de1.d;
		de.e = de1.e + max1*EPS + DBL_MIN;
	} else {
		de.d = de2.d;
		de.e = de2.e + max2*EPS + DBL_MIN;
	}

	return de;
}

double_error product(double_error de1, double_error de2){
	double_error de;
	de.d = de1.d * de2.d;
	double p12 = fabs(de1.d)*de2.e;
	double p21 = fabs(de2.d)*de1.e;
	double p22 = de1.e+de2.e;
	double s = p12 + p21 + p22; 
	de.e = s + ( fabs(de.d) + s )*EPS + DBL_MIN;
	return de;
}

double_error ratio(double_error de1, double_error de2){
	double_error de2i;
	de2i.d = 1/de2.d;
	if(de2.e > de2.d/2) de2i.e = INFINITY;
	else de2i.e = 2*de2.e/(de2.d*de2.d) + fabs(de2i.d)*EPS + DBL_MIN; //2* is extra margin in case de2.e is near de2.d/2
	return product(de1,de2i);

}

double_error absolute(double_error de){
	double_error dea;
	dea.d = fabs(de.d);
	dea.e = de.e;
	return dea;
}

double_error power(double_error de, double_error p) {
	double_error dep;
	dep.d = pow(de.d,p.d);
	if(de.d > 0 && p.d > 0){
		dep.e = 2*fabs(dep.d)*( de.e/fabs(de.d) + p.e/fabs(p.d) + EPS ); //again random 2*
	}
	else {
		dep.e = de.e;
	}
	return dep;
}

double to_double(double_error de){ return de.d; }

bool to_string(char* s, double_error de){
	sprintf(s,"%f +/- %.16f",de.d,de.e);
	return true;
}