From 828654f79fa3e38e99319250de25350109786320 Mon Sep 17 00:00:00 2001
From: Julian Weigt <juw@posteo.de>
Date: Sun, 28 Dec 2025 11:59:45 +0000
Subject: Handle error estimate for maximum properly.

---
 charf.c        |  2 +-
 double-error.c | 15 +++++++++++++++
 double-error.h |  2 ++
 double.c       |  5 +++++
 double.h       |  2 ++
 ratio.c        |  5 +++++
 ratio.h        |  2 ++
 7 files changed, 32 insertions(+), 1 deletion(-)

diff --git a/charf.c b/charf.c
index 8911ffe..6e319ec 100644
--- a/charf.c
+++ b/charf.c
@@ -106,7 +106,7 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf){
 	for(int i=0;i<D;i++){
 		Mf[i] = Af[i][i];
 		for(int j=0;j<D;j++){
-			if(is_greater(Af[i][j], Mf[i])) Mf[i] = Af[i][j];
+			Mf[i] = maximum(Af[i][j], Mf[i]);
 		}
 	}
 
diff --git a/double-error.c b/double-error.c
index aaaa405..4cd1406 100644
--- a/double-error.c
+++ b/double-error.c
@@ -37,6 +37,21 @@ double_error difference(double_error de1, double_error de2){
 
 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;
+	} else {
+	de.d = de2.d;
+	de.e = de2.e + max2*EPS;
+	}
+
+	return de;
+}
+
 double_error product(double_error de1, double_error de2){
 	double_error de;
 	de.d = de1.d * de2.d;
diff --git a/double-error.h b/double-error.h
index 413d0dc..c4dadb9 100644
--- a/double-error.h
+++ b/double-error.h
@@ -6,6 +6,8 @@ double_error to_exptype(double);
 
 bool is_greater(double_error,double_error);
 
+double_error maximum(double_error,double_error);
+
 double_error sum(double_error,double_error);
 
 double_error difference(double_error,double_error);
diff --git a/double.c b/double.c
index a587a99..cf4796c 100644
--- a/double.c
+++ b/double.c
@@ -12,6 +12,11 @@ double difference(double d1, double d2){ return d1-d2; }
 
 bool is_greater(double d1, double d2){ return (d1 > d2); }
 
+double maximum(double d1, double d2){
+	if(d1>d2) return d1;
+	else return d2;
+}
+
 double product(double d1, double d2){ return d1*d2; }
 
 double ratio(double d1, double d2){ return d1/d2; }
diff --git a/double.h b/double.h
index 56065ac..e0d0a21 100644
--- a/double.h
+++ b/double.h
@@ -4,6 +4,8 @@ double to_exptype(double);
 
 bool is_greater(double,double);
 
+double maximum(double,double);
+
 double sum(double,double);
 
 double difference(double,double);
diff --git a/ratio.c b/ratio.c
index 384f944..243426d 100644
--- a/ratio.c
+++ b/ratio.c
@@ -86,6 +86,11 @@ rational difference(rational r1, rational r2){
 
 bool is_greater(rational r1, rational r2){ return !(difference(r1,r2).s); }
 
+rational maximum(rational r1, rational r2){
+	if(is_greater(r1,r2)) return r1;
+	else return r2;
+}
+
 rational product(rational r1, rational r2){
 	rational r;
 	rational s1 = {r1.s, r1.n, r2.d};
diff --git a/ratio.h b/ratio.h
index 56701f3..5899cfa 100644
--- a/ratio.h
+++ b/ratio.h
@@ -10,6 +10,8 @@ unsigned int to_exptype(unsigned int);
 
 bool is_greater(rational,rational);
 
+rational maximum(rational,rational);
+
 rational sum(rational,rational);
 
 rational difference(rational,rational);
-- 
cgit v1.2.3