diff --git a/posts/stereo/1/stereo_complex.c b/posts/stereo/1/stereo_complex.c new file mode 100644 index 0000000..a3a42f9 --- /dev/null +++ b/posts/stereo/1/stereo_complex.c @@ -0,0 +1,143 @@ +#include +#include +#include +#include +#include + +#define STRRED "\x1b[31m" +#define STRGREEN "\x1b[32m" +#define STRNORM "\x1b[m" + +#define NUM_LOOPS 100000 + +double complex complex_turn(double turn) +{ + return cexp(I*M_PI*turn); +} + +double complex approx_turn(double turn) +{ + const static double a = 8*M_SQRT2/3 - 3; + const static double b = 4 - 8*M_SQRT2/3; + + double p = turn*(b * turn * turn + a); + double q = p*p; + double r = 1 + q; + double c = (1 - q) / r, s = 2*p / r; + + 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 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])); +} + +//time the length of the computation `f` in nanoseconds +long time_computation( double complex (*f)(double), + const double *inputs, + double complex *results, + int n) +{ + size_t i; + long tick; + struct timespec tp; + + clock_gettime(CLOCK_MONOTONIC, &tp); + tick = tp.tv_nsec; + for (i = 0; i < n; i++) { + results[i] = f(inputs[i]); + } + //this isn't quite proper, since the clock may have ticked over a second + clock_gettime(CLOCK_MONOTONIC, &tp); + + return tp.tv_nsec - tick; +} + +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; +} diff --git a/posts/stereo/1/stereo_math.c b/posts/stereo/1/stereo_math.c new file mode 100644 index 0000000..d841204 --- /dev/null +++ b/posts/stereo/1/stereo_math.c @@ -0,0 +1,128 @@ +#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) +{ + double arg = M_PI * turn; + ret->c = cos(arg); + ret->s = sin(arg); +} + +void rational(double turn, struct circle *ret) +{ + double p = turn*(b * turn * turn + a); + double q = p*p; + double r = 1 + q; + double c = (1 - q) / r, s = 2*p / r; + + ret->c = c*c - s*s; + ret->s = 2*c*s; +} + +double errors(int n, const struct circle *circles1, const struct circle *circles2) +{ + 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); + } +}