diff --git a/posts/stereo/1/.gitignore b/posts/stereo/1/.gitignore new file mode 100644 index 0000000..07337da --- /dev/null +++ b/posts/stereo/1/.gitignore @@ -0,0 +1,2 @@ +compare_complex +compare_math diff --git a/posts/stereo/1/Makefile b/posts/stereo/1/Makefile new file mode 100644 index 0000000..5cfffc1 --- /dev/null +++ b/posts/stereo/1/Makefile @@ -0,0 +1,10 @@ +default: compare_complex compare_math + +clean: + rm -rf compare_complex compare_math + +compare_complex: + gcc main.c stereo_complex.c -lm -DUSE_COMPLEX -o compare_complex + +compare_math: + gcc main.c stereo_math.c -lm -o compare_math diff --git a/posts/stereo/1/main.c b/posts/stereo/1/main.c new file mode 100644 index 0000000..074d565 --- /dev/null +++ b/posts/stereo/1/main.c @@ -0,0 +1,168 @@ +#include +#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 + +struct circle { + 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 CIRCLE_TYPE double complex +#else +#define EXTRACT_REAL(a) (a.c) +#define EXTRACT_IMAG(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; + + double total_c_error = 0, total_s_error = 0; + + size_t i; + CIRCLE_TYPE ideal, approx; + + for (i = 0; i < n; i++) { + ideal = ideals[i]; + approx = approxs[i]; + + // 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; + + if (largest_c_error < c_error) { + largest_c_error = c_error; + largest_c_index = i; + } + + 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: \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]) + ); +} + +// 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; + + clock_gettime(CLOCK_MONOTONIC, &tp1); + for (i = 0; i < n; i++) { + f(inputs[i], results + i); + } + clock_gettime(CLOCK_MONOTONIC, &tp2); + + return SECONDS_PER_NANOSECOND * (tp2.tv_sec - tp1.tv_sec) + + (tp2.tv_nsec - tp1.tv_sec); +} + +int main(int argn, char** args) +{ + long trig_time, rat_time; + + 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( +#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 + ); + + rat_time = time_computation(&approx_turn_update, rands, rats, NUM_LOOPS); + 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 + ); + } 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); + } + + return 0; +} diff --git a/posts/stereo/1/stereo_complex.c b/posts/stereo/1/stereo_complex.c index b898262..c63e4b3 100644 --- a/posts/stereo/1/stereo_complex.c +++ b/posts/stereo/1/stereo_complex.c @@ -1,20 +1,7 @@ #include #include -#include -#include -#include -#define STRRED "\x1b[31m" -#define STRGREEN "\x1b[32m" -#define STRNORM "\x1b[m" - -#define SECONDS_PER_NANOSECOND 1000000000 - -#define NUM_LOOPS 100000 - -double complex complex_turn(double turn) { return cexp(I * M_PI * turn); } - -double complex approx_turn(double turn) +double complex complex_approx_turn(double turn) { const static double a = 8 * M_SQRT2 / 3 - 3; const static double b = 4 - 8 * M_SQRT2 / 3; @@ -27,130 +14,14 @@ double complex approx_turn(double turn) return (c * c - s * s) + I * (2 * c * s); } -void print_errors( - const double* inputs, - const double complex* ideals, - const double complex* approxs, - int n -) +double complex complex_turn(double turn) { return cexp(I * M_PI * turn); } + +void turn_update(double turn, double complex* result) { - double c_error, s_error; - double largest_c_error, largest_s_error; - size_t largest_c_index, largest_s_index; - - double total_c_error = 0, total_s_error = 0; - - size_t i; - double complex ideal, approx; - - for (i = 0; i < n; i++) { - ideal = ideals[i]; - approx = approxs[i]; - - // squared error in c components - c_error = creal(ideal) - creal(approx); - c_error *= c_error; - // squared error in s components - s_error = cimag(ideal) - cimag(approx); - s_error *= s_error; - - if (largest_c_error < c_error) { - largest_c_error = c_error; - largest_c_index = i; - } - - 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: \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) * 100, largest_c_error, - sqrt(largest_c_error) * 100, inputs[largest_c_index], - creal(ideals[largest_c_index]), creal(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) * 100, largest_s_error, - sqrt(largest_s_error) * 100, inputs[largest_c_index], - cimag(ideals[largest_s_index]), cimag(approxs[largest_s_index]) - ); + *result = complex_turn(turn); } -// time the length of the computation `f` in nanoseconds -long time_computation( - double complex (*f)(double), - const double* inputs, - double complex* results, - int n -) +void approx_turn_update(double turn, double complex* result) { - size_t i; - long tick_s; - long tick_ns; - struct timespec tp; - - clock_gettime(CLOCK_MONOTONIC, &tp); - tick_ns = tp.tv_nsec; - tick_s = tp.tv_sec; - for (i = 0; i < n; i++) { - results[i] = f(inputs[i]); - } - clock_gettime(CLOCK_MONOTONIC, &tp); - - return SECONDS_PER_NANOSECOND * (tp.tv_sec - tick_s) + - (tp.tv_nsec - tick_ns); -} - -int main(int argn, char** args) -{ - long trig_time, rat_time; - - double rands[NUM_LOOPS]; - double complex trigs[NUM_LOOPS]; - double complex rats[NUM_LOOPS]; - - size_t i; - for (i = 0; i < NUM_LOOPS; i++) { - rands[i] = rand() / (double)RAND_MAX; - } - - trig_time = time_computation(&complex_turn, rands, trigs, NUM_LOOPS); - printf("Timing for %d math.h sin and cos:\t%ldns\n", NUM_LOOPS, trig_time); - - rat_time = time_computation(&approx_turn, rands, rats, NUM_LOOPS); - 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( - STRRED "math.h" STRNORM " faster, speedup: %ldns (%2.2fx)\n", diff, - frac_speed - ); - } else { - frac_speed = trig_time / (double)rat_time; - printf( - STRGREEN "Approximation" STRNORM - " faster, speedup: %ldns (%2.2fx)\n", - -diff, frac_speed - ); - print_errors(rands, trigs, rats, NUM_LOOPS); - } - - return 0; + *result = complex_approx_turn(turn); } diff --git a/posts/stereo/1/stereo_math.c b/posts/stereo/1/stereo_math.c index 2086c1f..a7025f6 100644 --- a/posts/stereo/1/stereo_math.c +++ b/posts/stereo/1/stereo_math.c @@ -1,29 +1,15 @@ #include -#include -#include -#include - -#define STRRED "\x1b[31m" -#define STRGREEN "\x1b[32m" -#define STRNORM "\x1b[m" struct circle { double c; double s; }; -double a = 8 * M_SQRT2 / 3 - 3; -double b = 4 - 8 * M_SQRT2 / 3; - -void trig(double turn, struct circle* ret) +void approx_turn_update(double turn, struct circle* ret) { - double arg = M_PI * turn; - ret->c = cos(arg); - ret->s = sin(arg); -} + static const double a = 8 * M_SQRT2 / 3 - 3; + static const double b = 4 - 8 * M_SQRT2 / 3; -void rational(double turn, struct circle* ret) -{ double p = turn * (b * turn * turn + a); double q = p * p; double r = 1 + q; @@ -33,108 +19,9 @@ void rational(double turn, struct circle* ret) ret->s = 2 * c * s; } -double -errors(int n, const struct circle* circles1, const struct circle* circles2) +void turn_update(double turn, struct circle* ret) { - double c_error, s_error; - double largest_c_error, largest_s_error; - - double total_c_error = 0, total_s_error = 0; - - int i; - struct circle circle1, circle2; - - for (i = 0; i < n; i++) { - circle1 = circles1[i]; - circle2 = circles2[i]; - - // squared error in c components - c_error = circle1.c - circle2.c; - c_error *= c_error; - // squared error in s components - s_error = circle1.s - circle2.s; - s_error *= s_error; - - if (largest_c_error < c_error) - largest_c_error = c_error; - - if (largest_s_error < s_error) - largest_s_error = s_error; - - 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: \n\tAverage: %f (%f%% error)\n\tLargest: %f " - "(%f%% error)\n", - total_c_error, sqrt(total_c_error) * 100, largest_c_error, - sqrt(largest_c_error) * 100 - ); - printf( - "Squared error in sines: \n\tAverage: %f (%f%% error)\n\tLargest: %f " - "(%f%% error)\n", - total_s_error, sqrt(total_s_error) * 100, largest_s_error, - sqrt(largest_s_error) * 100 - ); - - return 0; -} - -int main(int argn, char** args) -{ - // struct circle ret; - int i; - - struct timespec tp; - long tick; - long trig_time, rat_time; - - double rands[10000]; - struct circle trigs[10000]; - struct circle rats[10000]; - - for (i = 0; i < 10000; i++) { - rands[i] = rand() / (double)RAND_MAX; - } - - clock_gettime(CLOCK_MONOTONIC, &tp); - tick = tp.tv_nsec; - for (i = 0; i < 10000; i++) { - trig(rands[i], trigs + i); - } - clock_gettime(CLOCK_MONOTONIC, &tp); - // this isn't quite proper, since the clock may have ticked over a second - trig_time = tp.tv_nsec - tick; - printf("Timing for 10000 math.h sin and cos:\t%ldns\n", trig_time); - - clock_gettime(CLOCK_MONOTONIC, &tp); - tick = tp.tv_nsec; - for (i = 0; i < 10000; i++) { - rational(rands[i], rats + i); - } - clock_gettime(CLOCK_MONOTONIC, &tp); - rat_time = tp.tv_nsec - tick; - printf("Timing for 10000 approximations:\t%ldns\n", rat_time); - - double fracSpeed; - long linSpeed = rat_time - trig_time; - if (linSpeed > 0) { - fracSpeed = rat_time / (double)trig_time; - printf( - STRRED "math.h" STRNORM " faster, speedup: %ldns (%2.2fx)\n", - linSpeed, fracSpeed - ); - } else { - fracSpeed = trig_time / (double)rat_time; - printf( - STRGREEN "Approximation" STRNORM - " faster, speedup: %ldns (%2.2fx)\n", - -linSpeed, fracSpeed - ); - errors(10000, rats, trigs); - } + double arg = M_PI * turn; + ret->c = cos(arg); + ret->s = sin(arg); }