1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
/**************************************************************************
* *
* 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.3 $
* $Date: 1998/10/09 06:10:53 $
* $Author: has $
* $Source: /exdisk2/cvs/N64OS/Master/cvsmdev2/PR/libultra/monegi/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 );
}
|
