Do computations on circle, run through domains of all lengths, and add small sleep to threads to avoid overlap for small domains.

1 files changed, 64 insertions, 65 deletions

diff --git a/charf.c b/charf.c
index 670e4d0..7d12037 100644
--- a/charf.c
+++ b/charf.c
@@ -1,8 +1,11 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+//for multithreading
 #include <pthread.h>
 #include <assert.h>
+//for sleep
+#include <unistd.h>
 
 #define DOUBLEMODE 0
 #define DOUBLEERRORMODE 1
@@ -29,31 +32,23 @@
 #endif
 
 /*length of the support of f*/
-static const int N=16;
+static const int N=24;
 
-/*order of the derivative to consider. Should not be larger than (D-N)/2 because then the support of f^{(K)} reaches outside of our domain.*/
-static const int K=12;
+/*order of the derivative to consider.*/
+static const int K=3;
 
-/*length of the domain, minimum N+2*K to make sure the support of f^{(k)} belongs to the domain*/
-static const int D=N+2*K;
-
-
-/*given function df[0] on domain [0,D-1], compute derivatives f' until f^{(K)} and store them in df[1] to df[K]*/
-void differentiate(VALUETYPE* f, VALUETYPE* df){
+/*given function df[0] on domain [0,M-1], compute derivatives f' until f^{(K)} and store them in df[1] to df[K]*/
+void differentiate(VALUETYPE* f, VALUETYPE* df, int D){
+	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.*/
-		/*Only compute derivatives at arguments i < D-k, because for larger i we would need data from outside the domain of f.*/
-		for(int i=0; i<D-k; i++){
-			df[i] = difference(df[i+1], df[i]);
-		}
-		/*Make them zero for arguments where their dependence reaches outside the domain.*/
-		for(int i=D-k; i<D; i++){
-			df[i] = convert_int(0);
-		}
+		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];
 	}
 
 	/*
@@ -64,7 +59,7 @@ void differentiate(VALUETYPE* f, VALUETYPE* df){
 }
 
 /*given function f on domain [0,D-1] compute pth root of integral of |f|^p*/
-VALUETYPE integratep(VALUETYPE* f, EXPTYPE p){
+VALUETYPE integratep(VALUETYPE* f, EXPTYPE p, int D){
 	VALUETYPE integralp = convert_int(0);
 	for(int i=0;i<D;i++){
 		VALUETYPE padd = power(absolute(f[i]),p);
@@ -73,40 +68,38 @@ VALUETYPE integratep(VALUETYPE* f, EXPTYPE p){
 	return integralp;
 }
 
-void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf){
+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[D][D];
+	//VALUETYPE Sf[N][N];
 	/*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/
-	//VALUETYPE Af[D][D];
+	//VALUETYPE Af[N][N];
 	/*Apparently may become too big for stack or something.*/
-	VALUETYPE* Sf[D];
-	VALUETYPE* Af[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(N*sizeof(VALUETYPE));
+		Af[i] = malloc(N*sizeof(VALUETYPE));
 	}
 
-	for(int i=0;i<D;i++) {
+	for(int i=0; i<D; i++) {
 		Sf[i][i] = f[i];
 		Af[i][i] = f[i];
 	}
 	
-	for(int n=1;n<D;n++){
+	for(int n=1; n<D; n++){
 		/*Recursively compute all integrals and averages over intervals of increasing length*/
-		for(int i=0;i+n<D;i++){
-			Sf[i][i+n] = sum(Sf[i][i+n-1], f[i+n]);
-			Sf[i+n][i] = Sf[i][i+n];
-	
-			Af[i][i+n] = ratio(Sf[i][i+n],convert_int(n+1));
-			Af[i+n][i] = Af[i][i+n];
+		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],convert_int(n+1));
 		}
 	}
 	
 	/*Compute maximal function by picking the largest average*/
-	for(int i=0;i<D;i++){
+	for(int i=0; i<D; i++){
 		Mf[i] = Af[i][i];
-		for(int j=0;j<D;j++){
-			Mf[i] = maximum(Af[i][j], Mf[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]);
 		}
 	}
 
@@ -127,45 +120,42 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf){
 	}
 }
 
-void compute(EXPTYPE p, int t){
+void compute(EXPTYPE p, int t, int D){
 	/*allocate memory for f*/
-	VALUETYPE f[D];
+	VALUETYPE f[N];
 	/*Initiate f to be zero everywhere.*/
-	for(int i=0;i<D;i++) f[i] = convert_int(0);
-	//for(int i=0;i<N;i++) f[N+i] = rand() %2;
+	for(int i=0; i<D; i++) f[i] = convert_int(0);
+	//for(int i=0; i<N; i++) f[N+i] = rand() %2;
 	/*In the middle of the domain set f to the values encoded in bit string t*/
-	for(int i=0;i<N;i++){
-		/*Since we care about the Kth derivative, which in i depends on f on [i,i+K], shift f to the right by K so that the support of f^{(K)} stays within the domain.*/
-		f[(D-(N+2*K))/2+K+i] = convert_int((t >> i) & 1);
-		//if(i%3==0) f[2*N+i+K/2] = 1;
-	}
+	for(int i=0; i<D; i++) f[i] = convert_int((t >> i) & 1);
+	//if(i%3==0) f[2*N+i+K/2] = 1;
 
-	VALUETYPE Mf[D];
-	compute_maximalfunction(f, Mf);
+	VALUETYPE Mf[N];
+	compute_maximalfunction(f,Mf,D);
 
 	/*Allocate memory for derivatives.*/
-	VALUETYPE df[D];
-	VALUETYPE dMf[D];
+	VALUETYPE df[N];
+	VALUETYPE dMf[N];
 
 	/*Compute Kth derivative of f and Mf*/
-	differentiate(f,df);
-	differentiate(Mf,dMf);
+	differentiate(f,df,D);
+	differentiate(Mf,dMf,D);
 
 	/*Print derivatives*/
 	/*
-	for(int k=0;k<=K;k++){
+	for(int k=0; k<=K; k++){
 		printf("f  %d: ",k);
-		for(int i=0;i<D;i++) printf("%+0.1f ",df[k][i]);
+		for(int i=0; i<D; i++) printf("%+0.1f ",df[k][i]);
 		printf("\n");
 		printf("Mf %d: ",k);
-		for(int i=0;i<D;i++) printf("%+0.1f ",dMf[k][i]);
+		for(int i=0; i<D; i++) printf("%+0.1f ",dMf[k][i]);
 		printf("\n");
 	}
 	*/
 
 	/*Compute L^p norm of derivatives*/
-	VALUETYPE intdfp = integratep(df,p);
-	VALUETYPE intdMfp = integratep(dMf,p);
+	VALUETYPE intdfp = integratep(df,p,D);
+	VALUETYPE intdMfp = integratep(dMf,p,D);
 
 	//printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]);
 
@@ -176,11 +166,11 @@ void compute(EXPTYPE p, int t){
 	//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)
-	if(to_double(r)>.65)
-	//if(to_double(r)>.30)
+	//if(to_double(r)>=.58)
+	if(to_double(r)>.53)
 	{
 		printf("f: ");
-		for(int i=0;i<D;i++) printf("%1.0f ",to_double(f[i]));
+		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);
@@ -188,17 +178,24 @@ void compute(EXPTYPE p, int t){
 	}
 }
 
-struct Args {EXPTYPE* p; int* t; };
+struct Args {EXPTYPE* p; int* D; int* t; };
 
 void* perform_work(void* arguments){
 	struct Args *args = arguments;
 	EXPTYPE p = *(args->p);
+	int* D = args->D;
 	int* s = args->t;
 
-	while(*s <= (1 << N)-1){
-		int t = *s;
-		*s = t+1;
-		compute(p,t);
+	while(*D <= N){
+		if(*s <= (1 << *D)-2){
+			int t = *s;
+			*s = t+1;
+			compute(p,t,*D);
+		}
+		else {
+			*s = 0;
+			*D = *D+1;
+		}
 	}
 	return NULL;
 }
@@ -209,12 +206,13 @@ int main() {
 
 	/*Iterate over all strings of 0s and 1s with length N. Those will represent f.*/
 	int t = 1;
+	int D = 1;
 
 	//for(int t=1; t<=(1 << N)-1; t++){
 		//compute(N,K,D,p,t);
 	//}
 
-	struct Args args = {&p, &t};
+	struct Args args = {&p, &D, &t};
 
 	pthread_t threads[NUM_THREADS];
 	int result_code;
@@ -223,6 +221,7 @@ int main() {
 		printf("In main: Creating thread %d.\n", i);
 		result_code = pthread_create(&threads[i], NULL, perform_work, &args);
 		assert(!result_code);
+		usleep(1000*100);
 	}
 
 	// wait for each thread to complete