FENIX_libc/math/cos.c
2020-12-01 17:40:03 -06:00

91 lines
3.7 KiB
C

#include <math.h>
/*
Approximates the cosine of x.
Uses a taylor polynomial (roughly) 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^20, so it
should be pretty accurate within that range.
*/
double cos(double x) {
double pi = M_PI; /* Really, me?! -Kat */
int temp;
double deg_2, deg_4, deg_6, deg_8, deg_10;
double deg_12, deg_14, deg_16, deg_18, deg_20;
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);
deg_18 = deg_16 * deg_2 * 2 / (18 * 17);
deg_20 = deg_18 * deg_2 * 2 / (20 * 19);
/*
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?
(Of course, I've not done the radius of convergence work
for this to find out how many degrees we need to get a
radius of convergence of pi, so I may have used too many
terms, anyways. Not that I can't reduce that...)
*/
cosine = -1 + deg_2 - deg_4 + deg_6 - deg_8 + deg_10;
cosine = cosine - deg_12 + deg_14 - deg_16 + deg_18 + deg_20;
return cosine;
}