1 /*---------------------------------------------------------------------------+
2 | poly_sin.c |
3 | |
4 | Computation of an approximation of the sin function and the cosine |
5 | function by a polynomial. |
6 | |
7 | Copyright (C) 1992,1993,1994,1997,1999 |
8 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
9 | E-mail billm@melbpc.org.au |
10 | |
11 | |
12 +---------------------------------------------------------------------------*/
13
14 #include "exception.h"
15 #include "reg_constant.h"
16 #include "fpu_emu.h"
17 #include "fpu_system.h"
18 #include "control_w.h"
19 #include "poly.h"
20
21 #define N_COEFF_P 4
22 #define N_COEFF_N 4
23
24 static const unsigned long long pos_terms_l[N_COEFF_P] = {
25 0xaaaaaaaaaaaaaaabLL,
26 0x00d00d00d00cf906LL,
27 0x000006b99159a8bbLL,
28 0x000000000d7392e6LL
29 };
30
31 static const unsigned long long neg_terms_l[N_COEFF_N] = {
32 0x2222222222222167LL,
33 0x0002e3bc74aab624LL,
34 0x0000000b09229062LL,
35 0x00000000000c7973LL
36 };
37
38 #define N_COEFF_PH 4
39 #define N_COEFF_NH 4
40 static const unsigned long long pos_terms_h[N_COEFF_PH] = {
41 0x0000000000000000LL,
42 0x05b05b05b05b0406LL,
43 0x000049f93edd91a9LL,
44 0x00000000c9c9ed62LL
45 };
46
47 static const unsigned long long neg_terms_h[N_COEFF_NH] = {
48 0xaaaaaaaaaaaaaa98LL,
49 0x001a01a01a019064LL,
50 0x0000008f76c68a77LL,
51 0x0000000000d58f5eLL
52 };
53
54 /*--- poly_sine() -----------------------------------------------------------+
55 | |
56 +---------------------------------------------------------------------------*/
poly_sine(FPU_REG * st0_ptr)57 void poly_sine(FPU_REG *st0_ptr)
58 {
59 int exponent, echange;
60 Xsig accumulator, argSqrd, argTo4;
61 unsigned long fix_up, adj;
62 unsigned long long fixed_arg;
63 FPU_REG result;
64
65 exponent = exponent(st0_ptr);
66
67 accumulator.lsw = accumulator.midw = accumulator.msw = 0;
68
69 /* Split into two ranges, for arguments below and above 1.0 */
70 /* The boundary between upper and lower is approx 0.88309101259 */
71 if ((exponent < -1)
72 || ((exponent == -1) && (st0_ptr->sigh <= 0xe21240aa))) {
73 /* The argument is <= 0.88309101259 */
74
75 argSqrd.msw = st0_ptr->sigh;
76 argSqrd.midw = st0_ptr->sigl;
77 argSqrd.lsw = 0;
78 mul64_Xsig(&argSqrd, &significand(st0_ptr));
79 shr_Xsig(&argSqrd, 2 * (-1 - exponent));
80 argTo4.msw = argSqrd.msw;
81 argTo4.midw = argSqrd.midw;
82 argTo4.lsw = argSqrd.lsw;
83 mul_Xsig_Xsig(&argTo4, &argTo4);
84
85 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
86 N_COEFF_N - 1);
87 mul_Xsig_Xsig(&accumulator, &argSqrd);
88 negate_Xsig(&accumulator);
89
90 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
91 N_COEFF_P - 1);
92
93 shr_Xsig(&accumulator, 2); /* Divide by four */
94 accumulator.msw |= 0x80000000; /* Add 1.0 */
95
96 mul64_Xsig(&accumulator, &significand(st0_ptr));
97 mul64_Xsig(&accumulator, &significand(st0_ptr));
98 mul64_Xsig(&accumulator, &significand(st0_ptr));
99
100 /* Divide by four, FPU_REG compatible, etc */
101 exponent = 3 * exponent;
102
103 /* The minimum exponent difference is 3 */
104 shr_Xsig(&accumulator, exponent(st0_ptr) - exponent);
105
106 negate_Xsig(&accumulator);
107 XSIG_LL(accumulator) += significand(st0_ptr);
108
109 echange = round_Xsig(&accumulator);
110
111 setexponentpos(&result, exponent(st0_ptr) + echange);
112 } else {
113 /* The argument is > 0.88309101259 */
114 /* We use sin(st(0)) = cos(pi/2-st(0)) */
115
116 fixed_arg = significand(st0_ptr);
117
118 if (exponent == 0) {
119 /* The argument is >= 1.0 */
120
121 /* Put the binary point at the left. */
122 fixed_arg <<= 1;
123 }
124 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
125 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
126 /* There is a special case which arises due to rounding, to fix here. */
127 if (fixed_arg == 0xffffffffffffffffLL)
128 fixed_arg = 0;
129
130 XSIG_LL(argSqrd) = fixed_arg;
131 argSqrd.lsw = 0;
132 mul64_Xsig(&argSqrd, &fixed_arg);
133
134 XSIG_LL(argTo4) = XSIG_LL(argSqrd);
135 argTo4.lsw = argSqrd.lsw;
136 mul_Xsig_Xsig(&argTo4, &argTo4);
137
138 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
139 N_COEFF_NH - 1);
140 mul_Xsig_Xsig(&accumulator, &argSqrd);
141 negate_Xsig(&accumulator);
142
143 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
144 N_COEFF_PH - 1);
145 negate_Xsig(&accumulator);
146
147 mul64_Xsig(&accumulator, &fixed_arg);
148 mul64_Xsig(&accumulator, &fixed_arg);
149
150 shr_Xsig(&accumulator, 3);
151 negate_Xsig(&accumulator);
152
153 add_Xsig_Xsig(&accumulator, &argSqrd);
154
155 shr_Xsig(&accumulator, 1);
156
157 accumulator.lsw |= 1; /* A zero accumulator here would cause problems */
158 negate_Xsig(&accumulator);
159
160 /* The basic computation is complete. Now fix the answer to
161 compensate for the error due to the approximation used for
162 pi/2
163 */
164
165 /* This has an exponent of -65 */
166 fix_up = 0x898cc517;
167 /* The fix-up needs to be improved for larger args */
168 if (argSqrd.msw & 0xffc00000) {
169 /* Get about 32 bit precision in these: */
170 fix_up -= mul_32_32(0x898cc517, argSqrd.msw) / 6;
171 }
172 fix_up = mul_32_32(fix_up, LL_MSW(fixed_arg));
173
174 adj = accumulator.lsw; /* temp save */
175 accumulator.lsw -= fix_up;
176 if (accumulator.lsw > adj)
177 XSIG_LL(accumulator)--;
178
179 echange = round_Xsig(&accumulator);
180
181 setexponentpos(&result, echange - 1);
182 }
183
184 significand(&result) = XSIG_LL(accumulator);
185 setsign(&result, getsign(st0_ptr));
186 FPU_copy_to_reg0(&result, TAG_Valid);
187
188 #ifdef PARANOID
189 if ((exponent(&result) >= 0)
190 && (significand(&result) > 0x8000000000000000LL)) {
191 EXCEPTION(EX_INTERNAL | 0x150);
192 }
193 #endif /* PARANOID */
194
195 }
196
197 /*--- poly_cos() ------------------------------------------------------------+
198 | |
199 +---------------------------------------------------------------------------*/
poly_cos(FPU_REG * st0_ptr)200 void poly_cos(FPU_REG *st0_ptr)
201 {
202 FPU_REG result;
203 long int exponent, exp2, echange;
204 Xsig accumulator, argSqrd, fix_up, argTo4;
205 unsigned long long fixed_arg;
206
207 #ifdef PARANOID
208 if ((exponent(st0_ptr) > 0)
209 || ((exponent(st0_ptr) == 0)
210 && (significand(st0_ptr) > 0xc90fdaa22168c234LL))) {
211 EXCEPTION(EX_Invalid);
212 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
213 return;
214 }
215 #endif /* PARANOID */
216
217 exponent = exponent(st0_ptr);
218
219 accumulator.lsw = accumulator.midw = accumulator.msw = 0;
220
221 if ((exponent < -1)
222 || ((exponent == -1) && (st0_ptr->sigh <= 0xb00d6f54))) {
223 /* arg is < 0.687705 */
224
225 argSqrd.msw = st0_ptr->sigh;
226 argSqrd.midw = st0_ptr->sigl;
227 argSqrd.lsw = 0;
228 mul64_Xsig(&argSqrd, &significand(st0_ptr));
229
230 if (exponent < -1) {
231 /* shift the argument right by the required places */
232 shr_Xsig(&argSqrd, 2 * (-1 - exponent));
233 }
234
235 argTo4.msw = argSqrd.msw;
236 argTo4.midw = argSqrd.midw;
237 argTo4.lsw = argSqrd.lsw;
238 mul_Xsig_Xsig(&argTo4, &argTo4);
239
240 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
241 N_COEFF_NH - 1);
242 mul_Xsig_Xsig(&accumulator, &argSqrd);
243 negate_Xsig(&accumulator);
244
245 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
246 N_COEFF_PH - 1);
247 negate_Xsig(&accumulator);
248
249 mul64_Xsig(&accumulator, &significand(st0_ptr));
250 mul64_Xsig(&accumulator, &significand(st0_ptr));
251 shr_Xsig(&accumulator, -2 * (1 + exponent));
252
253 shr_Xsig(&accumulator, 3);
254 negate_Xsig(&accumulator);
255
256 add_Xsig_Xsig(&accumulator, &argSqrd);
257
258 shr_Xsig(&accumulator, 1);
259
260 /* It doesn't matter if accumulator is all zero here, the
261 following code will work ok */
262 negate_Xsig(&accumulator);
263
264 if (accumulator.lsw & 0x80000000)
265 XSIG_LL(accumulator)++;
266 if (accumulator.msw == 0) {
267 /* The result is 1.0 */
268 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
269 return;
270 } else {
271 significand(&result) = XSIG_LL(accumulator);
272
273 /* will be a valid positive nr with expon = -1 */
274 setexponentpos(&result, -1);
275 }
276 } else {
277 fixed_arg = significand(st0_ptr);
278
279 if (exponent == 0) {
280 /* The argument is >= 1.0 */
281
282 /* Put the binary point at the left. */
283 fixed_arg <<= 1;
284 }
285 /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
286 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
287 /* There is a special case which arises due to rounding, to fix here. */
288 if (fixed_arg == 0xffffffffffffffffLL)
289 fixed_arg = 0;
290
291 exponent = -1;
292 exp2 = -1;
293
294 /* A shift is needed here only for a narrow range of arguments,
295 i.e. for fixed_arg approx 2^-32, but we pick up more... */
296 if (!(LL_MSW(fixed_arg) & 0xffff0000)) {
297 fixed_arg <<= 16;
298 exponent -= 16;
299 exp2 -= 16;
300 }
301
302 XSIG_LL(argSqrd) = fixed_arg;
303 argSqrd.lsw = 0;
304 mul64_Xsig(&argSqrd, &fixed_arg);
305
306 if (exponent < -1) {
307 /* shift the argument right by the required places */
308 shr_Xsig(&argSqrd, 2 * (-1 - exponent));
309 }
310
311 argTo4.msw = argSqrd.msw;
312 argTo4.midw = argSqrd.midw;
313 argTo4.lsw = argSqrd.lsw;
314 mul_Xsig_Xsig(&argTo4, &argTo4);
315
316 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
317 N_COEFF_N - 1);
318 mul_Xsig_Xsig(&accumulator, &argSqrd);
319 negate_Xsig(&accumulator);
320
321 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
322 N_COEFF_P - 1);
323
324 shr_Xsig(&accumulator, 2); /* Divide by four */
325 accumulator.msw |= 0x80000000; /* Add 1.0 */
326
327 mul64_Xsig(&accumulator, &fixed_arg);
328 mul64_Xsig(&accumulator, &fixed_arg);
329 mul64_Xsig(&accumulator, &fixed_arg);
330
331 /* Divide by four, FPU_REG compatible, etc */
332 exponent = 3 * exponent;
333
334 /* The minimum exponent difference is 3 */
335 shr_Xsig(&accumulator, exp2 - exponent);
336
337 negate_Xsig(&accumulator);
338 XSIG_LL(accumulator) += fixed_arg;
339
340 /* The basic computation is complete. Now fix the answer to
341 compensate for the error due to the approximation used for
342 pi/2
343 */
344
345 /* This has an exponent of -65 */
346 XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
347 fix_up.lsw = 0;
348
349 /* The fix-up needs to be improved for larger args */
350 if (argSqrd.msw & 0xffc00000) {
351 /* Get about 32 bit precision in these: */
352 fix_up.msw -= mul_32_32(0x898cc517, argSqrd.msw) / 2;
353 fix_up.msw += mul_32_32(0x898cc517, argTo4.msw) / 24;
354 }
355
356 exp2 += norm_Xsig(&accumulator);
357 shr_Xsig(&accumulator, 1); /* Prevent overflow */
358 exp2++;
359 shr_Xsig(&fix_up, 65 + exp2);
360
361 add_Xsig_Xsig(&accumulator, &fix_up);
362
363 echange = round_Xsig(&accumulator);
364
365 setexponentpos(&result, exp2 + echange);
366 significand(&result) = XSIG_LL(accumulator);
367 }
368
369 FPU_copy_to_reg0(&result, TAG_Valid);
370
371 #ifdef PARANOID
372 if ((exponent(&result) >= 0)
373 && (significand(&result) > 0x8000000000000000LL)) {
374 EXCEPTION(EX_INTERNAL | 0x151);
375 }
376 #endif /* PARANOID */
377
378 }
379