Fast and accurate sine/cosine approximation

July 18, 2007 on 2:31 pm | In Actionscript | 44 Comments

Trigonometric functions are costly operations and can slow down your application if they are extensively used. There are two reasons why: First, Math.sin() is a function, and thus needs a function call which simple eats up some time. Second, the result is computed with much more precision than you would ever need in most situations.

Most often you just want the periodic wave-like characteristics of the sine or cosine, which can be approximated in various ways. One common way of making it faster is to create a lookup-table by computing the sine at discrete steps and storing the result in an array. For example:

var sineTable:Array = [];
for (var i:int = 0; i < 90; i++)
{
    sineTable[i] = Math.sin(Math.PI/180 * i)
}

Due to the symmetry of the sine wave, it’s sufficient to compute one quadrant only (0..pi/2), and the other 3/4’s of the circle can be computed by shifting and wrapping the input value. The biggest drawback is that the values are stored at a fixed resolution and so the result is not very accurate. This can be enhanced with linear interpolation:

x = 22.5;
y = sineTable[int(x)] + (sineTable[int(x + .5)] - sineTable[int(x)]) / 2;

Much better, but yet the error exists. It also involves accessing array elements which makes the code rather slow. Another technique uses taylor series approximation:

sin(x) = x - (x^3)/3! + (x^5)/5! - (x^7)/7! + ...

Like with the lookup-table, evaluating this term is costly.

After searching for alternatives, I finally found a fantastic solution using a quadratic curve which blows everything away in terms of performance and accuracy. For a detailed derivation, please follow the link because I won’t go into it.

I did minor optimizations to figure out what AS3 likes most, and arrived at some code that can be up to 14x faster, while still being very accurate. However, you have to use it directly - do not place the code inside a function, because the additional function call sweeps out the performance gain, and you are left with an approximation that is actually slower compared to a native Math.sin() or Math.cos() call. Also note that cos(x) = sin(x + pi/2) or cos(x - pi/2) = sin(x), so computing the cosine is just of matter adding pi/2 to the input value.

Download source: fastTrig.as.

Below is a simple visualization to show you the quality of the approximation. The high precision version can replace the Math.sin() and Math.cos() calls in nearly all situations.

Low precision sine/cosine (~14x faster)

//always wrap input angle to -PI..PI
if (x < -3.14159265)
    x += 6.28318531;
else
if (x >  3.14159265)
    x -= 6.28318531;

//compute sine
if (x < 0)
    sin = 1.27323954 * x + .405284735 * x * x;
else
    sin = 1.27323954 * x - 0.405284735 * x * x;

//compute cosine: sin(x + PI/2) = cos(x)
x += 1.57079632;
if (x >  3.14159265)
    x -= 6.28318531;

if (x < 0)
    cos = 1.27323954 * x + 0.405284735 * x * x
else
    cos = 1.27323954 * x - 0.405284735 * x * x;
}

High precision sine/cosine (~8x faster)

//always wrap input angle to -PI..PI
if (x < -3.14159265)
    x += 6.28318531;
else
if (x >  3.14159265)
    x -= 6.28318531;

//compute sine
if (x < 0)
{
    sin = 1.27323954 * x + .405284735 * x * x;

    if (sin < 0)
        sin = .225 * (sin *-sin - sin) + sin;
    else
        sin = .225 * (sin * sin - sin) + sin;
}
else
{
    sin = 1.27323954 * x - 0.405284735 * x * x;

    if (sin < 0)
        sin = .225 * (sin *-sin - sin) + sin;
    else
        sin = .225 * (sin * sin - sin) + sin;
}

//compute cosine: sin(x + PI/2) = cos(x)
x += 1.57079632;
if (x >  3.14159265)
    x -= 6.28318531;

if (x < 0)
{
    cos = 1.27323954 * x + 0.405284735 * x * x;

    if (cos < 0)
        cos = .225 * (cos *-cos - cos) + cos;
    else
        cos = .225 * (cos * cos - cos) + cos;
}
else
{
    cos = 1.27323954 * x - 0.405284735 * x * x;

    if (cos < 0)
        cos = .225 * (cos *-cos - cos) + cos;
    else
        cos = .225 * (cos * cos - cos) + cos;
}

Updates

July 16, 2007 on 1:21 am | In Actionscript, Physics | 13 Comments

Motor Physics

Finally I added support for arbitrary convex polygons, now supporting contact ids for managing persistent contacts. This is needed for a feature called warm starting which improves the quality of the solver. Furthermore the complexity of the polygon intersection algorithm was reduced from [(N0 + N1)^2] to [2*N0*N1] (N = number of separation axis to be tested) making it nearly twice as fast. Adding lines and circles is pretty straight-forward and shouldn’t take much time.

Contact point calculation also has been drastically improved. Computation now adds almost no overhead because most of the information needed is already computed during the collision detection phase. After all the simulation is very stable, and scales nicely with the amount of bodies since most algorithms are O(n) or O(log(n)). Chris is meanwhile working on different constraint types while I am enhancing the API for user interaction. We hope to get a new demo finished soon.

Data Structures

Besides some minor bug fixes, I included a new HashMap class, which is a replacement to the HashTable class (it’s much faster). The HashMap supports perfect hashing (direct lookup) using dictionaries. Think of it as an associative array with the additional capability of using object instances as keys as well.
The existing HashTable on the other hand is a ‘classic’ hash table using an array and linked overflow for resolving collisions. I leave it in the package for now, because it has one advantage over the HashMap - you can directly traverse the value/key pairs using an iterator without making a copy of the complete data set first.

Actually the HashMap class is so simple that you should considering using a Dictionary directly when speed is crucial, but it’s nice to have a common interface for all structures.