diff options
| -rw-r--r-- | charf.c | 114 |
1 files changed, 76 insertions, 38 deletions
@@ -29,7 +29,7 @@ #endif /*maximum length of the support of f*/ -static const int N=16; +static const int N=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){ @@ -114,52 +114,91 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ } } -void compute(int K, EXPTYPE p, int D, VALUETYPE* f){ +void compute(EXPTYPE p, int D, VALUETYPE* f){ 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); /*Allocate memory for derivatives.*/ - VALUETYPE df[N]; - VALUETYPE dMf[N]; + VALUETYPE df[2][N]; + VALUETYPE dMf[2][N]; - /*Compute Kth derivative of f and Mf*/ - differentiate(f,df,D,K); - differentiate(Mf,dMf,D,K); - - /*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 i=0; i<D; i++){ + df[1][i] = f[i]; + df[0][i] = f[i]; + dMf[0][i] = Mf[i]; + dMf[1][i] = Mf[i]; } - */ - /*Compute L^p norm of derivatives*/ - VALUETYPE intdfp = integratep(df,p,D); - VALUETYPE intdMfp = integratep(dMf,p,D); + VALUETYPE tmp[N]; - //printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]); + for(int k=1; k<=16; k++){ + /*Compute kth derivative of f and Mf*/ + differentiate(df[(k+1)%2],df[k%2],D,1); + 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"); + } + */ - /*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/ + /*Compute L^p norm of derivatives*/ + VALUETYPE intdfp = integratep(df[k%2],p,D); + VALUETYPE intdMfp = integratep(dMf[k%2],p,D); - VALUETYPE r = ratio(intdMfp, intdfp); + //printf("%d: %f / %f = %f\n",k,intdMfp[k],intdfp[k],intdMfp[k]/intdfp[k]); - //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)>=.58) - if(to_double(r)>.53) - { - 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("%s\n",s); + /*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==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) + ){ + 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); + } } } @@ -188,13 +227,12 @@ int generate_function(VALUETYPE* f, int i){ /*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); - int K = 3; EXPTYPE p = int_to_exptype(1); VALUETYPE f[N]; int d = generate_function(f,i); while(d >= 0){ - compute(K,p,d,f); + compute(p,d,f); i += NUM_THREADS; d = generate_function(f,i); } |