summaryrefslogtreecommitdiff
path: root/charf.c
diff options
context:
space:
mode:
Diffstat (limited to 'charf.c')
-rw-r--r--charf.c129
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
generated by cgit v1.2.3 (git 2.25.1) at 2026-02-07 00:49:49 +0000