projects:devry:sinf.c [Ilive4Unix]

/*                          sin.c
 * Version History:
 *   2006Jun14 : Travis Sidelinger : Copied from uC-Linux and modified for DEVMS
 *
 * License:
 *   This program is free software; you can redistribute it and/or
 *   modify it under the terms of the GNU General Public License
 *   as published by the Free Software Foundation; either version 2
 *   of the License, or (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 *
 *  Circular sine
 *
 * SYNOPSIS:
 *
 * float x, y, sin();
 *
 * y = sin( x );
 *
 * DESCRIPTION:
 *
 * Range reduction is into intervals of pi/4.  The reduction
 * error is nearly eliminated by contriving an extended precision
 * modular arithmetic.
 *
 * Two polynomial approximating functions are employed.
 * Between 0 and pi/4 the sine is approximated by
 *      x  +  x**3 P(x**2).
 * Between pi/4 and pi/2 the cosine is represented as
 *      1  -  x**2 Q(x**2).
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain      # trials      peak       rms
 *    IEEE    -4096,+4096   100,000      1.2e-7     3.0e-8
 *    IEEE    -8192,+8192   100,000      3.0e-7     3.0e-8
 *
 * ERROR MESSAGES:
 *
 *   message           condition        value returned
 * sin total loss      x > 2^24              0.0
 *
 * Partial loss of accuracy begins to occur at x = 2^13
 * = 8192. Results may be meaningless for x >= 2^24
 * The routine as implemented flags a TLOSS error
 * for x >= 2^24 and returns 0.0.
 */
/*                          cos.c
 *
 *  Circular cosine
 *
 * SYNOPSIS:
 *
 * float x, y, cos();
 *
 * y = cos( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Range reduction is into intervals of pi/4.  The reduction
 * error is nearly eliminated by contriving an extended precision
 * modular arithmetic.
 *
 * Two polynomial approximating functions are employed.
 * Between 0 and pi/4 the cosine is approximated by
 *      1  -  x**2 Q(x**2).
 * Between pi/4 and pi/2 the sine is represented as
 *      x  +  x**3 P(x**2).
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain      # trials      peak         rms
 *    IEEE    -8192,+8192   100,000      3.0e-7     3.0e-8
 */
/*
Cephes Math Library Release 2.2:  June, 1992
Copyright 1985, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
 
 
/* Single precision circular sine
 * test interval: [-pi/4, +pi/4]
 * trials: 10000
 * peak relative error: 6.8e-8
 * rms relative error: 2.6e-8
 */
 
/* #include "mconf.h" */
 
 
static float FOPI = 1.27323954473516;
 
static float PIO4F;
/* Note, these constants are for a 32-bit significand: */
/*
static float DP1 =  0.7853851318359375;
static float DP2 =  1.30315311253070831298828125e-5;
static float DP3 =  3.03855025325309630e-11;
static float lossth = 65536.;
*/
 
/* These are for a 24-bit significand: */
static float DP1 = 0.78515625;
static float DP2 = 2.4187564849853515625e-4;
static float DP3 = 3.77489497744594108e-8;
static float lossth = 8192.;
static float T24M1 = 16777215.;
 
static float sincof[] = {
-1.9515295891E-4,
 8.3321608736E-3,
-1.6666654611E-1
};
static float coscof[] = {
 2.443315711809948E-005,
-1.388731625493765E-003,
 4.166664568298827E-002
};
 
float sin( float xx )
{
float *p;
float x, y, z;
register unsigned long j;
register int sign;
 
sign = 1;
x = xx;
if( xx < 0 )
    {
    sign = -1;
    x = -xx;
    }
if( x > T24M1 )
    {
    /* mtherr( "sin", TLOSS ); */
    return(0.0f);
    }
j = FOPI * x; /* integer part of x/(PI/4) */
y = j;
/* map zeros to origin */
if( j & 1 )
    {
    j += 1;
    y += 1.0f;
    }
j &= 7; /* octant modulo 360 degrees */
/* reflect in x axis */
if( j > 3)
    {
    sign = -sign;
    j -= 4;
    }
 
if( x > lossth )
    {
    /* mtherr( "sin", PLOSS ); */
    x = x - y * PIO4F;
    }
else
    {
/* Extended precision modular arithmetic */
    x = ((x - y * DP1) - y * DP2) - y * DP3;
    }
/*einits();*/
z = x * x;
if( (j==1) || (j==2) )
    {
/* measured relative error in +/- pi/4 is 7.8e-8 */
/*
    y = ((  2.443315711809948E-005f * z
      - 1.388731625493765E-003f) * z
      + 4.166664568298827E-002f) * z * z;
*/
    p = coscof;
    y = *p++;
    y = y * z + *p++;
    y = y * z + *p++;
    y *= z * z;
    y -= 0.5f * z;
    y += 1.0f;
    }
else
    {
/* Theoretical relative error = 3.8e-9 in [-pi/4, +pi/4] */
/*
    y = ((-1.9515295891E-4f * z
         + 8.3321608736E-3f) * z
         - 1.6666654611E-1f) * z * x;
    y += x;
*/
    p = sincof;
    y = *p++;
    y = y * z + *p++;
    y = y * z + *p++;
    y *= z * x;
    y += x;
    }
/*einitd();*/
if(sign < 0)
    y = -y;
return( y);
}
 
 
/* Single precision circular cosine
 * test interval: [-pi/4, +pi/4]
 * trials: 10000
 * peak relative error: 8.3e-8
 * rms relative error: 2.2e-8
 */
 
float cos( float xx )
{
float x, y, z;
int j, sign;
 
/* make argument positive */
sign = 1;
x = xx;
if( x < 0 )
    x = -x;
 
if( x > T24M1 )
    {
    /* mtherr( "cos", TLOSS ); */
    return(0.0f);
    }
 
j = FOPI * x; /* integer part of x/PIO4 */
y = j;
/* integer and fractional part modulo one octant */
if( j & 1 ) /* map zeros to origin */
    {
    j += 1;
    y += 1.0f;
    }
j &= 7;
if( j > 3)
    {
    j -=4;
    sign = -sign;
    }
 
if( j > 1 )
    sign = -sign;
 
if( x > lossth )
    {
        /* mtherr( "cos", PLOSS ); */
    x = x - y * PIO4F;
    }
else
/* Extended precision modular arithmetic */
    x = ((x - y * DP1) - y * DP2) - y * DP3;
 
z = x * x;
 
if( (j==1) || (j==2) )
    {
    y = (((-1.9515295891E-4f * z
         + 8.3321608736E-3f) * z
         - 1.6666654611E-1f) * z * x)
         + x;
    }
else
    {
    y = ((  2.443315711809948E-005f * z
      - 1.388731625493765E-003f) * z
      + 4.166664568298827E-002f) * z * z;
    y -= 0.5f * z;
    y += 1.0f;
    }
if(sign < 0)
    y = -y;
return( y );
}