When implementing sine and cosine for the bcmath version of my Gauss-Newton class, I used the Maclaurin series. Unfortunately, this ends up being slow to converge for numbers close to π. The slow convergence had required setting the bcmath scale to twice the current value in order to get accurate calculations. This made the performance of the function very slow.

I wrote about using some algebraic tricks to speed the convergence of natural log and the exponential function, and I began to think about ways I could do this with my trigonometric functions. I had investigated trying to use a modified series:

It would have allowed me to do this:

So I could just multiply the series result by the exponential of the input. However, the series just seemed like a mess and I didn't get very far with it.

My next line of thought was to see if I could use the Taylor series speed up the convergence. I had done this in the past using a look-up table for some known points to quickly calculate the Taylor series for the exponential function. Here is the general expression for the Taylor series of a given function:

In this algorithm f⁽ⁿ⁾ is the n^th derivative of the function f( x ). The variable b is some point we pick. What b allows us to do is pick a point close to x and compute the series starting there. The Maclaurin series is just a Taylor series where b is always zero. However, using the Taylor series with the value of b closer to x than x is to zero will result in the Taylor series converging faster than the Maclaurin series.

So let us apply the Taylor series for cosine with any b:

Keep in mind that cos⁽ⁿ⁾ is the n^th derivative of cosine—not the power—as we are using Lagrange's notation. At first it doesn't look like this will be too useful. Now we need cosine, all it's derivatives, and those values calculated for some b. However, this is where the properties of sine become useful. Recall the following about the sine function:

In addition, there are some special points were we get an exact real number for these trigonometric functions.

So we have four values that can be used for b and we know how the derivatives will change. Time to put this together.

Here is the PHP implementation:

The very first step is to get the input into a suitable range. This can be done quickly using the modulo operation which represents how much is left over after an integer division has been preformed.

This is used to put the input into the ±2π range. All trigonometric functions repeat themselves in this period, so it makes sense simply to narrow the input down to this range.

Then the closest b point must be chosen. This is done starting at 7/4 π and dropping by π/2 increments if the current value is greater than the input value. The 7/4 π = 2π – π/4 because we are looking for a cutoff point midway between 2π and 3/2 π. Subtracting by π/2 each time because our b points are separated by π/2. The number of subtractions are counted. In this way, b can be calculated. If n subtractions were made, then b = 2π – πn/2 = π(2 – n/2).

The next part is the derivative table.

The cells in the columns represent the function of the in the top row where x is the value in the row's left most column. Notice how the values simply shift to the left on each subsequent row, and the table repeats every 4 rows. So the number of subtractions used above to calculate b can be used as a starting index into this table.

Now the summation loop itself. Recall the Taylor series summation:

We know that cos⁽ⁿ⁾( b ) is just the one of the values from the derivative table. Calculating n! in the loop is also easy, because a_n! = a_n * a_n-1 where a₀ = 1. The same kind of trick can be used for (x – b)ⁿ. Here p_n = (x - b) * p_n-1. That makes the loop itself small and simple.

The very last trick employed in the code is for computing sine. Cosine and sine are 90º out of phase. Simply by starting at the next location in the derivative table, the exact same code for cosine can be used for sine.

I found an interesting side-effect that took a while to trace down. The resolution for the cosine operation is increased to account for rounding error during the function. This can result in it's own rounding error. The easiest example is cos( 2/3 π ). The correct answer is -0.5. However, we can not exactly calculate Pi to an infinite number of places. So the function is always passed some approximate of Pi which results in something like cos( 2/3 * 3.14159 ) = -0.499998467···. For a more accurate representation of Pi, the more decimal places until the calculation stops being 9. For cos( 2/3 * 3.1415926535 ) = -0.4999999999481578···. But the answer will only ever equal -0.5 when exactly 2/3 Pi is used. The only way I found to get around this is to per-calculate Pi to greater precision then necessary, and that does get rid of the round errors. By using 10 decimal digits of Pi, we can expect 9 digits of cosine. Then rounding after 9 digits results in cos( 2/3 * 3.1415926535 ) = -0.5 and that is what we're looking for.

Tests show that the average speed improvement of almost 240 times using the Taylor series method over the straight Maclaurin series. The function is now as fast as the others I implemented. I have to give a shout out to my calculus professor, Dr. Stredulinsky who taught me Taylor series, and an other to my discrete math and differential equations professor Dr. Fuller. These two professors have really made a difference in my understanding of mathematics. I didn't just learn the subjects to pass a test—I understood the subjects enough to apply them here.