#include <math.h>

/*
  Approximates the cosine of x.

  Uses a taylor polynomial accurate between 0 and 2*pi to approximate the 
  cosine. Any values of x outside this range are modulo'd to be within this 
  range (you know, 'cos cosine is cyclic).

  The taylor poly is centered at pi, with a radius of
  convergence of no less than pi, making it roughly accurate 
  between 0 and 2 * pi. The polynomial is up to x^16, which I've
  confirmed as being accurate being 0 and 2 * pi.
*/
double cos(double x) {
  double pi = M_PI; /* Really, me?! -Kat */
  int temp;
  double deg_2, deg_4, deg_6, deg_8, deg_10, deg_12, deg_14, deg_16;
  double cosine;
  
  if(x < 0) x = -x;

  /* Hey look, an fmod implementation! -Kat */
  if(x >= 2 * pi || x <= -2 * pi) {
    temp = x / (2 * pi);
    x -= temp * (2 * pi);
  }

  /* 
     Terms of our taylor poly, split up to make writing the
     actual polynomial less hellish. If this were a Maclaurin
     polynomial, I could maybe just write it, but having to deal
     this (x-pi) nonse, this is easier, I feel.
     
     This might be worth a rewrite once I get double pow(double)
     written, to use pow(x - pi, y), where y is the degree. Also,
     if math.h or something has a factorial function, it might
     be good to use that for the denominator.
  */
  deg_2 = (x - pi) * (x - pi) / 2;
  deg_4 = deg_2 * deg_2 * 2 / (4 * 3);
  deg_6 = deg_4 * deg_2 * 2 / (6 * 5);
  deg_8 = deg_6 * deg_2 * 2 / (8 * 7);
  deg_10 = deg_8 * deg_2 * 2 / (10 * 9);
  deg_12 = deg_10 * deg_2 * 2 / (12 * 11);
  deg_14 = deg_12 * deg_2 * 2 / (14 * 13);
  deg_16 = deg_14 * deg_2 * 2 / (16 * 15);

  /*
     In case you aren't familiar with the theory of a Taylor
     polynomial, the basic idea is that *any* differentiable
     function can be written as a polynomial of some length.
     I won't go too much into how we get to this polynomial
     (it involves a lot of calculus), but if we carry out this
     process for cos(x), we get a polynomial of roughly the
     form of the sum from k=0 to infinity of: 
         -1^(k+1) * x^(2*k) / (2k)!
     So, for the 0th term, it's -1^1 * x^0 / 0! or -1/1 = -1.
     For the 1st term, it's -1^2 * x^2 / 2! or x^2 / 2.
     For term 2: -1^3 * x^4 / 4! or -(x^4) / 24.
     So on, so forth. If we want it centered at pi, as we do
     for this, we need to instead calculate that sum for
     (x - pi) instead of just x. Also, notice that we don't
     calculate an infinite sum. Just 11 terms. This is a result
     of cosine's cyclical nature. Every possible value can be
     given by just modulo-ing whatever we want to know the
     cosine of to be between 0 and 2 * pi. So, for instance,
     the cosine of 5 * pi / 2 is equal to the cosine of pi / 2.
     So, we only need it to be accurate between 0 and 2 * pi.
     So, why center on pi? To reduce the number of terms. If
     we center on zero, we have to stretch the radius of
     convergence from 0 to 2 * pi. But if we center at pi,
     that radius only needs to be from 0 to pi and pi to 2 * pi,
     giving a radius of convergence of only pi instead of 2 *pi.
     Plus, since we're taking advantage of cosine being an
     even function, we don't need any terms less than 0, so
     why *center* at zero and include negatives?
  */
  cosine = -1 + deg_2 - deg_4 + deg_6 - deg_8 + deg_10 - deg_12 + deg_14 - deg_16;

  return cosine;
}