Print if and only if it beats previous record.

3 files changed, 70 insertions, 27 deletions

diff --git a/charf.c b/charf.c
index c6a32c3..d6da18b 100644
--- a/charf.c
+++ b/charf.c
@@ -31,16 +31,19 @@
 /*maximum length of the support of f*/
 static const int N=24;
 
+/*maximal order of derivative*/
+static const int K=24;
+
 /*given function df[0] on domain [0,M-1], compute derivatives f' until f^{(K)} and store f^{(K)} in df*/
-void differentiate(VALUETYPE* f, VALUETYPE* df, int D, int K){
+void differentiate(VALUETYPE* f, VALUETYPE* df, int D, int k){
 	VALUETYPE df0[D];
 	/*Set zeroth derivative to be f.*/
 	for(int i=0; i<D; i++){
 		df0[i] = f[i];
 		df[i] = f[i];
 	}
-	for(int k=1; k<=K; k++){
-		/*Compute kth derivative of f from (k-1)th.*/
+	for(int l=1; l<=k; l++){
+		/*Compute lth derivative of f from (l-1)th.*/
 		for(int i=0; i<D; i++) df[i] = difference(df0[(i+1)%D], df0[i]);
 		for(int i=0; i<D; i++) df0[i] = df[i];
 	}
@@ -114,7 +117,7 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){
 	}
 }
 
-void compute(EXPTYPE p, int D, VALUETYPE* f){
+void compute(EXPTYPE p, int D, VALUETYPE* f, VALUETYPE* records){
 	VALUETYPE Mf[N];
 	/*This is the only O(D^2) operation in here so makes a lot of sense to only compute once and avoid repeating it.*/
 	compute_maximalfunction(f,Mf,D);
@@ -132,8 +135,8 @@ void compute(EXPTYPE p, int D, VALUETYPE* f){
 
 	VALUETYPE tmp[N];
 
-	for(int k=1; k<=16; k++){
-		/*Compute kth derivative of f and Mf*/
+	for(int k=1; k<=K; k++){
+		/*Compute kth derivative of f and Mf from (k-1)th derivative*/
 		differentiate(df[(k+1)%2],df[k%2],D,1);
 		differentiate(dMf[(k+1)%2],dMf[k%2],D,1);
 		
@@ -158,15 +161,17 @@ void compute(EXPTYPE p, int D, VALUETYPE* f){
 		/*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/
 
 		VALUETYPE r = ratio(intdMfp, intdfp);
+
 		double t = to_double(r);
 		if(
-		//(k==1 && t>.9)
-		//||
-		//(k==2 && t>=.5)
-		//||
+		(
+		(k==1 && t>=1)
+		||
+		(k==2 && t>=.5)
+		||
 		(k==3 && t>.53)
-		//||
-		//(k==4 && t>=.5)
+		||
+		(k==4 && t>=.5)
 		||
 		(k==5 && t>=.58)
 		||
@@ -191,13 +196,36 @@ void compute(EXPTYPE p, int D, VALUETYPE* f){
 		(k==15 && t>=.9817)
 		||
 		(k==16 && t>=.9817)
+		||
+		(k==17 && t>=.987)
+		||
+		(k==18 && t>=.991)
+		||
+		(k==19 && t>=.994)
+		||
+		(k==20 && t>=1.001)
+		||
+		(k==21 && t>=1.009)
+		||
+		(k==22 && t>=1.009)
+		||
+		(k==23 && t>=1.003)
+		||
+		(k==24 && t>=1.174)
+		)
+		&&
+		(is_greater(r,records[k]))
 		){
-			printf("f: ");
-			for(int i=0; i<D; i++) printf("%1.0f ",to_double(f[i]));
-			printf("\n");
-			char s[128];
-			root_to_string(s,r,p);
-			printf("ratio for %dth derivative: %s\n",k,s);
+			char s[256];
+			sprintf(s,"f: ");
+			int l = 3;
+			for(int i=0; i<D; i++) l += sprintf(s+l,"%1.0f ",to_double(f[i]));
+			for(int i=D; i<N; i++) l += sprintf(s+l,"  ");
+			char rts[128];
+			root_to_string(rts,r,p);
+			l += sprintf(s+l,"%dth der.: %s\n",k,rts);
+			printf("%s",s);
+			records[k] = maximum(records[k],r);
 		}
 	}
 }
@@ -224,31 +252,43 @@ int generate_function(VALUETYPE* f, int i){
 
 }
 
+typedef struct { VALUETYPE* records; int num_thread; int* domain_current; } Args;
+
 /*Goes through all functions indexed by those numbers which equal the thread number module NUM_THREADS.*/
 void* compute_chunk(void* arguments){
-	int i = *((int*)arguments);
+	Args* args = arguments;
+	int i = args -> num_thread;
+	VALUETYPE* records = args -> records;
+	int* domain_current = args -> domain_current;
+	/*exponent p of the L^p norm to consider*/
 	EXPTYPE p = int_to_exptype(1);
 
 	VALUETYPE f[N];
 	int d = generate_function(f,i);
 	while(d >= 0){
-		compute(p,d,f);
+		if(d>*domain_current){
+			*domain_current=d;
+			printf("Start considering length: %d\n",d);
+		}
+		compute(p,d,f,records);
 		i += NUM_THREADS;
 		d = generate_function(f,i);
 	}	
 }
 
 int main() {
-	/*exponent p of the L^p norm to consider*/
-	EXPTYPE p = int_to_exptype(1);
-
 	pthread_t threads[NUM_THREADS];
-	int args[NUM_THREADS];
+	Args args[NUM_THREADS];
+	VALUETYPE records[K+1];
+	int domain_current = 0;
+	for(int k=0; k<=K; k++) records[k] = int_to_valuetype(0);
 	int result_code;
 
 	for (int i = 0; i < NUM_THREADS; i++) {
 		printf("In main: Creating thread %d.\n", i);
-		args[i] = i;
+		args[i].records = records;
+		args[i].num_thread = i;
+		args[i].domain_current = &domain_current;
 		result_code = pthread_create(&threads[i], NULL, compute_chunk, args+i);
 	}
 
diff --git a/double-error.c b/double-error.c
index 7d11a92..672956f 100644
--- a/double-error.c
+++ b/double-error.c
@@ -40,7 +40,7 @@ double_error difference(double_error de1, double_error de2){
 	return sum(de1,de2m);
 }
 
-bool is_greater(double_error de1, double_error de2){ return (de1.v > de2.v); }
+bool is_greater(double_error de1, double_error de2){ return ((1+EPS)*(de1.v+de1.e) > de2.v-de2.e); }
 
 double_error maximum(double_error de1, double_error de2){
 	double_error de;
diff --git a/ratio.c b/ratio.c
index ad2b0ea..4b5c56b 100644
--- a/ratio.c
+++ b/ratio.c
@@ -84,7 +84,10 @@ rational difference(rational r1, rational r2){
 	return sum(r1,r2);
 }
 
-bool is_greater(rational r1, rational r2){ return !(difference(r1,r2).s); }
+bool is_greater(rational r1, rational r2){
+	rational diff = difference(r1,r2);
+	return !diff.s && diff.n > 0;
+}
 
 rational maximum(rational r1, rational r2){
 	if(is_greater(r1,r2)) return r1;