zenzicubi.co/posts/stereo/1/main.c

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>

#define STR_RED "\x1b[31m"
#define STR_GREEN "\x1b[32m"
#define STR_NORM "\x1b[m"

#define SECONDS_PER_NANOSECOND 1000000000
#define NUM_STATS 100000
#define SQRT_NUM_STATS 100
#define NUM_LOOPS 10000000

struct circle {
  double c;
  double s;
};

#ifdef USE_COMPLEX
#define EXTRACT_COSINE(a) creal(a)
#define EXTRACT_SINE(a) cimag(a)
#define CIRCLE_TYPE double complex
#else
#define EXTRACT_COSINE(a) (a.c)
#define EXTRACT_SINE(a) (a.s)
#define CIRCLE_TYPE struct circle
#endif

void turn_update(double turn, CIRCLE_TYPE *result);
void approx_turn_update(double turn, CIRCLE_TYPE *result);

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;

  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_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_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 (%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
// 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)

int main(int argn, char **args) {
  long trig_time, rat_time;

  double rands[NUM_STATS];
  CIRCLE_TYPE trigs[NUM_STATS];
  CIRCLE_TYPE rats[NUM_STATS];

  size_t i;
  for (i = 0; i < NUM_STATS; i++) {
    rands[i] = rand() / (double)RAND_MAX;
  }

  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);

  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);

  // 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 {
      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);
    }

    print_errors(rands, trigs, rats, NUM_STATS);
  }

  return 0;
}