diff --git a/posts/stereo/1/main.c b/posts/stereo/1/main.c index 074d565..28d2768 100644 --- a/posts/stereo/1/main.c +++ b/posts/stereo/1/main.c @@ -3,166 +3,172 @@ #include #include #include +#include #define STR_RED "\x1b[31m" #define STR_GREEN "\x1b[32m" #define STR_NORM "\x1b[m" #define SECONDS_PER_NANOSECOND 1000000000 -#define NUM_LOOPS 100000 -#define SQRT_NUM_LOOPS 100 +#define NUM_STATS 100000 +#define SQRT_NUM_STATS 100 +#define NUM_LOOPS 10000000 struct circle { - double c; - double s; + double c; + double s; }; -void turn_update(double turn, void* result); -void approx_turn_update(double turn, void* result); - #ifdef USE_COMPLEX -#define EXTRACT_REAL(a) creal(a) -#define EXTRACT_IMAG(a) cimag(a) +#define EXTRACT_COSINE(a) creal(a) +#define EXTRACT_SINE(a) cimag(a) #define CIRCLE_TYPE double complex #else -#define EXTRACT_REAL(a) (a.c) -#define EXTRACT_IMAG(a) (a.s) +#define EXTRACT_COSINE(a) (a.c) +#define EXTRACT_SINE(a) (a.s) #define CIRCLE_TYPE struct circle #endif -void print_errors( - const double* inputs, - const CIRCLE_TYPE* ideals, - const CIRCLE_TYPE* approxs, - int n -) -{ - double c_error, s_error; - double largest_c_error, largest_s_error; - size_t largest_c_index, largest_s_index; +void turn_update(double turn, CIRCLE_TYPE *result); +void approx_turn_update(double turn, CIRCLE_TYPE *result); - double total_c_error = 0, total_s_error = 0; +void print_errors(const double *inputs, const CIRCLE_TYPE *ideals, + const CIRCLE_TYPE *approxs, int n) { + double c_error, s_error; + double largest_c_error = 0, largest_s_error = 0; + size_t largest_c_index, largest_s_index; - size_t i; - CIRCLE_TYPE ideal, approx; + double total_c_error = 0, total_s_error = 0; - for (i = 0; i < n; i++) { - ideal = ideals[i]; - approx = approxs[i]; + size_t i; + CIRCLE_TYPE ideal, approx; - // squared error in c components - c_error = EXTRACT_REAL(ideal) - EXTRACT_REAL(approx); - c_error *= c_error; - // squared error in s components - s_error = EXTRACT_IMAG(ideal) - EXTRACT_IMAG(approx); - s_error *= s_error; + for (i = 0; i < n; i++) { + ideal = ideals[i]; + approx = approxs[i]; - if (largest_c_error < c_error) { - largest_c_error = c_error; - largest_c_index = i; - } + // squared error in c components + c_error = EXTRACT_COSINE(ideal) - EXTRACT_COSINE(approx); + c_error *= c_error; + // squared error in s components + s_error = EXTRACT_SINE(ideal) - EXTRACT_SINE(approx); + s_error *= s_error; - if (largest_s_error < s_error) { - largest_s_error = s_error; - largest_s_index = i; - } - - total_c_error += c_error; - total_s_error += s_error; + if (largest_c_error < c_error) { + largest_c_error = c_error; + largest_c_index = i; } - // these now contain the *average* squared error - total_c_error /= (double)n; - total_s_error /= (double)n; - printf( - "Squared error in cosines: \n" - "\tAverage: %f (%f%% error)\n" - "\tLargest: %f (%f%% error)\n" - "\t\tInput:\t\t%f\n" - "\t\tValue:\t\t%f\n" - "\t\tApproximation:\t%f\n", - total_c_error, sqrt(total_c_error) * SQRT_NUM_LOOPS, largest_c_error, - sqrt(largest_c_error) * SQRT_NUM_LOOPS, inputs[largest_c_index], - EXTRACT_REAL(ideals[largest_c_index]), - EXTRACT_REAL(approxs[largest_c_index]) - ); - printf( - "Squared error in sines: \n" - "\tAverage: %f (%f%% error)\n" - "\tLargest: %f (%f%% error)\n" - "\t\tInput:\t\t%f\n" - "\t\tValue:\t\t%f\n" - "\t\tApproximation:\t%f\n", - total_s_error, sqrt(total_s_error) * SQRT_NUM_LOOPS, largest_s_error, - sqrt(largest_s_error) * SQRT_NUM_LOOPS, inputs[largest_c_index], - EXTRACT_IMAG(ideals[largest_s_index]), - EXTRACT_IMAG(approxs[largest_s_index]) - ); + if (largest_s_error < s_error) { + largest_s_error = s_error; + largest_s_index = i; + } + + total_c_error += c_error; + total_s_error += s_error; + } + // these now contain the *average* squared error + total_c_error /= (double)n; + total_s_error /= (double)n; + + printf("Squared error in cosines (%d runs): \n" + "\tAverage: %f (%f%% error)\n" + "\tLargest: %f (%f%% error)\n" + "\t\tInput:\t\t%f\n" + "\t\tValue:\t\t%f\n" + "\t\tApproximation:\t%f\n", + NUM_STATS, total_c_error, sqrt(total_c_error) * SQRT_NUM_STATS, + largest_c_error, sqrt(largest_c_error) * SQRT_NUM_STATS, + inputs[largest_c_index], EXTRACT_COSINE(ideals[largest_c_index]), + EXTRACT_COSINE(approxs[largest_c_index])); + printf("Squared error in sines (%d runs): \n" + "\tAverage: %f (%f%% error)\n" + "\tLargest: %f (%f%% error)\n" + "\t\tInput:\t\t%f\n" + "\t\tValue:\t\t%f\n" + "\t\tApproximation:\t%f\n", + NUM_STATS, total_s_error, sqrt(total_s_error) * SQRT_NUM_STATS, + largest_s_error, sqrt(largest_s_error) * SQRT_NUM_STATS, + inputs[largest_s_index], EXTRACT_SINE(ideals[largest_s_index]), + EXTRACT_SINE(approxs[largest_s_index])); } // time the length of the computation `f` in nanoseconds -long time_computation( - void (*f)(double, void*), const double* inputs, CIRCLE_TYPE* results, int n -) -{ - size_t i; - struct timespec tp1; - struct timespec tp2; +// using a macro rather than a function taking a function pointer for +// more accurate time +#define TIME_COMPUTATION(f, inputs, results, n_stats, n_loop, time) \ + do { \ + size_t i; \ + struct timespec tp1; \ + struct timespec tp2; \ + \ + for (i = 0; i < n_stats; i++) { \ + f(inputs[i], results + i); \ + } \ + CIRCLE_TYPE temp; \ + clock_gettime(CLOCK_MONOTONIC, &tp1); \ + for (i = 0; i < n_stats; i++) { \ + f(rand() / (double)RAND_MAX, &temp); \ + } \ + clock_gettime(CLOCK_MONOTONIC, &tp2); \ + \ + time = SECONDS_PER_NANOSECOND * (tp2.tv_sec - tp1.tv_sec) + \ + (tp2.tv_nsec - tp1.tv_nsec); \ + } while (0) - clock_gettime(CLOCK_MONOTONIC, &tp1); - for (i = 0; i < n; i++) { - f(inputs[i], results + i); - } - clock_gettime(CLOCK_MONOTONIC, &tp2); +int main(int argn, char **args) { + long trig_time, rat_time; - return SECONDS_PER_NANOSECOND * (tp2.tv_sec - tp1.tv_sec) + - (tp2.tv_nsec - tp1.tv_sec); -} + double rands[NUM_STATS]; + CIRCLE_TYPE trigs[NUM_STATS]; + CIRCLE_TYPE rats[NUM_STATS]; -int main(int argn, char** args) -{ - long trig_time, rat_time; + size_t i; + for (i = 0; i < NUM_STATS; i++) { + rands[i] = rand() / (double)RAND_MAX; + } - double rands[NUM_LOOPS]; - CIRCLE_TYPE trigs[NUM_LOOPS]; - CIRCLE_TYPE rats[NUM_LOOPS]; - - size_t i; - for (i = 0; i < NUM_LOOPS; i++) { - rands[i] = rand() / (double)RAND_MAX; - } - - trig_time = time_computation(&turn_update, rands, trigs, NUM_LOOPS); - printf( + TIME_COMPUTATION(turn_update, rands, trigs, NUM_STATS, NUM_LOOPS, trig_time); + printf( #ifdef USE_COMPLEX "Timing for %d complex.h cexp:\t%ldns\n", #else "Timing for %d math.h sin and cos:\t%ldns\n", #endif - NUM_LOOPS, - trig_time - ); + NUM_LOOPS, trig_time); - rat_time = time_computation(&approx_turn_update, rands, rats, NUM_LOOPS); - printf("Timing for %d approximations:\t%ldns\n", NUM_LOOPS, rat_time); + TIME_COMPUTATION(approx_turn_update, rands, rats, NUM_STATS, NUM_LOOPS, + rat_time); + printf("Timing for %d approximations:\t%ldns\n", NUM_LOOPS, rat_time); - long diff = rat_time - trig_time; - double frac_speed; - if (diff > 0) { - frac_speed = rat_time / (double)trig_time; - printf( - STR_RED "stdlib" STR_NORM " faster, speedup: %ldns (%2.2fx)\n", - diff, frac_speed - ); + // Report results + long diff = rat_time - trig_time; + double frac_speed; + if (diff > 0) { + frac_speed = rat_time / (double)trig_time; + + // Disable colors for non-terminal output + if (isatty(STDOUT_FILENO)) { + printf(STR_RED "stdlib" STR_NORM " faster, speedup: %ldns (%2.2fx)\n", + diff, frac_speed); } else { - frac_speed = trig_time / (double)rat_time; - printf( - STR_GREEN "Approximation" STR_NORM - " faster, speedup: %ldns (%2.2fx)\n", - -diff, frac_speed - ); - print_errors(rands, trigs, rats, NUM_LOOPS); + printf("stdlib faster, speedup: %ldns (%2.2fx)\n", diff, frac_speed); + } + } else { + frac_speed = trig_time / (double)rat_time; + + // Disable colors for non-terminal output + if (isatty(STDOUT_FILENO)) { + printf(STR_GREEN "Approximation" STR_NORM + " faster, speedup: %ldns (%2.2fx)\n", + -diff, frac_speed); + } else { + printf("Approximation faster, speedup: %ldns (%2.2fx)\n", -diff, + frac_speed); } - return 0; + print_errors(rands, trigs, rats, NUM_STATS); + } + + return 0; }