What is it really?

I would like to think that I am a pretty decent person looking at the code of other peoples ... but I'm at a loss. This comes from the Doom3 Math library. I believe this probably existed in the GPL since Quake 1. Note that this makes reference to math.h. I guess there is a way that this actually calculates the cosine ... but I can't figure it out. Does anyone explain?

    ID_INLINE float idMath::Cos16( float a ) {
    float s, d;

    if ( ( a < 0.0f ) || ( a >= TWO_PI ) ) {
        a -= floorf( a / TWO_PI ) * TWO_PI;
    }
#if 1
    if ( a < PI ) {
        if ( a > HALF_PI ) {
            a = PI - a;
            d = -1.0f;
        } else {
            d = 1.0f;
        }
    } else {
        if ( a > PI + HALF_PI ) {
            a = a - TWO_PI;
            d = 1.0f;
        } else {
            a = PI - a;
            d = -1.0f;
        }
    }
#else
    a = PI - a;
    if ( fabs( a ) >= HALF_PI ) {
        a = ( ( a < 0.0f ) ? -PI : PI ) - a;
        d = 1.0f;
    } else {
        d = -1.0f;
    }
#endif
    s = a * a;
    return d * ( ( ( ( ( -2.605e-07f * s + 2.47609e-05f ) * s - 1.3888397e-03f ) * s + 4.16666418e-02f ) * s - 4.999999963e-01f ) * s + 1.0f );
}

      

+3


source to share


1 answer


All conventional things look like they are rolling the angle down to one quadrant (or maybe an octant, I can't figure out how to figure it out ...;)) and writing down the correction factor ( d

).

The last line does some 10th order polynomial approximation (maybe Taylor or Chebyshev ?). * But it only works on even powers, and cos

is an even function. He also uses Horner's method to avoid calculating high powers directly multiple times.

Then reapply the correct sign with d

.




* To make sure this is "perfectly" accurate, try running the following code in Octave (you can do it online, for example here ):
% One quadrant
a = (-pi/2):0.1:(+pi/2);

% Exact result
y_exact = cos(a);

% Our approximation
s = a .* a;
y_approx = (((((-2.605e-07 .* s + 2.47609e-05) .* s - 1.3888397e-03) .* s + 4.16666418e-02) .* s - 4.999999963e-01) .* s + 1);

% Plot
hold on
plot(a, y_exact,  'b')
plot(a, y_approx, 'r')

      

sub>

+5


source







All Articles