Change double to long double because numbers get too big for high derivatives.

5 files changed, 79 insertions, 60 deletions

diff --git a/charf.c b/charf.c
index d2a42c8..0028dc5 100644
--- a/charf.c
+++ b/charf.c
@@ -34,8 +34,8 @@ Floating point without error bounds are the default.
 #define EXPTYPE unsigned int
 #else
 #include "double.h"
-#define VALUETYPE double
-#define EXPTYPE double
+#define VALUETYPE long double
+#define EXPTYPE long double
 #endif
 
 #define STRING_SIZE 65536
diff --git a/double-error.c b/double-error.c
index dfbb34e..4e52b18 100644
--- a/double-error.c
+++ b/double-error.c
@@ -4,13 +4,18 @@
 #include <float.h>
 #include <stdbool.h>
 
-typedef double vtype;
-typedef double etype;
+typedef long double vtype;
+typedef long double etype;
 #define EPS 2*DBL_EPSILON 
 #define MIN 2*DBL_MIN
 
 typedef struct {vtype v; etype e;} double_error;
 
+vtype vabsolute(vtype v){
+	//return fabs(v);
+	return fabsl(v);
+}
+
 double_error int_to_valuetype(int i){
 	double_error de;
 	de.v = (vtype) i;
@@ -33,7 +38,7 @@ 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;
+	de.e = de1.e + de2.e + ( vabsolute(de1.v) + vabsolute(de2.v) + de1.e + de2.e )*EPS + MIN;
 	return de;
 }
 
@@ -45,11 +50,11 @@ double_error difference(double_error de1, double_error de2){
 }
 
 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;
+	return de1.v-de1.e > de2.v+de2.e + EPS*(vabsolute(de1.v)+vabsolute(de1.e)+vabsolute(de2.v)+vabsolute(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;
+	return de1.v+de1.e + EPS*(vabsolute(de1.v)+vabsolute(de1.e)+vabsolute(de2.v)+vabsolute(de2.e)) + MIN > de2.v-de2.e;
 }
 
 double_error maximum(double_error de1, double_error de2){
@@ -68,26 +73,26 @@ double_error maximum(double_error de1, double_error de2){
 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 p12 = vabsolute(de1.v)*de2.e*(1+EPS) + MIN;
+	vtype p21 = vabsolute(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;
+	de.e = s + ( vabsolute(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
+	if(de2.e > vabsolute(de2.v)/2) de2i.e = INFINITY;
+	else de2i.e = 2*de2.e/(de2.v*de2.v) + vabsolute(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.v = vabsolute(de.v);
 	dea.e = de.e;
 	return dea;
 }
@@ -134,9 +139,11 @@ int valuetype_to_string(char* s, double_error de){
 	if(de.e == 0.0){
 		/*If number really is an int, print accordingly to give cleaner output for example when printing the function itself, which often will be integer valued without errors.*/
 		if((vtype) ((int) de.v) == de.v) sprintf(s,"%d", (int)de.v);
-		else sprintf(s,"%4.3f…",de.v);
+		//else sprintf(s,"%4.3f…",de.v);
+		else sprintf(s,"%4.3Lf…",de.v);
 	}
-	else sprintf(s,"%4.3f… +/- %6.1e",de.v,de.e);
+	//else sprintf(s,"%4.3f… +/- %6.1e",de.v,de.e);
+	else sprintf(s,"%4.3Lf… +/- %6.1e",de.v,de.e);
 	return 0;
 }
 
@@ -149,9 +156,11 @@ int valuetype_to_latex(char* s, double_error de){
 	if(de.e == 0.0){
 		/*If number really is an int, print accordingly to give cleaner output for example when printing the function itself, which often will be integer valued without errors.*/
 		if((vtype) ((int) de.v) == de.v) sprintf(s,"%d", (int)de.v);
-		else sprintf(s,"%4.3f\\ldots",de.v);
+		//else sprintf(s,"%4.3f\\ldots",de.v);
+		else sprintf(s,"%4.3Lf\\ldots",de.v);
 	}
-	else sprintf(s,"%4.3f\\ldots\\pm\\texttt{%1.0e}",de.v,de.e);
+	//else sprintf(s,"%4.3f\\ldots\\pm\\texttt{%1.0e}",de.v,de.e);
+	else sprintf(s,"%4.3Lf\\ldots\\pm\\texttt{%1.0e}",de.v,de.e);
 	return 0;
 }
 
diff --git a/double-error.h b/double-error.h
index 3f87ab2..76ebc31 100644
--- a/double-error.h
+++ b/double-error.h
@@ -3,8 +3,8 @@
 
 #include <stdbool.h>
 
-typedef double vtype;
-typedef double etype;
+typedef long double vtype;
+typedef long double etype;
 typedef struct {vtype v; etype e;} double_error;
 
 double_error int_to_valuetype(int);
diff --git a/double.c b/double.c
index 0799e7d..5ab6178 100644
--- a/double.c
+++ b/double.c
@@ -4,46 +4,53 @@
 #include <stdio.h>
 #include <stdbool.h>
 
-double int_to_valuetype(int i){ return (double) i; }
+typedef long double vtype;
+typedef long double etype;
 
-double int_to_exptype(double d){ return d; }
+vtype int_to_valuetype(int i){ return (vtype)i; }
 
-double infinity_to_exptype(){ return INFINITY; }
+etype int_to_exptype(int d){ return (etype)d; }
 
-bool exptype_is_infinite(double d){ return d==INFINITY; }
+etype infinity_to_exptype(){ return INFINITY; }
 
-double sum(double d1, double d2){ return d1+d2; }
+bool exptype_is_infinite(etype d){ return d==INFINITY; }
 
-double difference(double d1, double d2){ return d1-d2; }
+vtype sum(vtype d1, vtype d2){ return d1+d2; }
 
-bool is_greater_certainly(double d1, double d2){ return (d1 > d2); }
+vtype difference(vtype d1, vtype d2){ return d1-d2; }
 
-bool is_greater_possibly(double d1, double d2){ return is_greater_certainly(d1,d2); }
+bool is_greater_certainly(vtype d1, vtype d2){ return (d1 > d2); }
 
-double maximum(double d1, double d2){
+bool is_greater_possibly(vtype d1, vtype d2){ return is_greater_certainly(d1,d2); }
+
+vtype maximum(vtype d1, vtype d2){
 	if(d1>d2) return d1;
 	else return d2;
 }
 
-double product(double d1, double d2){ return d1*d2; }
+vtype product(vtype d1, vtype d2){ return d1*d2; }
 
-double ratio(double d1, double d2){ return d1/d2; }
+vtype ratio(vtype d1, vtype d2){ return d1/d2; }
 
-double absolute(double d){ return fabs(d); }
+vtype absolute(vtype d){
+	//return fabs(d);
+	return fabsl(d);
+}
 
-double power(double d, double p) { return pow(d,p); }
+vtype power(vtype d, vtype p) { return pow(d,p); }
 
-double valuetype_to_double(double d){ return d; }
+double valuetype_to_double(vtype d){ return (double)d; }
 
-int valuetype_to_string(char* s, double d){
+int valuetype_to_string(char* s, vtype d){
 	if( d == 0.0 || d == 1.0 || d == 2.0 || d == 4.0 || d == 8.0 || d == 16.0 ) sprintf(s,"%d", (int)d);
-	else sprintf(s,"%4.3f",d);
+	//else sprintf(s,"%4.3f",d);
+	else sprintf(s,"%4.3Lf",d);
 	return 0;
 }
 
-int valuetype_to_latex(char* s, double d){ return valuetype_to_string(s,d); }
+int valuetype_to_latex(char* s, vtype d){ return valuetype_to_string(s,d); }
 
-int exptype_to_string(char *s, double p){ 
+int exptype_to_string(char *s, etype p){ 
 	if(exptype_is_infinite(p)){
 		sprintf(s,"inf");
 		return 0;
@@ -51,7 +58,7 @@ int exptype_to_string(char *s, double p){
 	else return valuetype_to_string(s,p);
 }
 
-int exptype_to_latex(char *s, double p){
+int exptype_to_latex(char *s, etype p){
 	if(exptype_is_infinite(p)){
 		sprintf(s,"\\infty");
 		return 0;
@@ -59,13 +66,13 @@ int exptype_to_latex(char *s, double p){
 	else return valuetype_to_latex(s,p);
 }
 
-int root_to_string(char* s, double d, double p){
+int root_to_string(char* s, vtype d, etype p){
 	if(exptype_is_infinite(p)) valuetype_to_string(s,d);
 	else valuetype_to_string(s,pow(d,1.0/p));
 	return 0;
 }
 
-int root_to_latex(char* s, double d, double p){
+int root_to_latex(char* s, vtype d, etype p){
 	if(exptype_is_infinite(p)) valuetype_to_latex(s,d);
 	else valuetype_to_latex(s,pow(d,1.0/p));
 	return 0;
diff --git a/double.h b/double.h
index 032d994..6678ba6 100644
--- a/double.h
+++ b/double.h
@@ -3,44 +3,47 @@
 
 #include <stdbool.h>
 
-double int_to_valuetype(int);
+typedef long double vtype;
+typedef long double etype;
 
-double int_to_exptype(double);
+vtype int_to_valuetype(int);
 
-double infinity_to_exptype();
+etype int_to_exptype(int);
 
-bool exptype_is_infinite(double);
+vtype infinity_to_exptype();
 
-bool is_greater_certainly(double,double);
+bool exptype_is_infinite(vtype);
 
-bool is_greater_possibly(double,double);
+bool is_greater_certainly(vtype,vtype);
 
-double maximum(double,double);
+bool is_greater_possibly(vtype,vtype);
 
-double sum(double,double);
+vtype maximum(vtype,vtype);
 
-double difference(double,double);
+vtype sum(vtype,vtype);
 
-double product(double,double);
+vtype difference(vtype,vtype);
 
-double ratio(double,double);
+vtype product(vtype,vtype);
 
-double absolute(double);
+vtype ratio(vtype,vtype);
 
-double power(double,double);
+vtype absolute(vtype);
 
-double valuetype_to_double(double);
+vtype power(vtype,vtype);
 
-int exptype_to_string(char*,double);
+double valuetype_to_double(vtype);
 
-int exptype_to_latex(char*,double);
+int exptype_to_string(char*,vtype);
 
-int valuetype_to_string(char*,double);
+int exptype_to_latex(char*,vtype);
 
-int valuetype_to_latex(char*,double);
+int valuetype_to_string(char*,vtype);
 
-int root_to_string(char*,double,double);
+int valuetype_to_latex(char*,vtype);
 
-int root_to_latex(char*,double,double);
+int root_to_string(char*,vtype,etype);
+
+int root_to_latex(char*,vtype,etype);
 
 #endif