diff options
| -rw-r--r-- | charf.c | 66 |
1 files changed, 42 insertions, 24 deletions
@@ -12,8 +12,10 @@ #define DOUBLEERRORMODE 1 #define RATIOMODE 2 -/*Determines if computations are done with floating point type, floating point with error bounds or fractions. -Floating point without error bounds are the default.*/ +/* +Determines if computations are done with floating point type, floating point with error bounds or fractions. +Floating point without error bounds are the default. +*/ #ifndef MODE #define MODE DOUBLEMODE #endif @@ -40,17 +42,23 @@ Floating point without error bounds are the default.*/ typedef unsigned long index_t; -/*maximum length of the circle -Adjust as desired.*/ +/* +maximum length of the circle +Adjust as desired. +*/ #define N 36 -/*maximal order of derivative -Adjust as desired.*/ +/* +maximal order of derivative +Adjust as desired. +*/ #define K 64 -/*maximal number of exponents +/* +maximal number of exponents Adjust as desired. -Which exponents are considered is set using the array exponents in the beginning of the main routine.*/ +Which exponents are considered is set using the array exponents in the beginning of the main routine. +*/ #define P 5 /*Given function f on domain [0,D-1], compute derivative of order k and store result in df.*/ @@ -114,7 +122,8 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ /*Because it simplifies code later on we also assign the average on twice the circle, starting in 0.*/ Af[0][2*D] = Af[0][D]; - /*Finding the maximal average actually looks like the most costly computation, in the whole algorithm, being of order D^2. + /* + Finding the maximal average actually looks like the most costly computation, in the whole algorithm, being of order D^2. Hence we put effort into making it efficient. The strategy is to go trough all possible intervals, starting with those of largest length. For each [i,i+l) that we encounter we check for each n∈[i,i+l) if the average on [i,i+l) beats the current best average that we have so far computed for intervals containing n. @@ -137,8 +146,7 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ int intervalend = fmin(i+l,i+D); int n=i; while(n < intervalend){ - /*If the average of the current best interval for n is surely larger than A[i][l] then the same will be true for all subsequent points in that interval. - Thus, to save time just skip ahead to the end of that interval.*/ + /*If the average of the current best interval for n is surely larger than A[i][l] then the same will be true for all subsequent points in that interval. Thus, to save time just skip ahead to the end of that interval.*/ if(is_greater_certainly(Af[Pbest[n%D]][Lbest[n%D]],Af[i][l])) n += (Pbest[n%D]-n-D)%D + Lbest[n%D]; else{ /*Since the new average might win set the new best average to be the maximum of the last maximum and the new.*/ @@ -155,8 +163,7 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ n=k; } } - /*Since we may be working with errors, the above action of taking maxima may have increased the error. - That means we may have to update the average of the current interval to carry that increased error.*/ + /*Since we may be working with errors, the above action of taking maxima may have increased the error. That means we may have to update the average of the current interval to carry that increased error.*/ Af[i][l] = A; } } @@ -181,20 +188,26 @@ void compute_maximalfunction(VALUETYPE* f, VALUETYPE* Mf, int D){ } } -/*Generates all characteristic functions of length up N, indexed by i, but skip trivial ones, and pick only one representative from each class of translaties. -The return value is the length of the domain, with special values 0 which indicates we skip this function, and -1 which indicates that the index i is too large to represent a function.*/ +/* +Generates all characteristic functions of length up N, indexed by i, but skip trivial ones, and pick only one representative from each class of translaties. +The return value is the length of the domain, with special values 0 which indicates we skip this function, and -1 which indicates that the index i is too large to represent a function. +*/ int generate_each_charf(VALUETYPE* f, index_t i){ /*In order to find which function corresponds to index i, we efficiently go through all indeces up to i starting at 0.*/ index_t s=0; /*The smallest domain we consider is the circle of length 2.*/ int d=2; - /*We skip the functions that are all 0 or all 1s. + /* + We skip the functions that are all 0 or all 1s. That means every function has a point n such that f[n]=1 and f[n+1]=0. Since we also don't want to consider different translates, we may translate each function such that this value n equals d-1. - That means for each d we have to consider only arrays of length d that begin with 1 and end with 0, of which there exist 2^{d-2} many.*/ + That means for each d we have to consider only arrays of length d that begin with 1 and end with 0, of which there exist 2^{d-2} many. + */ index_t powd = 1UL << d-2; - /*At this point s is the first index that corresponds to functions on the circle of length d. - Check if our index i goes past the last index that corresponds to a function on the circle of length d.*/ + /* + At this point s is the first index that corresponds to functions on the circle of length d. + Check if our index i goes past the last index that corresponds to a function on the circle of length d. + */ while(i-s >= powd){ /*Increment s,d,powd to correspond to functions of length d+1.*/ s += powd; @@ -213,8 +226,10 @@ int generate_each_charf(VALUETYPE* f, index_t i){ Next we will perform test to determine if it should be set to false.*/ bool is_representative = true; - /*Check if t encodes a function that is made up of multiple copies of a shorter function. - We skip those because they and their maximal function are periodic and thus should be considered on a smaller circle instead.*/ + /* + Check if t encodes a function that is made up of multiple copies of a shorter function. + We skip those because they and their maximal function are periodic and thus should be considered on a smaller circle instead. + */ /*Go through all possible lenghts r of that shorter function.*/ for(int r = 1; r < d; r++){ /*r has to be a divisor of d.*/ @@ -231,7 +246,8 @@ int generate_each_charf(VALUETYPE* f, index_t i){ } } - /*We aim to filter out equivalent translates. + /* + We aim to filter out equivalent translates. The representative we pick is the one which represents the smallest integer. Every function which is a nontrivial translate of itself must in fact consist of multiplie copies of a shorter array (right?). Those we already filtered out in the previous step, so every nontrivial translate of that representative indeed represents a strictly smaller integer. @@ -345,8 +361,10 @@ void format_result(char* s, index_t index, int k, EXPTYPE p, VALUETYPE r, int fo } } -/*Given an index compute the ratio of the L^p norms of derivatives up to order k of the maximal function and the function for a given range of exponents p. -If for any order of derivative or exponent p the latest record for that ratio is broken, update the value of the record and the index of the witnessing function.*/ +/* +Given an index compute the ratio of the L^p norms of derivatives up to order k of the maximal function and the function for a given range of exponents p. +If for any order of derivative or exponent p the latest record for that ratio is broken, update the value of the record and the index of the witnessing function. +*/ int compute(index_t index, EXPTYPE exponents[P], VALUETYPE (*records_ratio)[K+1][P], index_t (*records_index)[K+1][P]){ VALUETYPE f[N]; int D = generate_function(f,index); |