From a821b9953f9386614fb7657919c431d0c90eece3 Mon Sep 17 00:00:00 2001
From: Julian Weigt <juw@posteo.de>
Date: Sat, 27 Dec 2025 09:40:19 +0000
Subject: Add error estimates.

---
 charf.c | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/charf.c b/charf.c
index 1508c61..6e860ca 100644
--- a/charf.c
+++ b/charf.c
@@ -132,9 +132,15 @@ void compute(EXPTYPE p, int t){
 		//if(i%3==0) f[2*N+i+K/2] = 1;
 	}
 
+	double max = 1.0;
+	unsigned int max_operations_per_variable = 1;
+
 	VALUETYPE Mf[D];
 	compute_maximalfunction(f, Mf);
 
+	max_operations_per_variable += (D-1)*(1+max_operations_per_variable); //sum of nominators and denominators
+	max_operations_per_variable += 2*(1+max_operations_per_variable)+1; //division of nominator and denominator
+
 	/*Allocate memory for derivatives.*/
 	VALUETYPE df[D];
 	VALUETYPE dMf[D];
@@ -143,6 +149,9 @@ void compute(EXPTYPE p, int t){
 	differentiate(f,df);
 	differentiate(Mf,dMf);
 
+	max = max*(1<<K); //times 2^K
+	for(int i=1; i<= K; i++) max_operations_per_variable = max_operations_per_variable*2+1;
+
 	
 	/*Print derivatives*/
 	/*
@@ -159,11 +168,19 @@ void compute(EXPTYPE p, int t){
 	/*Compute L^p norm of derivatives*/
 	VALUETYPE intdfp = integratep(df,p);
 	VALUETYPE intdMfp = integratep(dMf,p);
+
+	max = max*pow(D,to_double(ratio(convert_int(1),convert_int(p))));
+	max_operations_per_variable = D*(max_operations_per_variable+1) + D-1 + 1;
 	//printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]);
 
 	/*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/
+
+	double divisor = to_double(intdfp);
+
 	VALUETYPE r = ratio(intdMfp, intdfp);
 
+	max_operations_per_variable = 2*max_operations_per_variable+1;
+
 	//printf("%.3d: %.3f  \n",t,r);
 	/*Print f and ||Mf^{(k)}||_p/||f^{(k)}||_p if the latter is close to 1/2.*/
 	//if(to_double(r)>.4997)
@@ -173,6 +190,7 @@ void compute(EXPTYPE p, int t){
 		for(int i=0;i<D;i++) printf("%1.0f ",to_double(f[i]));
 		printf("\n");
 		printf("%.4f\n",to_double(r));
+		if(!EXACT) printf("max error: %u * %.1f * 2^{-51} / %f = %.12f\n", max_operations_per_variable, max, divisor, max_operations_per_variable*max/divisor*pow(2.0,-51));
 	}
 }
 
-- 
cgit v1.2.3