This source file includes following definitions.
- poly_sine
1
2
3
4
5
6
7
8
9
10
11
12
13 #include "exception.h"
14 #include "reg_constant.h"
15 #include "fpu_emu.h"
16
17
18 #define HIPOWER 5
19 static unsigned short lterms[HIPOWER][4] =
20 {
21 { 0x846a, 0x42d1, 0xb544, 0x921f},
22 { 0xe110, 0x75aa, 0xbc67, 0x1466},
23 { 0x503d, 0xa43f, 0x83c1, 0x000a},
24 { 0x8f9d, 0x7a19, 0x00f4, 0x0000},
25 { 0xda03, 0x06aa, 0x0000, 0x0000},
26 };
27
28 static unsigned short negterms[HIPOWER][4] =
29 {
30 { 0x95ed, 0x2df2, 0xe731, 0xa55d},
31 { 0xd159, 0xe62b, 0xd2cc, 0x0132},
32 { 0x6342, 0xe9fb, 0x3c60, 0x0000},
33 { 0x6256, 0xdf5a, 0x0002, 0x0000},
34 { 0xf279, 0x000b, 0x0000, 0x0000},
35 };
36
37
38
39
40
41 void poly_sine(FPU_REG *arg, FPU_REG *result)
42 {
43 short exponent;
44 FPU_REG Xx, Xx2, Xx4, accum, negaccum;
45
46
47 exponent = arg->exp - EXP_BIAS;
48
49 if ( arg->tag == TW_Zero )
50 {
51
52 reg_move(&CONST_Z, result);
53 return;
54 }
55
56 #ifdef PARANOID
57 if ( arg->sign != 0 )
58 {
59 EXCEPTION(EX_Invalid);
60 reg_move(&CONST_QNaN, result);
61 return;
62 }
63
64 if ( exponent >= 0 )
65 {
66 if ( (exponent == 0) && (arg->sigl == 0) && (arg->sigh == 0x80000000) )
67 {
68 reg_move(&CONST_1, result);
69 return;
70 }
71 EXCEPTION(EX_Invalid);
72 reg_move(&CONST_QNaN, result);
73 return;
74 }
75 #endif PARANOID
76
77 Xx.sigl = arg->sigl;
78 Xx.sigh = arg->sigh;
79 if ( exponent < -1 )
80 {
81
82 if ( shrx(&(Xx.sigl), -1-exponent) >= 0x80000000U )
83 (*((long long *)(&(Xx.sigl))))++;
84 }
85
86 mul64((long long *)&(Xx.sigl), (long long *)&(Xx.sigl),
87 (long long *)&(Xx2.sigl));
88 mul64((long long *)&(Xx2.sigl), (long long *)&(Xx2.sigl),
89 (long long *)&(Xx4.sigl));
90
91
92 *(short *)&(accum.sign) = 0;
93 accum.exp = 0;
94
95
96 polynomial(&(accum.sigl), &(Xx4.sigl), lterms, HIPOWER-1);
97
98
99 *(short *)&(negaccum.sign) = 0;
100 negaccum.exp = 0;
101
102
103 polynomial(&(negaccum.sigl), &(Xx4.sigl), negterms, HIPOWER-1);
104 mul64((long long *)&(Xx2.sigl), (long long *)&(negaccum.sigl),
105 (long long *)&(negaccum.sigl));
106
107
108 *((long long *)(&(accum.sigl))) -= *((long long *)(&(negaccum.sigl)));
109
110
111 accum.exp = EXP_BIAS - 1 + accum.exp;
112
113 *(short *)&(result->sign) = *(short *)&(accum.sign);
114 result->exp = accum.exp;
115 result->sigl = accum.sigl;
116 result->sigh = accum.sigh;
117
118 normalize(result);
119
120 reg_mul(result, arg, result);
121 reg_u_add(result, arg, result);
122
123
124 if ( result->exp >= EXP_BIAS )
125 {
126 if ( (result->exp > EXP_BIAS)
127 || (result->sigl > 1)
128 || (result->sigh != 0x80000000)
129 )
130 {
131 #ifdef DEBUGGING
132 RE_ENTRANT_CHECK_OFF
133 printk("\nEXP=%d, MS=%08x, LS=%08x\n", result->exp,
134 result->sigh, result->sigl);
135 RE_ENTRANT_CHECK_ON
136 #endif DEBUGGING
137 EXCEPTION(EX_INTERNAL|0x103);
138 }
139
140 #ifdef DEBUGGING
141 RE_ENTRANT_CHECK_OFF
142 printk("\n***CORRECTING ILLEGAL RESULT*** in poly_sin() computation\n");
143 printk("EXP=%d, MS=%08x, LS=%08x\n", result->exp,
144 result->sigh, result->sigl);
145 RE_ENTRANT_CHECK_ON
146 #endif DEBUGGING
147
148 result->sigl = 0;
149 }
150 }