From a57aacd8116a2a0ba3f56f41883c9c5f30130927 Mon Sep 17 00:00:00 2001
From: Julian Weigt <juw@posteo.de>
Date: Mon, 19 Jan 2026 15:08:17 +0100
Subject: Fix computation of maximal function.

---
 charf.c | 54 ++++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 40 insertions(+), 14 deletions(-)

(limited to 'charf.c')

diff --git a/charf.c b/charf.c
index 0d11f22..33220ea 100644
--- a/charf.c
+++ b/charf.c
@@ -83,31 +83,57 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
 	VALUETYPE* Sf[N];
 	VALUETYPE* Af[N];
 	for(int i=0; i<D; i++){
-		Sf[i] = malloc(D*sizeof(VALUETYPE));
-		Af[i] = malloc(D*sizeof(VALUETYPE));
+		Sf[i] = malloc((2*D+1)*sizeof(VALUETYPE));
+		Af[i] = malloc((2*D+1)*sizeof(VALUETYPE));
 	}
 
 	for(int i=0; i<D; i++) {
-		Sf[i][i] = f[i];
-		Af[i][i] = f[i];
+		Sf[i][1] = f[i];
+		Af[i][1] = f[i];
 	}
-	
-	for(int n=1; n<D; n++){
+
+	for(int l=2; l < 2*D; l++){
 		/*Recursively compute all integrals and averages over intervals of increasing length*/
 		for(int i=0; i<D; i++){
-			Sf[i][(i+n)%D] = sum(Sf[i][(i+n-1)%D], f[(i+n)%D]);
-			Af[i][(i+n)%D] = ratio(Sf[i][(i+n)%D],int_to_valuetype(n+1));
+			Sf[i][l] = sum(Sf[i][l-1], f[(i+l-1)%D]);
+			Af[i][l] = ratio(Sf[i][l],int_to_valuetype(l));
 		}
 	}
-	
-	/*Compute maximal function by picking the largest average*/
+
+	Af[0][2*D] = Af[0][D];
+
+	int Pbest[N];
+	int Lbest[N];
 	for(int i=0; i<D; i++){
-		Mf[i] = Af[i][i];
-		for(int j=1; j<=D; j++){
-			Mf[i] = maximum(Af[i][(i+j)%D], Mf[i]);
-			Mf[i] = maximum(Af[(i-j+D)%D][i], Mf[i]);
+		Pbest[i]=0;
+		Lbest[i]=2*D;
+	}
+
+	for(int l=2*D-1; l>0; l--) {
+		for(int i=0; i<D; i++){
+			VALUETYPE A = Af[i][l];
+			int n=i;
+			while(n<i+l && n<i+D){
+				if(is_greater_certainly(Af[Pbest[n%D]][Lbest[n%D]],Af[i][l])) n += (Pbest[n%D]-n-D)%D + Lbest[n%D];
+				else{
+					A = maximum(A,Af[Pbest[n%D]][Lbest[n%D]]);
+					int k = n+1;
+					while((Pbest[k%D] == Pbest[n%D]) && (Lbest[k%D] == Lbest[n%D]) && k<i+l && k<i+D){
+						Pbest[k%D] = i;
+						Lbest[k%D] = l;
+						k++;
+					}
+					Pbest[n%D] = i;
+					Lbest[n%D] = l;
+					n=k;
+				}
+			}
+			Af[i][l] = A;
 		}
 	}
+	
+	/*Compute maximal function by picking the largest average*/
+	for(int i=0; i<D; i++) Mf[i] = Af[Pbest[i]][Lbest[i]];
 
 	/*Print computed functions and averages*/
 	/*
-- 
cgit v1.2.3