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
|
/**************************************************************************
* *
* 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 <libultra/types.h>
#include "PR/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);
}
|
