138 lines
2.8 KiB
C
138 lines
2.8 KiB
C
/**************************************************************************
|
||
* *
|
||
* Copyright (C) 1994, Silicon Graphics, Inc. *
|
||
* *
|
||
* These coded instructions, statements, and computer programs contain *
|
||
* unpublished proprietary information of Silicon Graphics, Inc., and *
|
||
* are protected by Federal copyright law. They may not be disclosed *
|
||
* to third parties or copied or duplicated in any form, in whole or *
|
||
* in part, without the prior written consent of Silicon Graphics, Inc. *
|
||
* *
|
||
**************************************************************************/
|
||
|
||
#include "guint.h"
|
||
|
||
/* ====================================================================
|
||
* ====================================================================
|
||
*
|
||
* Module: fcos.c
|
||
* $Revision: 1.2 $
|
||
* $Date: 1995/07/12 17:47:57 $
|
||
* $Author: jeffd $
|
||
* $Source: /disk6/Master/cvsmdev2/PR/libultra/gu/cosf.c,v $
|
||
*
|
||
* Revision history:
|
||
* 09-Jun-93 - Original Version
|
||
*
|
||
* Description: source code for fcos function
|
||
*
|
||
* ====================================================================
|
||
* ====================================================================
|
||
*/
|
||
|
||
#pragma weak fcos = __cosf
|
||
#pragma weak cosf = __cosf
|
||
#define fcos __cosf
|
||
|
||
/* coefficients for polynomial approximation of cos on +/- pi/2 */
|
||
|
||
static const du P[] =
|
||
{
|
||
{0x3ff00000, 0x00000000},
|
||
{0xbfc55554, 0xbc83656d},
|
||
{0x3f8110ed, 0x3804c2a0},
|
||
{0xbf29f6ff, 0xeea56814},
|
||
{0x3ec5dbdf, 0x0e314bfe},
|
||
};
|
||
|
||
static const du rpi =
|
||
{0x3fd45f30, 0x6dc9c883};
|
||
|
||
static const du pihi =
|
||
{0x400921fb, 0x50000000};
|
||
|
||
static const du pilo =
|
||
{0x3e6110b4, 0x611a6263};
|
||
|
||
static const fu zero = {0x00000000};
|
||
|
||
|
||
/* ====================================================================
|
||
*
|
||
* FunctionName fcos
|
||
*
|
||
* Description computes cosine of arg
|
||
*
|
||
* ====================================================================
|
||
*/
|
||
|
||
float
|
||
fcos( float x )
|
||
{
|
||
float absx;
|
||
double dx, xsq, poly;
|
||
double dn;
|
||
int n;
|
||
double result;
|
||
int ix, xpt;
|
||
|
||
|
||
ix = *(int *)&x;
|
||
xpt = (ix >> 22);
|
||
xpt &= 0x1ff;
|
||
|
||
/* xpt is exponent(x) + 1 bit of mantissa */
|
||
|
||
|
||
if ( xpt < 0x136 )
|
||
{
|
||
/* |x| < 2^28 */
|
||
|
||
/* use the standard algorithm from Cody and Waite, doing
|
||
the computations in double precision
|
||
*/
|
||
|
||
absx = ABS(x);
|
||
|
||
dx = absx;
|
||
|
||
dn = dx*rpi.d + 0.5;
|
||
n = ROUND(dn);
|
||
dn = n;
|
||
|
||
dn -= 0.5;
|
||
|
||
dx = dx - dn*pihi.d;
|
||
dx = dx - dn*pilo.d; /* dx = x - (n - 0.5)*pi */
|
||
|
||
xsq = dx*dx;
|
||
|
||
poly = ((P[4].d*xsq + P[3].d)*xsq + P[2].d)*xsq + P[1].d;
|
||
|
||
result = dx + (dx*xsq)*poly;
|
||
|
||
/* negate result if n is odd */
|
||
|
||
if ( (n & 1) == 0 )
|
||
return ( (float)result );
|
||
|
||
return ( -(float)result );
|
||
}
|
||
|
||
if ( x != x )
|
||
{
|
||
/* x is a NaN; return a quiet NaN */
|
||
|
||
#ifdef _IP_NAN_SETS_ERRNO
|
||
|
||
*__errnoaddr = EDOM;
|
||
#endif
|
||
|
||
return ( __libm_qnan_f );
|
||
}
|
||
|
||
/* just give up and return 0.0 */
|
||
|
||
return ( zero.f );
|
||
}
|