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(n) is the nth 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(n) is the nth 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:
/*
* Computes cosine of x.
* @author Andrew Que.
* @license Public domain.
*/
function bccos( $x, $isSine = false )
{
// Cosine and sine repeat in intervals of 2*pi >= x >= -2*pi. So reduce
// input down to this range.
// To do this, get the whole number of how many times two pi divides into the
// value. Then subtract off the whole number part times two pi. This is a
// modulus remainder, but bcmod only returns the integer part of the
// remainder.
$twoPi = bcmul( bcpi(), 2 );
$mod = bcdiv( $x, $twoPi );
list( $whole, $fractional ) = explode( '.', $mod );
$x = bcsub( $x, bcmul( $whole, $twoPi ) );
$correction = bcpi();
// Modify the scale for additional accuracy.
// Done after the interval has been reduced to avoid rounds errors with pi.
$scale = bcGetScale();
bcscale( $scale + 5 );
//-----------------------------------
// The Taylor series for cosine is as follows:
// inf
// ---
// \ d^n (x - b)^n
// / ---- cos(x) ---------
// --- dx^n n!
// n=0
// We use the Taylor series to compute cosine. We could use the Maclaurin
// series, but this converges slowly at points close to pi. So we center
// on one of four points: 0, pi/2, pi, and 3/2 pi. This is done because the
// nth derivative is easy to compute at these points: it simply alternates
// between 1, 0 and -1. Where it starts depends on what value the input is
// closest to.
//-----------------------------------
// First, figure out which point the input value is closest to.
// NOTE: Regular numbers used here, not BC numbers. No reason to be that
// that accurate yet.
$taylorIndex = 0;
$taylorPoint = pi() * 7/4;
$halfPi = pi() / 2;
while ( $taylorPoint >= $x )
{
$taylorPoint -= $halfPi;
$taylorIndex += 1;
}
// Figure out what point we are starting.
// Now the value is needed as a BC number.
$taylorPoint = 2 - $taylorIndex / 2;
$taylorPoint = bcmul( $taylorPoint, bcpi() );
// The index for what the nth derivative will be.
// Note that $taylorIndex is subtracted by one in order to iterate it. This
// is the same as adding 3 and doing a modulus by 4.
$taylorIndex %= 4;
// To do sine rather than cosine, simply add one to the starting index.
if ( $isSine )
$taylorIndex = ( $taylorIndex + 1 ) % 4;
// This is a look-up table for the derivatives.
$taylorIndexMap = array( 1, 0, -1, 0 );
// Setup the series variables.
$x = bcsub( $x, $taylorPoint );
$power = 1;
$newResult = 0;
$result = -1; // <- Something not equal to $newResult
$n = 0;
$fact = 1;
while ( bccomp( $newResult, $result ) )
{
// Get the nth derivative of the Taylor point.
$derivative = $taylorIndexMap[ $taylorIndex ];
// Don't bother doing any addition of derivative is zero.
if ( 0 != $derivative )
{
$result = $newResult;
if ( -1 == $derivative )
$accumulator = bcneg( $power );
else
$accumulator = $power;
$accumulator = bcdiv( $accumulator, $fact );
$newResult = bcadd( $result, $accumulator );
}
$n += 1;
$power = bcmul( $power, $x );
$fact = bcmul( $fact, $n );
$taylorIndex = ( $taylorIndex + 3 ) % 4;
}
// Cut off the last few used to account for accumulated error--they are
// incorrect anyway.
$result = bcround( $result, $scale - 1 );
bcscale( $scale );
return $result;
}
/*
* Computes cosine of x.
* @author Andrew Que.
* @license Public domain.
*/
function bcsin( $x )
{
return bccos( $x, true );
}
/*
* Return the BC math scale.
*
* This function will return the last value passed to bcscale.
*
* @author Andrew Que.
* @license Public domain.
*/
function bcGetScale()
{
return strlen( bcadd( 0, 0 ) ) - 2;
}
/*
* Rounds a value to the number of places requested.
* @author Andrew Que.
* @license Public domain.
*/
function bcRound( $value, $places )
{
// Round result by adding 5x10^(-places).
if ( -1 == bccomp( $value, "0" ) )
$value = bcadd( $value, "-0." . str_pad( "", $places, "0" ) . "5" );
else
$value = bcadd( $value, "0." . str_pad( "", $places, "0" ) . "5" );
// Split number into the whole and fractional parts.
list( $integer, $fractional ) = explode( ".", $value );
// If are result uses all the digits, cut off the last few used to account
// for rounding error--they are incorrect anyway.
if ( strlen( $fractional ) > $places )
{
$fractional = substr( $fractional, 0, $places );
$value = $integer . "." . $fractional;
}
return $value;
}
/*
* Negate value.
* @author Andrew Que.
* @license Public domain.
*/
function bcneg( $value )
{
// Already negative?
if ( '-' == $value[ 0 ] )
$value = substr( $value, 1 );
else
$value = "-$value";
return $value;
}
/*
* Computes constant pi.
*
* Implementation uses Spigot algorithms. Should converge rapidly.
* Verified to 1000 decimal places.
*
* @author Andrew Que.
* @license Public domain.
*/
function bcpi( $desiredScale = 0 )
{
// Cache value of Pi so it is not recalculated if the scale isn't changed.
static $pi = NAN;
static $piScale = 0;
// Turn up the scale for the duration of this function to avoid rounding
// error.
$scale = bcGetScale();
// Has the scale increased?
if ( ( $scale > $piScale )
|| ( $desiredScale > $piScale ) )
{
$piScale = max( $desiredScale, $scale );
// $$$FUTURE - Arbitrary. Test and determine what is actually needed.
bcscale( $piScale + 5 );
$index = 0;
$newResult = 0;
$result = -1;
while ( bccomp( $newResult, $result ) )
{
$result = $newResult;
$accumulator = bcdiv( 4, ( 8 * $index + 1 ) );
$accumulator = bcsub( $accumulator, bcdiv( 2, ( 8 * $index + 4 ) ) );
$accumulator = bcsub( $accumulator, bcdiv( 1, ( 8 * $index + 5 ) ) );
$accumulator = bcsub( $accumulator, bcdiv( 1, ( 8 * $index + 6 ) ) );
$accumulator = bcmul( $accumulator, bcdiv( 1, bcpow( 16, $index ) ) );
$newResult = bcadd( $newResult, $accumulator );
$index += 1;
}
// If are result uses all the digits, cut off the last few used to account
// for rounding error--they are incorrect anyway.
if ( strlen( $result ) > $piScale )
$result = substr( $result, 0, $piScale );
bcscale( $scale );
$pi = $result;
}
return $pi;
}
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(n)( b ) is just the one of the values from the derivative table. Calculating n! in the loop is also easy, because an! = an * an-1 where a0 = 1. The same kind of trick can be used for (x – b)n. Here pn = (x - b) * pn-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.
