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

#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;
}