This source file includes following definitions.
- poly_sine
- poly_cos
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 #include "exception.h"
16 #include "reg_constant.h"
17 #include "fpu_emu.h"
18 #include "control_w.h"
19 #include "poly.h"
20
21
22 #define N_COEFF_P 4
23 #define N_COEFF_N 4
24
25 static const unsigned long long pos_terms_l[N_COEFF_P] =
26 {
27 0xaaaaaaaaaaaaaaabLL,
28 0x00d00d00d00cf906LL,
29 0x000006b99159a8bbLL,
30 0x000000000d7392e6LL
31 };
32
33 static const unsigned long long neg_terms_l[N_COEFF_N] =
34 {
35 0x2222222222222167LL,
36 0x0002e3bc74aab624LL,
37 0x0000000b09229062LL,
38 0x00000000000c7973LL
39 };
40
41
42
43 #define N_COEFF_PH 4
44 #define N_COEFF_NH 4
45 static const unsigned long long pos_terms_h[N_COEFF_PH] =
46 {
47 0x0000000000000000LL,
48 0x05b05b05b05b0406LL,
49 0x000049f93edd91a9LL,
50 0x00000000c9c9ed62LL
51 };
52
53 static const unsigned long long neg_terms_h[N_COEFF_NH] =
54 {
55 0xaaaaaaaaaaaaaa98LL,
56 0x001a01a01a019064LL,
57 0x0000008f76c68a77LL,
58 0x0000000000d58f5eLL
59 };
60
61
62
63
64
65 void poly_sine(FPU_REG const *arg, FPU_REG *result)
66 {
67 int exponent, echange;
68 Xsig accumulator, argSqrd, argTo4;
69 unsigned long fix_up, adj;
70 unsigned long long fixed_arg;
71
72
73 #ifdef PARANOID
74 if ( arg->tag == TW_Zero )
75 {
76
77 reg_move(&CONST_Z, result);
78 return;
79 }
80 #endif PARANOID
81
82 exponent = arg->exp - EXP_BIAS;
83
84 accumulator.lsw = accumulator.midw = accumulator.msw = 0;
85
86
87
88 if ( (exponent < -1) || ((exponent == -1) && (arg->sigh <= 0xe21240aa)) )
89 {
90
91
92 argSqrd.msw = arg->sigh; argSqrd.midw = arg->sigl; argSqrd.lsw = 0;
93 mul64_Xsig(&argSqrd, &significand(arg));
94 shr_Xsig(&argSqrd, 2*(-1-exponent));
95 argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
96 argTo4.lsw = argSqrd.lsw;
97 mul_Xsig_Xsig(&argTo4, &argTo4);
98
99 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
100 N_COEFF_N-1);
101 mul_Xsig_Xsig(&accumulator, &argSqrd);
102 negate_Xsig(&accumulator);
103
104 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
105 N_COEFF_P-1);
106
107 shr_Xsig(&accumulator, 2);
108 accumulator.msw |= 0x80000000;
109
110 mul64_Xsig(&accumulator, &significand(arg));
111 mul64_Xsig(&accumulator, &significand(arg));
112 mul64_Xsig(&accumulator, &significand(arg));
113
114
115 exponent = 3*exponent + EXP_BIAS;
116
117
118 shr_Xsig(&accumulator, arg->exp - exponent);
119
120 negate_Xsig(&accumulator);
121 XSIG_LL(accumulator) += significand(arg);
122
123 echange = round_Xsig(&accumulator);
124
125 result->exp = arg->exp + echange;
126 }
127 else
128 {
129
130
131
132 fixed_arg = significand(arg);
133
134 if ( exponent == 0 )
135 {
136
137
138
139 fixed_arg <<= 1;
140 }
141
142 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
143
144 XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
145 mul64_Xsig(&argSqrd, &fixed_arg);
146
147 XSIG_LL(argTo4) = XSIG_LL(argSqrd); argTo4.lsw = argSqrd.lsw;
148 mul_Xsig_Xsig(&argTo4, &argTo4);
149
150 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
151 N_COEFF_NH-1);
152 mul_Xsig_Xsig(&accumulator, &argSqrd);
153 negate_Xsig(&accumulator);
154
155 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
156 N_COEFF_PH-1);
157 negate_Xsig(&accumulator);
158
159 mul64_Xsig(&accumulator, &fixed_arg);
160 mul64_Xsig(&accumulator, &fixed_arg);
161
162 shr_Xsig(&accumulator, 3);
163 negate_Xsig(&accumulator);
164
165 add_Xsig_Xsig(&accumulator, &argSqrd);
166
167 shr_Xsig(&accumulator, 1);
168
169 accumulator.lsw |= 1;
170 negate_Xsig(&accumulator);
171
172
173
174
175
176
177
178 fix_up = 0x898cc517;
179
180 if ( argSqrd.msw & 0xffc00000 )
181 {
182
183 mul_32_32(0x898cc517, argSqrd.msw, &adj);
184 fix_up -= adj/6;
185 }
186 mul_32_32(fix_up, LL_MSW(fixed_arg), &fix_up);
187
188 adj = accumulator.lsw;
189 accumulator.lsw -= fix_up;
190 if ( accumulator.lsw > adj )
191 XSIG_LL(accumulator) --;
192
193 echange = round_Xsig(&accumulator);
194
195 result->exp = EXP_BIAS - 1 + echange;
196 }
197
198 significand(result) = XSIG_LL(accumulator);
199 result->tag = TW_Valid;
200 result->sign = arg->sign;
201
202 #ifdef PARANOID
203 if ( (result->exp >= EXP_BIAS)
204 && (significand(result) > 0x8000000000000000LL) )
205 {
206 EXCEPTION(EX_INTERNAL|0x150);
207 }
208 #endif PARANOID
209
210 }
211
212
213
214
215
216
217 void poly_cos(FPU_REG const *arg, FPU_REG *result)
218 {
219 long int exponent, exp2, echange;
220 Xsig accumulator, argSqrd, fix_up, argTo4;
221 unsigned long adj;
222 unsigned long long fixed_arg;
223
224
225 #ifdef PARANOID
226 if ( arg->tag == TW_Zero )
227 {
228
229 reg_move(&CONST_1, result);
230 return;
231 }
232
233 if ( (arg->exp > EXP_BIAS)
234 || ((arg->exp == EXP_BIAS)
235 && (significand(arg) > 0xc90fdaa22168c234LL)) )
236 {
237 EXCEPTION(EX_Invalid);
238 reg_move(&CONST_QNaN, result);
239 return;
240 }
241 #endif PARANOID
242
243 exponent = arg->exp - EXP_BIAS;
244
245 accumulator.lsw = accumulator.midw = accumulator.msw = 0;
246
247 if ( (exponent < -1) || ((exponent == -1) && (arg->sigh <= 0xb00d6f54)) )
248 {
249
250
251 argSqrd.msw = arg->sigh; argSqrd.midw = arg->sigl; argSqrd.lsw = 0;
252 mul64_Xsig(&argSqrd, &significand(arg));
253
254 if ( exponent < -1 )
255 {
256
257 shr_Xsig(&argSqrd, 2*(-1-exponent));
258 }
259
260 argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
261 argTo4.lsw = argSqrd.lsw;
262 mul_Xsig_Xsig(&argTo4, &argTo4);
263
264 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_h,
265 N_COEFF_NH-1);
266 mul_Xsig_Xsig(&accumulator, &argSqrd);
267 negate_Xsig(&accumulator);
268
269 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_h,
270 N_COEFF_PH-1);
271 negate_Xsig(&accumulator);
272
273 mul64_Xsig(&accumulator, &significand(arg));
274 mul64_Xsig(&accumulator, &significand(arg));
275 shr_Xsig(&accumulator, -2*(1+exponent));
276
277 shr_Xsig(&accumulator, 3);
278 negate_Xsig(&accumulator);
279
280 add_Xsig_Xsig(&accumulator, &argSqrd);
281
282 shr_Xsig(&accumulator, 1);
283
284
285
286 negate_Xsig(&accumulator);
287
288 if ( accumulator.lsw & 0x80000000 )
289 XSIG_LL(accumulator) ++;
290 if ( accumulator.msw == 0 )
291 {
292
293 reg_move(&CONST_1, result);
294 }
295 else
296 {
297 significand(result) = XSIG_LL(accumulator);
298
299
300 *(short *)&(result->sign) = 0;
301 result->exp = EXP_BIAS - 1;
302 }
303 }
304 else
305 {
306 fixed_arg = significand(arg);
307
308 if ( exponent == 0 )
309 {
310
311
312
313 fixed_arg <<= 1;
314 }
315
316 fixed_arg = 0x921fb54442d18469LL - fixed_arg;
317
318 exponent = -1;
319 exp2 = -1;
320
321
322
323 if ( !(LL_MSW(fixed_arg) & 0xffff0000) )
324 {
325 fixed_arg <<= 16;
326 exponent -= 16;
327 exp2 -= 16;
328 }
329
330 XSIG_LL(argSqrd) = fixed_arg; argSqrd.lsw = 0;
331 mul64_Xsig(&argSqrd, &fixed_arg);
332
333 if ( exponent < -1 )
334 {
335
336 shr_Xsig(&argSqrd, 2*(-1-exponent));
337 }
338
339 argTo4.msw = argSqrd.msw; argTo4.midw = argSqrd.midw;
340 argTo4.lsw = argSqrd.lsw;
341 mul_Xsig_Xsig(&argTo4, &argTo4);
342
343 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), neg_terms_l,
344 N_COEFF_N-1);
345 mul_Xsig_Xsig(&accumulator, &argSqrd);
346 negate_Xsig(&accumulator);
347
348 polynomial_Xsig(&accumulator, &XSIG_LL(argTo4), pos_terms_l,
349 N_COEFF_P-1);
350
351 shr_Xsig(&accumulator, 2);
352 accumulator.msw |= 0x80000000;
353
354 mul64_Xsig(&accumulator, &fixed_arg);
355 mul64_Xsig(&accumulator, &fixed_arg);
356 mul64_Xsig(&accumulator, &fixed_arg);
357
358
359 exponent = 3*exponent;
360
361
362 shr_Xsig(&accumulator, exp2 - exponent);
363
364 negate_Xsig(&accumulator);
365 XSIG_LL(accumulator) += fixed_arg;
366
367
368
369
370
371
372
373 XSIG_LL(fix_up) = 0x898cc51701b839a2ll;
374 fix_up.lsw = 0;
375
376
377 if ( argSqrd.msw & 0xffc00000 )
378 {
379
380 mul_32_32(0x898cc517, argSqrd.msw, &adj);
381 fix_up.msw -= adj/2;
382 mul_32_32(0x898cc517, argTo4.msw, &adj);
383 fix_up.msw += adj/24;
384 }
385
386 exp2 += norm_Xsig(&accumulator);
387 shr_Xsig(&accumulator, 1);
388 exp2++;
389 shr_Xsig(&fix_up, 65 + exp2);
390
391 add_Xsig_Xsig(&accumulator, &fix_up);
392
393 echange = round_Xsig(&accumulator);
394
395 result->exp = exp2 + EXP_BIAS + echange;
396 *(short *)&(result->sign) = 0;
397 significand(result) = XSIG_LL(accumulator);
398 }
399
400 #ifdef PARANOID
401 if ( (result->exp >= EXP_BIAS)
402 && (significand(result) > 0x8000000000000000LL) )
403 {
404 EXCEPTION(EX_INTERNAL|0x151);
405 }
406 #endif PARANOID
407
408 }