From b54caf5830bc5d67d403dadd4903c5a54ec91eb0 Mon Sep 17 00:00:00 2001
From: Julian Weigt <juw@posteo.de>
Date: Mon, 9 Feb 2026 15:12:35 +0100
Subject: Add more comments to code.

---
 charf.c | 96 ++++++++++++++++++++++++++++++++++++++++++++++++-----------------
 1 file changed, 71 insertions(+), 25 deletions(-)

diff --git a/charf.c b/charf.c
index 1948971..6908832 100644
--- a/charf.c
+++ b/charf.c
@@ -12,12 +12,16 @@
 #define DOUBLEERRORMODE 1
 #define RATIOMODE 2
 
+/*Determines if computations are done with floating point type, floating point with error bounds or fractions.
+Floating point without error bounds are the default.*/
 #ifndef MODE
 #define MODE DOUBLEMODE
 #endif
 
+/*number of threads for parallel execution.*/
 #define NUM_THREADS 7
 
+/*Make sure the appropriate files are included for the type of computation and sets the data types accordingly.*/
 #if MODE == DOUBLEERRORMODE
 #include "double-error.h"
 #define VALUETYPE double_error
@@ -36,7 +40,7 @@
 
 typedef unsigned long index_t;
 
-/*maximum length of the support of f*/
+/*maximum length of the circle*/
 #define N 36
 
 /*maximal order of derivative*/
@@ -45,7 +49,7 @@ typedef unsigned long index_t;
 /*maximal number of exponents*/
 #define P 5
 
-/*given function df[0] on domain [0,M-1], compute derivatives f' until f^{(K)} and store f^{(K)} in df*/
+/*Given function f on domain [0,D-1], compute derivative of order k and store result in df.*/
 void differentiate(VALUETYPE* f, VALUETYPE* df, int D, int k){
 	VALUETYPE df0[D];
 	/*Set zeroth derivative to be f.*/
@@ -66,7 +70,8 @@ void differentiate(VALUETYPE* f, VALUETYPE* df, int D, int k){
 	*/
 }
 
-/*given function f on domain [0,D-1] compute pth root of integral of |f|^p*/
+/*Given function f on domain [0,D-1] compute pth root of integral of |f|^p.*/
+/*For p=∞ instead compute the maximum, i.e. ||f||_∞.*/
 VALUETYPE integratep(VALUETYPE* f, EXPTYPE p, int D){
 	VALUETYPE integralp = int_to_valuetype(0);
 	for(int i=0; i<D; i++){
@@ -77,48 +82,64 @@ VALUETYPE integratep(VALUETYPE* f, EXPTYPE p, int D){
 }
 
 void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
-	/*Sf[i][j] will be the integral of f on [min(i,j),max(i,j)]*/
-	//VALUETYPE Sf[N][N];
-	/*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/
-	//VALUETYPE Af[N][N];
-	/*Apparently may become too big for stack or something so have to use malloc instead.*/
+	/*Sf[i][j] will be the integral of f on [i,i+j).*/
 	VALUETYPE* Sf[N];
+	/*Af[i][j] will be the average of f on [i,i+j).*/
 	VALUETYPE* Af[N];
+	/*Would use C arrays, but apparently they may become too big for stack or something on my machine so have to use malloc instead.*/
+	//VALUETYPE Sf[N][2*N+1];
+	//VALUETYPE Af[N][2*N+1];
+
 	for(int i=0; i<D; i++){
 		Sf[i] = malloc((2*D+1)*sizeof(VALUETYPE));
 		Af[i] = malloc((2*D+1)*sizeof(VALUETYPE));
 	}
 
+	/*Recursively compute all integrals and averages over intervals of increasing length*/
+	/*We have to consider intervals of length up to twice the length of the circle, i.e. 2*D-1.*/
 	for(int i=0; i<D; i++) {
 		Sf[i][1] = f[i];
 		Af[i][1] = f[i];
 	}
-
 	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][l] = sum(Sf[i][l-1], f[(i+l-1)%D]);
 			Af[i][l] = ratio(Sf[i][l],int_to_valuetype(l));
 		}
 	}
-
+	/*Because it simplifies code later on we also assign the average on twice the circle, starting in 0.*/
 	Af[0][2*D] = Af[0][D];
 
+	/*Finding the maximal average actually looks like the most costly computation, in the whole algorithm, being of order D^3.
+	Hence we put effort into making it efficient.
+	The strategy is to go trough all possible intervals, starting with those of largest length.
+	For each [i,i+l) that we encounter we check for each n∈[i,i+l) if the average on [i,i+l) beats the current best average that we have so far computed for intervals containing n.
+	The starting point of the current best interval is stored in Pbest[n] and the length in Lbest[n].
+	*/
 	int Pbest[N];
 	int Lbest[N];
+	/*Start with the largest possible interval, [0,2*D), which contains each point n∈[0,D).*/
 	for(int i=0; i<D; i++){
 		Pbest[i]=0;
 		Lbest[i]=2*D;
 	}
-
+	/*Go through all lenghts l <= 2*D-1 in decreasing order.*/
 	for(int l=2*D-1; l>0; l--) {
+		/*Go through all intervals [i,i+l).*/
 		for(int i=0; i<D; i++){
 			VALUETYPE A = Af[i][l];
+			/*Go through all points n ∈ [i,i+l).
+			If l>D it suffices to go until i+D-1.*/
+			int intervalend = fmin(i+l,i+D);
 			int n=i;
-			while(n<i+l && n<i+D){
+			while(n < intervalend){
+				/*If the average of the current best interval for n is surely larger than A[i][l] then  the same will be true for all subsequent points in that interval.
+				Thus, to save time just skip ahead to the end of that interval.*/
 				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{
+					/*Since the new average might win set the new best average to be the maximum of the last maximum and the new.*/
 					A = maximum(A,Af[Pbest[n%D]][Lbest[n%D]]);
+					/*To avoid duplicating computations, look ahead how many subsequent points have as their latest best interval the same interval, and accordingly set their new current best average to be the same.*/
 					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;
@@ -130,11 +151,13 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
 					n=k;
 				}
 			}
+			/*Since we may be working with errors, the above action of taking maxima may have increased the error.
+			That means we may have to update the average of the current interval to carry that increased error.*/
 			Af[i][l] = A;
 		}
 	}
 	
-	/*Compute maximal function by picking the largest average*/
+	/*Finally store best averages in maximal function.*/
 	for(int i=0; i<D; i++) Mf[i] = Af[Pbest[i]][Lbest[i]];
 
 	/*Print computed functions and averages*/
@@ -154,45 +177,68 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
 	}
 }
 
-/*Generates all characteristic functions of length up N.*/
+/*Generates all characteristic functions of length up N, indexed by i, but skip trivial ones, and pick only one representative from each class of translaties.
+The return value is the length of the domain, with special values 0 which indicates we skip this function, and -1 which indicates that the index i is too large to represent a function.*/
 int generate_each_charf(VALUETYPE* f, index_t i){
+	/*In order to find which function corresponds to index i, we efficiently go through all indeces up to i starting at 0.*/
 	index_t s=0;
+	/*The smallest domain we consider is the circle of length 2.*/
 	int d=2;
-	/*number of strings of length d that begin with 1 and end with 0*/
+	/*We skip the functions that are all 0 or all 1s.
+	That means every function has a point n such that f[n]=1 and f[n+1]=0.
+	Since we also don't want to consider different translates, we may translate each function such that this value n equals d-1.
+	That means for each d we have to consider only arrays of length d that begin with 1 and end with 0, of which there exist 2^{d-2} many.*/
 	index_t powd = 1UL << d-2;
+	/*At this point s is the first index that corresponds to functions on the circle of length d.
+	Check if our index i goes past the last index that corresponds to a function on the circle of length d.*/
 	while(i-s >= powd){
+		/*Increment s,d,powd to correspond to functions of length d+1.*/
 		s += powd;
 		d++;
 		powd = powd << 1;
 	}
 
+	/*If our index i corresponds to a function on a circle of length larger than N then return -1 which will causo the program to stop.*/
 	if(d>N) return -1;
 
-	/*t starts with 1, then comes i-s and then 0.*/
+	/*t encodes a characteristic function on the circle of length d.
+	Since we only consider functions with f[0]=0 and f[d-1]=1, set its places accordingly.*/
 	index_t t = ( (i-s) << 1 ) + 1UL;
 
-	/*Indicates if we want to consider this t.*/
+	/*This variable indicates if we want to consider this t.
+	Next we will perform test to determine if it should be set to false.*/
 	bool is_representative = true;
 
-	/*Check if t encodes a function that is made up of translations of a shorter function.*/
+	/*Check if t encodes a function that is made up of multiple copies of a shorter function.
+	We skip those because they and their maximal function are periodic and thus should be considered on a smaller circle instead.*/
+	/*Go through all possible lenghts r of that shorter function.*/
 	for(int r = 1; r < d; r++){
+		/*r has to be a divisor of d.*/
 		if(d % r == 0){
 			bool is_r_copy = true;
+			/*bit mask of length r*/
 			index_t ones = 0;
 			for(int o = 0; o<r; o++) ones += 1UL <<o;
+			/*extract from t an array of length r*/
 			index_t tail = t & ones;
+			/*check if t is d/r many copies of that array*/
 			for(int n = 1; n < d/r; n++) if( ( (t>>n*r) & ones ) != tail ) is_r_copy = false;
 			if(is_r_copy) is_representative = false;
 		}
 	}
 
-	/*Check if t encodes the strictly largest of its equivalent translates.*/
-	/*Such t which may equal its translates have already been filtered out by the previous step because they are copies of at least two shorter functions.*/
-	index_t ones = 0;
-	for(int o = 0; o<d; o++) ones += 1UL << o;
-	if(is_representative) for(int n=1; n<d; n++) if( t >= ( (t<<n) & ones ) + ((t>>(d-n)) & ones) ) is_representative = false;
+	/*We aim to filter out equivalent translates.
+	The representative we pick is the one which represents the smallest integer.
+	Every function which is a nontrivial translate of itself must in fact consist of multiplie copies of a shorter array (right?).
+	Those we already filtered out in the previous step, so every nontrivial translate of that representative indeed represents a strictly smaller integer.
+	*/
+	if(is_representative){
+		index_t ones = 0;
+		for(int o = 0; o<d; o++) ones += 1UL << o;
+		for(int n=1; n<d; n++) if( t >= ( (t<<n) & ones ) + ((t>>(d-n)) & ones) ) is_representative = false;
+	}
 
-	/*Set f to the values encoded in bit string t which is a value between 1 and powd = (1<<d)-2.*/
+	/*Set f to the values encoded in bit array t.*/
 	if(is_representative){
 		for(int n=0; n<d; n++) f[n] = int_to_valuetype((t >> n) & 1UL);
 		return d;
-- 
cgit v1.2.3