diff options
| -rw-r--r-- | charf.c | 92 |
1 files changed, 47 insertions, 45 deletions
@@ -8,7 +8,7 @@ #define EXACT false #endif -#define NUM_THREADS 4 +#define NUM_THREADS 6 #if EXACT #include "ratio.h" @@ -20,9 +20,18 @@ #define EXPTYPE double #endif +/*length of the support of f*/ +static const int N=16; + +/*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; + +/*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, int D, int K){ +void differentiate(VALUETYPE* f, VALUETYPE* df){ /*Set zeroth derivative to be f.*/ for(int i=0; i<D; i++){ df[i] = f[i]; @@ -38,23 +47,33 @@ void differentiate(VALUETYPE* f, VALUETYPE* df, int D, int K){ df[i] = convert_int(0); } } + + /* + printf("df "); + for(int i=0;i<D;i++) printf("%+6.1f ",to_double(df[i])); + printf("\n"); + */ } /*given function f on domain [0,D-1] compute pth root of integral of |f|^p*/ -VALUETYPE integratep(VALUETYPE* f, EXPTYPE p, int D){ +VALUETYPE integratep(VALUETYPE* f, EXPTYPE p){ VALUETYPE integralp = convert_int(0); for(int i=0;i<D;i++){ - integralp = sum(integralp,power(absolute(f[i]),p)); + VALUETYPE padd = power(absolute(f[i]),p); + integralp = sum(integralp,padd); } return integralp; } -void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ +void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf){ /*Sf[i][j] will be the integral of f on [min(i,j),max(i,j)]*/ - VALUETYPE** Sf = malloc(D*sizeof(VALUETYPE*)); + //VALUETYPE Sf[D][D]; /*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/ - VALUETYPE** Af = malloc(D*sizeof(VALUETYPE*)); - for(int i=0;i<D;i++){ + //VALUETYPE Af[D][D]; + /*Apparently may become too big for stack or something.*/ + VALUETYPE* Sf[D]; + VALUETYPE* Af[D]; + for(int i=0; i<D;i++){ Sf[i] = malloc(D*sizeof(VALUETYPE)); Af[i] = malloc(D*sizeof(VALUETYPE)); } @@ -82,7 +101,7 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ if(is_greater(Af[i][j], Mf[i])) Mf[i] = Af[i][j]; } } - + /*Print computed functions and averages*/ /* printf("Mf "); @@ -93,37 +112,37 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ printf("\n"); } */ - for(int i=0;i<D;i++){ + + for(int i=0; i<D; i++){ free(Sf[i]); free(Af[i]); } - free(Sf); - free(Af); } -void compute(int N, int K, int D, EXPTYPE p, int t){ +void compute(EXPTYPE p, int t){ /*allocate memory for f*/ - VALUETYPE* f = malloc(D*sizeof(VALUETYPE)); + VALUETYPE f[D]; /*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; /*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+K))/2+K+i] = convert_int((t >> i) & 1); + f[(D-(N+2*K))/2+K+i] = convert_int((t >> i) & 1); //if(i%3==0) f[2*N+i+K/2] = 1; } - VALUETYPE* Mf = malloc(D*sizeof(VALUETYPE)); - compute_maximalfunction(f, Mf, D); + VALUETYPE Mf[D]; + compute_maximalfunction(f, Mf); /*Allocate memory for derivatives.*/ - VALUETYPE* df = malloc(D*sizeof(VALUETYPE)); - VALUETYPE* dMf = malloc(D*sizeof(VALUETYPE)); + VALUETYPE df[D]; + VALUETYPE dMf[D]; /*Compute Kth derivative of f and Mf*/ - differentiate(f,df,D,K); - differentiate(Mf,dMf,D,K); + differentiate(f,df); + differentiate(Mf,dMf); + /*Print derivatives*/ /* @@ -138,69 +157,52 @@ void compute(int N, int K, int D, EXPTYPE p, int t){ */ /*Compute L^p norm of derivatives*/ - VALUETYPE intdfp = integratep(df,p,D); - VALUETYPE intdMfp = integratep(dMf,p,D); + VALUETYPE intdfp = integratep(df,p); + VALUETYPE intdMfp = integratep(dMf,p); //printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]); /*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/ VALUETYPE r = ratio(intdMfp, intdfp); - free(df); - free(dMf); - //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)>.6) + if(to_double(r)>.65) { printf("f: "); for(int i=0;i<D;i++) printf("%1.0f ",to_double(f[i])); printf("\n"); printf("%.4f\n",to_double(r)); } - free(f); - free(Mf); } -struct Args {int* N; int* K; int* D; EXPTYPE* p; int* t; }; +struct Args {EXPTYPE* p; int* t; }; void* perform_work(void* arguments){ struct Args *args = arguments; - int N = *(args->N); - int K = *(args->K); - int D = *(args->D); EXPTYPE p = *(args->p); int* s = args->t; while(*s <= (1 << N)-1){ int t = *s; *s = t+1; - compute(N,K,D,p,t); + compute(p,t); } return NULL; } int main() { - /*length of the support of f*/ - int N=16; - - /*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.*/ - int K=10; - - /*length of the domain, minimum N+K to make sure the support of f^{(k)} belongs to the domain*/ - int D=N+K+N/2; - /*exponent p of the L^p norm to consider*/ EXPTYPE p = 1; - int t = 1; /*Iterate over all strings of 0s and 1s with length N. Those will represent f.*/ + int t = 1; //for(int t=1; t<=(1 << N)-1; t++){ //compute(N,K,D,p,t); //} - struct Args args = {&N, &K, &D, &p, &t}; + struct Args args = {&p, &t}; pthread_t threads[NUM_THREADS]; int result_code; |