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