From 944f90c47ffcde862dfe5f258de0b1ebf229c20e Mon Sep 17 00:00:00 2001
From: Julian Weigt <juw@posteo.de>
Date: Sun, 28 Dec 2025 11:20:32 +0000
Subject: Implement double arithmetic with error bounds.

---
 double-error.c | 83 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 83 insertions(+)
 create mode 100644 double-error.c

(limited to 'double-error.c')

diff --git a/double-error.c b/double-error.c
new file mode 100644
index 0000000..aaaa405
--- /dev/null
+++ b/double-error.c
@@ -0,0 +1,83 @@
+#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 = EPS*fabs(d);
+	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;
+	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 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;
+	return de;
+}
+
+double_error ratio(double_error de1, double_error de2){
+	double_error de2i;
+	de2i.d = 1/de2.d;
+	de2i.e = 2*de2.e/(de2.d*de2.d) + fabs(de2i.d)*EPS; //2* is random extra margin in case error is large
+	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 + EPS;
+	}
+	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;
+}
-- 
cgit v1.2.3