diff options
| -rw-r--r-- | charf.c | 129 |
1 files changed, 79 insertions, 50 deletions
@@ -33,10 +33,13 @@ #endif /*maximum length of the support of f*/ -static const int N=24; +#define N 24 /*maximal order of derivative*/ -static const int K=24; +#define K 24 + +/*maximal number of exponents*/ +#define P 1 /*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){ @@ -159,23 +162,28 @@ int generate_function(VALUETYPE* f, int i){ //return generate_triangle(f,i); } -void format_result(char* s, EXPTYPE p, int D, int index, int k, VALUETYPE r){ +void format_result(char* s, int index, int k, EXPTYPE p, VALUETYPE r){ VALUETYPE f[N]; int d = generate_function(f,index); sprintf(s,"f: "); int l = 3; - for(int i=0; i<D; i++) l += sprintf(s+l,"%2.0f ",to_double(f[i])); - for(int i=D; i<N; i++) l += sprintf(s+l," "); + for(int i=0; i<d; i++) l += sprintf(s+l,"%2.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); + l += sprintf(s+l,"%dth der.: %s",k,rts); } bool over_threshold_charf(double t, int k){ return (k==1 && t>=1) || (k==2 && t>=.5) || (k==3 && t>.53) || (k==4 && t>=.5) || (k==5 && t>=.58) || (k==6 && t>=.58) || (k==7 && t>=.69) || (k==8 && t>=.83) || (k==9 && t>=.8699) || (k==10 && t>=.919) || (k==11 && t>=.97) || (k==12 && t>=.97) || (k==13 && t>=.98) || (k==14 && t>=.98) || (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) ; } -bool compute(EXPTYPE p, int D, VALUETYPE* f, int index, VALUETYPE* records_ratio, int* records_index){ +int compute(int index, int num_exponents[K+1], EXPTYPE exponents[K+1][P], VALUETYPE (*records_ratio)[K+1][P], int (*records_index)[K+1][P]){ + VALUETYPE f[N]; + int D = generate_function(f,index); + /*Immediately abort if index is out of bounds.*/ + if(D <= 0) return D; + 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); @@ -199,93 +207,112 @@ bool compute(EXPTYPE p, int D, VALUETYPE* f, int index, VALUETYPE* records_ratio differentiate(dMf[(k+1)%2],dMf[k%2],D,1); } - /*Print derivatives*/ - /* - for(int k=0; k<=K; k++){ - printf("f %d: ",k); - 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]); - printf("\n"); - } - */ + for(int p=0; p<num_exponents[k]; p++){ + /*Print derivatives*/ + /* + for(int k=0; k<=K; k++){ + printf("f %d: ",k); + 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]); + printf("\n"); + } + */ - /*Compute L^p norm of derivatives*/ - VALUETYPE intdfp = integratep(df[k%2],p,D); - VALUETYPE intdMfp = integratep(dMf[k%2],p,D); + /*Compute L^p norm of derivatives*/ + VALUETYPE intdfp = integratep(df[k%2],exponents[k][p],D); + VALUETYPE intdMfp = integratep(dMf[k%2],exponents[k][p],D); - //printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]); + //printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]); - /*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/ + /*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/ - VALUETYPE r = ratio(intdMfp, intdfp); + VALUETYPE r = ratio(intdMfp, intdfp); - if( - //over_threshold_charf(to_double(r),k) && - is_greater_possibly(r,records_ratio[k])){ - /*extra check for printing only because in error mode for some reason floats randomly seem to increase by tiny amounts*/ - if(is_greater_certainly(r,records_ratio[k])){ - records_index[k] = index; - char s[1024]; - format_result(s, p, D, index, k, r); - printf("%s",s); + if( + //over_threshold_charf(to_double(r),k) && + is_greater_possibly(r,(*records_ratio)[k][p])){ + /*extra check for printing only because in error mode for some reason floats randomly seem to increase by tiny amounts*/ + if(is_greater_certainly(r,(*records_ratio)[k][p])){ + (*records_index)[k][p] = index; + char s[1024]; + format_result(s, index, k, exponents[k][p], r); + printf("%s\n",s); + } + (*records_ratio)[k][p] = maximum((*records_ratio)[k][p],r); } - records_ratio[k] = maximum(records_ratio[k],r); } } + return D; +} + +void print_records(int num_exponents[K+1], EXPTYPE exponents[K+1][P], VALUETYPE ratios[K+1][P], int indeces[K+1][P]){ + char s[1024]; + printf("Current records:\n"); + for(int k=0; k<=K; k++) for(int p=0; p<num_exponents[k]; p++){ + format_result(s, indeces[k][p], k, exponents[k][p], ratios[k][p]); + printf("%s\n",s); + } } -typedef struct { VALUETYPE* records_ratio; int* records_index; int num_thread; int* domain_current; } Args; +typedef struct { int num_exponents[K+1]; EXPTYPE exponents[K+1][P]; VALUETYPE (*records_ratio)[K+1][P]; int (*records_index)[K+1][P]; 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){ Args* args = arguments; int i = args -> num_thread; - VALUETYPE* records_ratio = args -> records_ratio; - int* records_index = args -> records_index; 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); + int d = 0; while(d >= 0){ + d = compute(i, args -> num_exponents, args -> exponents, args -> records_ratio, args -> records_index); if(d > *domain_current){ *domain_current = d; printf("Start considering length: %d\n",d); } - compute(p,d,f,i,records_ratio,records_index); i += NUM_THREADS; - d = generate_function(f,i); } } int main() { + int num_exponents[K+1]; + EXPTYPE exponents[K+1][P]; + for(int k=0; k<=K; k++){ + num_exponents[k] = 1; + exponents[k][0] = int_to_exptype(1); + } + pthread_t threads[NUM_THREADS]; Args args[NUM_THREADS]; - VALUETYPE records_ratio[K+1]; - int records_index[K+1]; + VALUETYPE records_ratio[K+1][P]; + int records_index[K+1][P]; + for(int k=0; k<=K; k++) for(int p=0; p<num_exponents[k]; p++) records_ratio[k][p] = int_to_valuetype(0); int domain_current = 0; - for(int k=0; k<=K; k++) records_ratio[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].records_ratio = records_ratio; - args[i].records_index = records_index; + memcpy(args[i].num_exponents,num_exponents,sizeof num_exponents); + memcpy(args[i].exponents,exponents,sizeof exponents); + args[i].records_ratio = &records_ratio; + args[i].records_index = &records_index; args[i].num_thread = i; args[i].domain_current = &domain_current; result_code = pthread_create(&threads[i], NULL, compute_chunk, args+i); } size_t sz = 1024; - char prmpt[] = "To exit enter q,quit or exit.\nTo print current records to terminal enter r.\nTo print current records to a file enter p.\nTo do so with tex formatting enter t.\n"; + char prmpt[] = "To exit enter q, quit or exit.\nTo print current records to terminal enter r.\nTo print current records to a file enter p.\nTo do so with tex formatting enter t.\n"; char buff[sz]; bool cont = true; do { getLine(prmpt, buff, sz); - if( 0 == strcmp(buff,"q") || 0 == strcmp(buff,"quit") || 0 == strcmp(buff,"exit") ) exit(EXIT_SUCCESS); + if( 0 == strcmp(buff,"q") || 0 == strcmp(buff,"quit") || 0 == strcmp(buff,"exit") ){ + print_records(num_exponents,exponents,records_ratio,records_index); + exit(EXIT_SUCCESS); + } + else if( 0 == strcmp(buff,"r")) print_records(num_exponents,exponents,records_ratio,records_index); } while(cont); // wait for each thread to complete @@ -294,5 +321,7 @@ int main() { printf("In main: Thread %d has ended.\n", i); } + print_records(num_exponents,exponents,records_ratio,records_index); + return 0; } |