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

typedef double vtype;
typedef double etype;
#define EPS 2*DBL_EPSILON 
#define MIN 2*DBL_MIN

typedef struct {vtype v; etype e;} double_error;

double_error int_to_valuetype(int i){
	double_error de;
	de.v = (vtype) i;
	if(abs(i) < 9999999) de.e = 0.0;
	else de.e = EPS*abs(i);
	return de;
}

double_error int_to_exptype(int i){ return int_to_valuetype(i); }

double_error infinity_to_exptype(){
	double_error de;
	de.v = INFINITY;
	de.e = 0;
	return de;
}

bool exptype_is_infinite(double_error d){ return d.v == INFINITY; }

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

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

bool is_greater_certainly(double_error de1, double_error de2){
	return de1.v-de1.e > de2.v+de2.e + EPS*(fabs(de1.v)+fabs(de1.e)+fabs(de2.v)+fabs(de2.e)) + MIN;
}

bool is_greater_possibly(double_error de1, double_error de2){
	return de1.v+de1.e + EPS*(fabs(de1.v)+fabs(de1.e)+fabs(de2.v)+fabs(de2.e)) + MIN > de2.v-de2.e;
}

double_error maximum(double_error de1, double_error de2){
	if(is_greater_certainly(de1,de2)) return de1;
	else if(is_greater_certainly(de2,de1)) return de2;
	else{
		double_error de;
		if( de1.v >= de2.v ) de.v = de1.v;
		else de.v = de2.v;
		if( de1.e >= de2.e ) de.e = de1.e;
		else de.e = de2.e;
		return de;
	}
}

double_error product(double_error de1, double_error de2){
	double_error de;
	de.v = de1.v * de2.v;
	vtype p12 = fabs(de1.v)*de2.e*(1+EPS) + MIN;
	vtype p21 = fabs(de2.v)*de1.e*(1+EPS) + MIN;
	vtype p22 = de1.e*de2.e*(1+EPS) + MIN;
	etype s = (p12 + p21 + p22)*(1+2*EPS);
	de.e = s + ( fabs(de.v) + s )*EPS + MIN;
	return de;
}

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

}

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

double_error power(double_error de, double_error p) {
	double_error dep;
	dep.v = pow(de.v,p.v);
	vtype dmin = (de.v - de.e)*(1-EPS);
	vtype dmax = (de.v + de.e)*(1+EPS);
	vtype pmin = (p.v - p.e)*(1-EPS);
	vtype pmax = (p.v + p.e)*(1+EPS);
	vtype rmin,rmax;
	if(dmin >= 1){
		rmin = pow(dmin,pmin);
		rmax = pow(dmax,pmax);
	} else if(dmax >= 1) {
		rmin = pow(dmin,pmax);
		rmax = pow(dmax,pmax);
	} else {
		rmin = pow(dmin,pmax);
		rmax = pow(dmax,pmin);
	}
	if(dmin >= 0) dep.e = rmax-rmin + 2*EPS*rmax;
	else dep.e = (1+2*EPS)*rmax;
	return dep;
}

double valuetype_to_double(double_error de){ return de.v; }

double exptype_to_double(double_error de){ return de.v; }

int valuetype_to_string(char* s, double_error de){
	if(de.e == 0.0){
		if(de.v == 0.0 || de.v == 1.0 || de.v == 2.0 || de.v == 4.0 || de.v == 8.0 || de.v == 16.0) sprintf(s,"%d", (int)de.v);
		else sprintf(s,"%4.3f…",de.v);
	}
	else sprintf(s,"%4.3f… +/- %6.1e",de.v,de.e);
	return 0;
}

int exptype_to_string(char* s, double_error de){
	if(exptype_is_infinite(de)) return sprintf(s,"inf");
	else return valuetype_to_string(s,de);
}

int valuetype_to_latex(char* s, double_error de){
	if(de.e == 0.0){
		if(de.v == 0.0 || de.v == 1.0 || de.v == 2.0 || de.v == 4.0 || de.v == 8.0 || de.v == 16.0) sprintf(s,"%d", (int)de.v);
		else sprintf(s,"%4.3f\\ldots",de.v);
	}
	else sprintf(s,"%4.3f\\ldots\\pm\\texttt{%1.0e}",de.v,de.e);
	return 0;
}

int exptype_to_latex(char* s, double_error de){ 
	if(exptype_is_infinite(de)) return sprintf(s,"\\infty");
	else return valuetype_to_latex(s,de);
}

int root_to_string(char* s, double_error de, double_error p){
	double_error de1p;
	if(exptype_is_infinite(p)) de1p = de;
	else de1p = power(de,ratio(int_to_valuetype(1),p));
	valuetype_to_string(s,de1p);
	return 0;
}

int root_to_latex(char* s, double_error de, double_error p){
	double_error de1p;
	if(exptype_is_infinite(p)) de1p = de;
	else de1p = power(de,ratio(int_to_valuetype(1),p));
	valuetype_to_latex(s,de1p);
	return 0;
}