diff options
| -rw-r--r-- | charf.c | 56 |
1 files changed, 46 insertions, 10 deletions
@@ -3,16 +3,20 @@ #include <math.h> +/*given function f on domain [0,D-1], compute derivatives f^{(0)} until f^{(K)} and store them in df*/ void differentiate(double* f, double** df, int D, int K){ df[0] = f; int 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(i=0; i<D-k; i++){ df[k][i] = df[k-1][i+1] - df[k-1][i]; } } } +/*given function f on domain [0,D-1] compute pth root of integral of |f|^p*/ double integrate(double* f, double p, int D){ double sum = 0.0; for(int i=0;i<D;i++){ @@ -21,16 +25,22 @@ double integrate(double* f, double p, int D){ return pow(sum,1/p); } -double compute_derivatives(int* f, int** Sf, double** Af, double* Mf, double** df, double** dMf, double p, double* intdf, double* intdMf, int N, int D, int K){ +/*given integer valued function f on domain D, compute ||Mf^{(K)}||_p/||f^{(K)}||_p*/ +/*All the other arguments are pointers to storage for intermediate variables. The purpose is that we do not have to allocate new storage with each invocation. Probably there's a more user friendly way but I don't know much about garbage collection.*/ +double compute_derivatives(int* f, int** Sf, double** Af, double* Mf, double** df, double** dMf, double p, double* intdf, double* intdMf, int D, int K){ +/*Convert integer valued f to float valued df[0]*/ for(int i=0;i<D;i++) df[0][i] = (double) f[i]; -int r=0; +/*Sf[i][j] will be the integral of f on [min(i,j),max(i,j)]*/ +/*Af[i][j] will be the average of f on [min(i,j),max(i,j)]*/ for(int i=0;i<D;i++) { Sf[i][i] = f[i]; Af[i][i] = 1.0*Sf[i][i]; } + 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] = Sf[i][i+n-1]+f[i+n]; Sf[i+n][i] = Sf[i][i+n]; @@ -40,6 +50,7 @@ for(int n=1;n<D;n++){ } } +/*Compute maximal function by picking the largest average*/ for(int i=0;i<D;i++){ Mf[i] = Af[i][i]; for(int j=0;j<D;j++){ @@ -47,17 +58,20 @@ for(int i=0;i<D;i++){ } } +/*Print computed functions and averages*/ //printf("Mf "); //for(int i=0;i<D;i++) printf("%+0.1f ",Mf[i]); //printf("\n"); -for(int i=0;i<D;i++){ +//for(int i=0;i<D;i++){ //for(int j=0;j<D;j++) printf("%0.1f ",Af[i][j]); //printf("\n"); -} +//} +/*Compute 0th until Kth derivative of f and Mf*/ differentiate(df[0],df,D,K); differentiate(dMf[0],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]); @@ -68,24 +82,37 @@ differentiate(dMf[0],dMf,D,K); //} //for(int k=0;k<=K;k++){ + +/*Compute L^p norm of derivatives*/ int k=K; - intdf[k] = integrate(df[k],p,D); - intdMf[k] = integrate(dMf[k],p,D); - //printf("%d: %f / %f = %f\n",k,intdMf[k],intdf[k],intdMf[k]/intdf[k]); +intdf[k] = integrate(df[k],p,D); +intdMf[k] = integrate(dMf[k],p,D); +//printf("%d: %f / %f = %f\n",k,intdMf[k],intdf[k],intdMf[k]/intdf[k]); +//} +/*Return ratio of L^p norms*/ return intdMf[k]/intdf[k]; -//} } int main() { +/*length of the support of f*/ int N=16; + +/*length of the domain*/ int D=5*N; + +/*order of the derivative to consider*/ int K=3; + +/*exponent p of the L^p norm to consider*/ double p = 1.0; +/*allocate memory for f*/ int* f = malloc(D*sizeof(int)); + +/*allocate memory for temporary variables, such as averages*/ int** Sf = malloc(D*sizeof(int*)); double** Af = malloc(D*sizeof(double*)); double* Mf = malloc(D*sizeof(double)); @@ -94,9 +121,9 @@ for(int i=0;i<D;i++){ Af[i] = malloc(D*sizeof(double)); } +/*Allocate memory for derivatives*/ double** df = malloc(D*sizeof(double*)); double** dMf = malloc(D*sizeof(double*)); - for(int k=0;k<=K;k++){ df[k] = malloc(D*sizeof(double)); dMf[k] = malloc(D*sizeof(double)); @@ -107,19 +134,28 @@ for(int k=0;k<=K;k++){ } dMf[0] = Mf; +/*Allocate memory for integrals*/ double* intdf = malloc(sizeof(double)*K); double* intdMf = malloc(sizeof(double)*K); + +/*Allocate memory for ||Mf^{(k)}||_p/||f^{(k)}||_p.*/ double r = 0.0; +/*Iterate over all strings of 0s and 1s with length N. Those will represent f.*/ for(int t=1; t<=(1 << N)-1; t++){ + /*Initiate f to be zero everywhere.*/ for(int i=0;i<D;i++) f[i]=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/2 so that the most interesting part of f^{(K)} and Mf^{(K)} will be around the center of the domain*/ f[2*N+i+K/2] = (t >> i) & 1; //if(i%3==0) f[2*N+i+K/2] = 1; } - r = compute_derivatives(f, Sf, Af, Mf, df, dMf, p, intdf, intdMf, N, D, K); + /*Compute ||Mf^{(k)}||_p/||f^{(k)}||_p.*/ + r = compute_derivatives(f, Sf, Af, Mf, df, dMf, p, intdf, intdMf, D, 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(r>.4997){ printf("f: "); for(int i=0;i<D;i++) printf("%.1d ",f[i]); |