1 // SPDX-License-Identifier: GPL-2.0
2 /*---------------------------------------------------------------------------+
3 | fpu_trig.c |
4 | |
5 | Implementation of the FPU "transcendental" functions. |
6 | |
7 | Copyright (C) 1992,1993,1994,1997,1999 |
8 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
9 | Australia. E-mail billm@melbpc.org.au |
10 | |
11 | |
12 +---------------------------------------------------------------------------*/
13
14 #include "fpu_system.h"
15 #include "exception.h"
16 #include "fpu_emu.h"
17 #include "status_w.h"
18 #include "control_w.h"
19 #include "reg_constant.h"
20
21 static void rem_kernel(unsigned long long st0, unsigned long long *y,
22 unsigned long long st1, unsigned long long q, int n);
23
24 #define BETTER_THAN_486
25
26 #define FCOS 4
27
28 /* Used only by fptan, fsin, fcos, and fsincos. */
29 /* This routine produces very accurate results, similar to
30 using a value of pi with more than 128 bits precision. */
31 /* Limited measurements show no results worse than 64 bit precision
32 except for the results for arguments close to 2^63, where the
33 precision of the result sometimes degrades to about 63.9 bits */
trig_arg(FPU_REG * st0_ptr,int even)34 static int trig_arg(FPU_REG *st0_ptr, int even)
35 {
36 FPU_REG tmp;
37 u_char tmptag;
38 unsigned long long q;
39 int old_cw = control_word, saved_status = partial_status;
40 int tag, st0_tag = TAG_Valid;
41
42 if (exponent(st0_ptr) >= 63) {
43 partial_status |= SW_C2; /* Reduction incomplete. */
44 return -1;
45 }
46
47 control_word &= ~CW_RC;
48 control_word |= RC_CHOP;
49
50 setpositive(st0_ptr);
51 tag = FPU_u_div(st0_ptr, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f,
52 SIGN_POS);
53
54 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't overflow
55 to 2^64 */
56 q = significand(&tmp);
57 if (q) {
58 rem_kernel(significand(st0_ptr),
59 &significand(&tmp),
60 significand(&CONST_PI2),
61 q, exponent(st0_ptr) - exponent(&CONST_PI2));
62 setexponent16(&tmp, exponent(&CONST_PI2));
63 st0_tag = FPU_normalize(&tmp);
64 FPU_copy_to_reg0(&tmp, st0_tag);
65 }
66
67 if ((even && !(q & 1)) || (!even && (q & 1))) {
68 st0_tag =
69 FPU_sub(REV | LOADED | TAG_Valid, (int)&CONST_PI2,
70 FULL_PRECISION);
71
72 #ifdef BETTER_THAN_486
73 /* So far, the results are exact but based upon a 64 bit
74 precision approximation to pi/2. The technique used
75 now is equivalent to using an approximation to pi/2 which
76 is accurate to about 128 bits. */
77 if ((exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64)
78 || (q > 1)) {
79 /* This code gives the effect of having pi/2 to better than
80 128 bits precision. */
81
82 significand(&tmp) = q + 1;
83 setexponent16(&tmp, 63);
84 FPU_normalize(&tmp);
85 tmptag =
86 FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
87 FULL_PRECISION, SIGN_POS,
88 exponent(&CONST_PI2extra) +
89 exponent(&tmp));
90 setsign(&tmp, getsign(&CONST_PI2extra));
91 st0_tag = FPU_add(&tmp, tmptag, 0, FULL_PRECISION);
92 if (signnegative(st0_ptr)) {
93 /* CONST_PI2extra is negative, so the result of the addition
94 can be negative. This means that the argument is actually
95 in a different quadrant. The correction is always < pi/2,
96 so it can't overflow into yet another quadrant. */
97 setpositive(st0_ptr);
98 q++;
99 }
100 }
101 #endif /* BETTER_THAN_486 */
102 }
103 #ifdef BETTER_THAN_486
104 else {
105 /* So far, the results are exact but based upon a 64 bit
106 precision approximation to pi/2. The technique used
107 now is equivalent to using an approximation to pi/2 which
108 is accurate to about 128 bits. */
109 if (((q > 0)
110 && (exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64))
111 || (q > 1)) {
112 /* This code gives the effect of having p/2 to better than
113 128 bits precision. */
114
115 significand(&tmp) = q;
116 setexponent16(&tmp, 63);
117 FPU_normalize(&tmp); /* This must return TAG_Valid */
118 tmptag =
119 FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
120 FULL_PRECISION, SIGN_POS,
121 exponent(&CONST_PI2extra) +
122 exponent(&tmp));
123 setsign(&tmp, getsign(&CONST_PI2extra));
124 st0_tag = FPU_sub(LOADED | (tmptag & 0x0f), (int)&tmp,
125 FULL_PRECISION);
126 if ((exponent(st0_ptr) == exponent(&CONST_PI2)) &&
127 ((st0_ptr->sigh > CONST_PI2.sigh)
128 || ((st0_ptr->sigh == CONST_PI2.sigh)
129 && (st0_ptr->sigl > CONST_PI2.sigl)))) {
130 /* CONST_PI2extra is negative, so the result of the
131 subtraction can be larger than pi/2. This means
132 that the argument is actually in a different quadrant.
133 The correction is always < pi/2, so it can't overflow
134 into yet another quadrant. */
135 st0_tag =
136 FPU_sub(REV | LOADED | TAG_Valid,
137 (int)&CONST_PI2, FULL_PRECISION);
138 q++;
139 }
140 }
141 }
142 #endif /* BETTER_THAN_486 */
143
144 FPU_settag0(st0_tag);
145 control_word = old_cw;
146 partial_status = saved_status & ~SW_C2; /* Reduction complete. */
147
148 return (q & 3) | even;
149 }
150
151 /* Convert a long to register */
convert_l2reg(long const * arg,int deststnr)152 static void convert_l2reg(long const *arg, int deststnr)
153 {
154 int tag;
155 long num = *arg;
156 u_char sign;
157 FPU_REG *dest = &st(deststnr);
158
159 if (num == 0) {
160 FPU_copy_to_regi(&CONST_Z, TAG_Zero, deststnr);
161 return;
162 }
163
164 if (num > 0) {
165 sign = SIGN_POS;
166 } else {
167 num = -num;
168 sign = SIGN_NEG;
169 }
170
171 dest->sigh = num;
172 dest->sigl = 0;
173 setexponent16(dest, 31);
174 tag = FPU_normalize(dest);
175 FPU_settagi(deststnr, tag);
176 setsign(dest, sign);
177 return;
178 }
179
single_arg_error(FPU_REG * st0_ptr,u_char st0_tag)180 static void single_arg_error(FPU_REG *st0_ptr, u_char st0_tag)
181 {
182 if (st0_tag == TAG_Empty)
183 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
184 else if (st0_tag == TW_NaN)
185 real_1op_NaN(st0_ptr); /* return with a NaN in st(0) */
186 #ifdef PARANOID
187 else
188 EXCEPTION(EX_INTERNAL | 0x0112);
189 #endif /* PARANOID */
190 }
191
single_arg_2_error(FPU_REG * st0_ptr,u_char st0_tag)192 static void single_arg_2_error(FPU_REG *st0_ptr, u_char st0_tag)
193 {
194 int isNaN;
195
196 switch (st0_tag) {
197 case TW_NaN:
198 isNaN = (exponent(st0_ptr) == EXP_OVER)
199 && (st0_ptr->sigh & 0x80000000);
200 if (isNaN && !(st0_ptr->sigh & 0x40000000)) { /* Signaling ? */
201 EXCEPTION(EX_Invalid);
202 if (control_word & CW_Invalid) {
203 /* The masked response */
204 /* Convert to a QNaN */
205 st0_ptr->sigh |= 0x40000000;
206 push();
207 FPU_copy_to_reg0(st0_ptr, TAG_Special);
208 }
209 } else if (isNaN) {
210 /* A QNaN */
211 push();
212 FPU_copy_to_reg0(st0_ptr, TAG_Special);
213 } else {
214 /* pseudoNaN or other unsupported */
215 EXCEPTION(EX_Invalid);
216 if (control_word & CW_Invalid) {
217 /* The masked response */
218 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
219 push();
220 FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
221 }
222 }
223 break; /* return with a NaN in st(0) */
224 #ifdef PARANOID
225 default:
226 EXCEPTION(EX_INTERNAL | 0x0112);
227 #endif /* PARANOID */
228 }
229 }
230
231 /*---------------------------------------------------------------------------*/
232
f2xm1(FPU_REG * st0_ptr,u_char tag)233 static void f2xm1(FPU_REG *st0_ptr, u_char tag)
234 {
235 FPU_REG a;
236
237 clear_C1();
238
239 if (tag == TAG_Valid) {
240 /* For an 80486 FPU, the result is undefined if the arg is >= 1.0 */
241 if (exponent(st0_ptr) < 0) {
242 denormal_arg:
243
244 FPU_to_exp16(st0_ptr, &a);
245
246 /* poly_2xm1(x) requires 0 < st(0) < 1. */
247 poly_2xm1(getsign(st0_ptr), &a, st0_ptr);
248 }
249 set_precision_flag_up(); /* 80486 appears to always do this */
250 return;
251 }
252
253 if (tag == TAG_Zero)
254 return;
255
256 if (tag == TAG_Special)
257 tag = FPU_Special(st0_ptr);
258
259 switch (tag) {
260 case TW_Denormal:
261 if (denormal_operand() < 0)
262 return;
263 goto denormal_arg;
264 case TW_Infinity:
265 if (signnegative(st0_ptr)) {
266 /* -infinity gives -1 (p16-10) */
267 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
268 setnegative(st0_ptr);
269 }
270 return;
271 default:
272 single_arg_error(st0_ptr, tag);
273 }
274 }
275
fptan(FPU_REG * st0_ptr,u_char st0_tag)276 static void fptan(FPU_REG *st0_ptr, u_char st0_tag)
277 {
278 FPU_REG *st_new_ptr;
279 int q;
280 u_char arg_sign = getsign(st0_ptr);
281
282 /* Stack underflow has higher priority */
283 if (st0_tag == TAG_Empty) {
284 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
285 if (control_word & CW_Invalid) {
286 st_new_ptr = &st(-1);
287 push();
288 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
289 }
290 return;
291 }
292
293 if (STACK_OVERFLOW) {
294 FPU_stack_overflow();
295 return;
296 }
297
298 if (st0_tag == TAG_Valid) {
299 if (exponent(st0_ptr) > -40) {
300 if ((q = trig_arg(st0_ptr, 0)) == -1) {
301 /* Operand is out of range */
302 return;
303 }
304
305 poly_tan(st0_ptr);
306 setsign(st0_ptr, (q & 1) ^ (arg_sign != 0));
307 set_precision_flag_up(); /* We do not really know if up or down */
308 } else {
309 /* For a small arg, the result == the argument */
310 /* Underflow may happen */
311
312 denormal_arg:
313
314 FPU_to_exp16(st0_ptr, st0_ptr);
315
316 st0_tag =
317 FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
318 FPU_settag0(st0_tag);
319 }
320 push();
321 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
322 return;
323 }
324
325 if (st0_tag == TAG_Zero) {
326 push();
327 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
328 setcc(0);
329 return;
330 }
331
332 if (st0_tag == TAG_Special)
333 st0_tag = FPU_Special(st0_ptr);
334
335 if (st0_tag == TW_Denormal) {
336 if (denormal_operand() < 0)
337 return;
338
339 goto denormal_arg;
340 }
341
342 if (st0_tag == TW_Infinity) {
343 /* The 80486 treats infinity as an invalid operand */
344 if (arith_invalid(0) >= 0) {
345 st_new_ptr = &st(-1);
346 push();
347 arith_invalid(0);
348 }
349 return;
350 }
351
352 single_arg_2_error(st0_ptr, st0_tag);
353 }
354
fxtract(FPU_REG * st0_ptr,u_char st0_tag)355 static void fxtract(FPU_REG *st0_ptr, u_char st0_tag)
356 {
357 FPU_REG *st_new_ptr;
358 u_char sign;
359 register FPU_REG *st1_ptr = st0_ptr; /* anticipate */
360
361 if (STACK_OVERFLOW) {
362 FPU_stack_overflow();
363 return;
364 }
365
366 clear_C1();
367
368 if (st0_tag == TAG_Valid) {
369 long e;
370
371 push();
372 sign = getsign(st1_ptr);
373 reg_copy(st1_ptr, st_new_ptr);
374 setexponent16(st_new_ptr, exponent(st_new_ptr));
375
376 denormal_arg:
377
378 e = exponent16(st_new_ptr);
379 convert_l2reg(&e, 1);
380 setexponentpos(st_new_ptr, 0);
381 setsign(st_new_ptr, sign);
382 FPU_settag0(TAG_Valid); /* Needed if arg was a denormal */
383 return;
384 } else if (st0_tag == TAG_Zero) {
385 sign = getsign(st0_ptr);
386
387 if (FPU_divide_by_zero(0, SIGN_NEG) < 0)
388 return;
389
390 push();
391 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
392 setsign(st_new_ptr, sign);
393 return;
394 }
395
396 if (st0_tag == TAG_Special)
397 st0_tag = FPU_Special(st0_ptr);
398
399 if (st0_tag == TW_Denormal) {
400 if (denormal_operand() < 0)
401 return;
402
403 push();
404 sign = getsign(st1_ptr);
405 FPU_to_exp16(st1_ptr, st_new_ptr);
406 goto denormal_arg;
407 } else if (st0_tag == TW_Infinity) {
408 sign = getsign(st0_ptr);
409 setpositive(st0_ptr);
410 push();
411 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
412 setsign(st_new_ptr, sign);
413 return;
414 } else if (st0_tag == TW_NaN) {
415 if (real_1op_NaN(st0_ptr) < 0)
416 return;
417
418 push();
419 FPU_copy_to_reg0(st0_ptr, TAG_Special);
420 return;
421 } else if (st0_tag == TAG_Empty) {
422 /* Is this the correct behaviour? */
423 if (control_word & EX_Invalid) {
424 FPU_stack_underflow();
425 push();
426 FPU_stack_underflow();
427 } else
428 EXCEPTION(EX_StackUnder);
429 }
430 #ifdef PARANOID
431 else
432 EXCEPTION(EX_INTERNAL | 0x119);
433 #endif /* PARANOID */
434 }
435
fdecstp(void)436 static void fdecstp(void)
437 {
438 clear_C1();
439 top--;
440 }
441
fincstp(void)442 static void fincstp(void)
443 {
444 clear_C1();
445 top++;
446 }
447
fsqrt_(FPU_REG * st0_ptr,u_char st0_tag)448 static void fsqrt_(FPU_REG *st0_ptr, u_char st0_tag)
449 {
450 int expon;
451
452 clear_C1();
453
454 if (st0_tag == TAG_Valid) {
455 u_char tag;
456
457 if (signnegative(st0_ptr)) {
458 arith_invalid(0); /* sqrt(negative) is invalid */
459 return;
460 }
461
462 /* make st(0) in [1.0 .. 4.0) */
463 expon = exponent(st0_ptr);
464
465 denormal_arg:
466
467 setexponent16(st0_ptr, (expon & 1));
468
469 /* Do the computation, the sign of the result will be positive. */
470 tag = wm_sqrt(st0_ptr, 0, 0, control_word, SIGN_POS);
471 addexponent(st0_ptr, expon >> 1);
472 FPU_settag0(tag);
473 return;
474 }
475
476 if (st0_tag == TAG_Zero)
477 return;
478
479 if (st0_tag == TAG_Special)
480 st0_tag = FPU_Special(st0_ptr);
481
482 if (st0_tag == TW_Infinity) {
483 if (signnegative(st0_ptr))
484 arith_invalid(0); /* sqrt(-Infinity) is invalid */
485 return;
486 } else if (st0_tag == TW_Denormal) {
487 if (signnegative(st0_ptr)) {
488 arith_invalid(0); /* sqrt(negative) is invalid */
489 return;
490 }
491
492 if (denormal_operand() < 0)
493 return;
494
495 FPU_to_exp16(st0_ptr, st0_ptr);
496
497 expon = exponent16(st0_ptr);
498
499 goto denormal_arg;
500 }
501
502 single_arg_error(st0_ptr, st0_tag);
503
504 }
505
frndint_(FPU_REG * st0_ptr,u_char st0_tag)506 static void frndint_(FPU_REG *st0_ptr, u_char st0_tag)
507 {
508 int flags, tag;
509
510 if (st0_tag == TAG_Valid) {
511 u_char sign;
512
513 denormal_arg:
514
515 sign = getsign(st0_ptr);
516
517 if (exponent(st0_ptr) > 63)
518 return;
519
520 if (st0_tag == TW_Denormal) {
521 if (denormal_operand() < 0)
522 return;
523 }
524
525 /* Fortunately, this can't overflow to 2^64 */
526 if ((flags = FPU_round_to_int(st0_ptr, st0_tag)))
527 set_precision_flag(flags);
528
529 setexponent16(st0_ptr, 63);
530 tag = FPU_normalize(st0_ptr);
531 setsign(st0_ptr, sign);
532 FPU_settag0(tag);
533 return;
534 }
535
536 if (st0_tag == TAG_Zero)
537 return;
538
539 if (st0_tag == TAG_Special)
540 st0_tag = FPU_Special(st0_ptr);
541
542 if (st0_tag == TW_Denormal)
543 goto denormal_arg;
544 else if (st0_tag == TW_Infinity)
545 return;
546 else
547 single_arg_error(st0_ptr, st0_tag);
548 }
549
f_sin(FPU_REG * st0_ptr,u_char tag)550 static int f_sin(FPU_REG *st0_ptr, u_char tag)
551 {
552 u_char arg_sign = getsign(st0_ptr);
553
554 if (tag == TAG_Valid) {
555 int q;
556
557 if (exponent(st0_ptr) > -40) {
558 if ((q = trig_arg(st0_ptr, 0)) == -1) {
559 /* Operand is out of range */
560 return 1;
561 }
562
563 poly_sine(st0_ptr);
564
565 if (q & 2)
566 changesign(st0_ptr);
567
568 setsign(st0_ptr, getsign(st0_ptr) ^ arg_sign);
569
570 /* We do not really know if up or down */
571 set_precision_flag_up();
572 return 0;
573 } else {
574 /* For a small arg, the result == the argument */
575 set_precision_flag_up(); /* Must be up. */
576 return 0;
577 }
578 }
579
580 if (tag == TAG_Zero) {
581 setcc(0);
582 return 0;
583 }
584
585 if (tag == TAG_Special)
586 tag = FPU_Special(st0_ptr);
587
588 if (tag == TW_Denormal) {
589 if (denormal_operand() < 0)
590 return 1;
591
592 /* For a small arg, the result == the argument */
593 /* Underflow may happen */
594 FPU_to_exp16(st0_ptr, st0_ptr);
595
596 tag = FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
597
598 FPU_settag0(tag);
599
600 return 0;
601 } else if (tag == TW_Infinity) {
602 /* The 80486 treats infinity as an invalid operand */
603 arith_invalid(0);
604 return 1;
605 } else {
606 single_arg_error(st0_ptr, tag);
607 return 1;
608 }
609 }
610
fsin(FPU_REG * st0_ptr,u_char tag)611 static void fsin(FPU_REG *st0_ptr, u_char tag)
612 {
613 f_sin(st0_ptr, tag);
614 }
615
f_cos(FPU_REG * st0_ptr,u_char tag)616 static int f_cos(FPU_REG *st0_ptr, u_char tag)
617 {
618 u_char st0_sign;
619
620 st0_sign = getsign(st0_ptr);
621
622 if (tag == TAG_Valid) {
623 int q;
624
625 if (exponent(st0_ptr) > -40) {
626 if ((exponent(st0_ptr) < 0)
627 || ((exponent(st0_ptr) == 0)
628 && (significand(st0_ptr) <=
629 0xc90fdaa22168c234LL))) {
630 poly_cos(st0_ptr);
631
632 /* We do not really know if up or down */
633 set_precision_flag_down();
634
635 return 0;
636 } else if ((q = trig_arg(st0_ptr, FCOS)) != -1) {
637 poly_sine(st0_ptr);
638
639 if ((q + 1) & 2)
640 changesign(st0_ptr);
641
642 /* We do not really know if up or down */
643 set_precision_flag_down();
644
645 return 0;
646 } else {
647 /* Operand is out of range */
648 return 1;
649 }
650 } else {
651 denormal_arg:
652
653 setcc(0);
654 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
655 #ifdef PECULIAR_486
656 set_precision_flag_down(); /* 80486 appears to do this. */
657 #else
658 set_precision_flag_up(); /* Must be up. */
659 #endif /* PECULIAR_486 */
660 return 0;
661 }
662 } else if (tag == TAG_Zero) {
663 FPU_copy_to_reg0(&CONST_1, TAG_Valid);
664 setcc(0);
665 return 0;
666 }
667
668 if (tag == TAG_Special)
669 tag = FPU_Special(st0_ptr);
670
671 if (tag == TW_Denormal) {
672 if (denormal_operand() < 0)
673 return 1;
674
675 goto denormal_arg;
676 } else if (tag == TW_Infinity) {
677 /* The 80486 treats infinity as an invalid operand */
678 arith_invalid(0);
679 return 1;
680 } else {
681 single_arg_error(st0_ptr, tag); /* requires st0_ptr == &st(0) */
682 return 1;
683 }
684 }
685
fcos(FPU_REG * st0_ptr,u_char st0_tag)686 static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
687 {
688 f_cos(st0_ptr, st0_tag);
689 }
690
fsincos(FPU_REG * st0_ptr,u_char st0_tag)691 static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
692 {
693 FPU_REG *st_new_ptr;
694 FPU_REG arg;
695 u_char tag;
696
697 /* Stack underflow has higher priority */
698 if (st0_tag == TAG_Empty) {
699 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
700 if (control_word & CW_Invalid) {
701 st_new_ptr = &st(-1);
702 push();
703 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
704 }
705 return;
706 }
707
708 if (STACK_OVERFLOW) {
709 FPU_stack_overflow();
710 return;
711 }
712
713 if (st0_tag == TAG_Special)
714 tag = FPU_Special(st0_ptr);
715 else
716 tag = st0_tag;
717
718 if (tag == TW_NaN) {
719 single_arg_2_error(st0_ptr, TW_NaN);
720 return;
721 } else if (tag == TW_Infinity) {
722 /* The 80486 treats infinity as an invalid operand */
723 if (arith_invalid(0) >= 0) {
724 /* Masked response */
725 push();
726 arith_invalid(0);
727 }
728 return;
729 }
730
731 reg_copy(st0_ptr, &arg);
732 if (!f_sin(st0_ptr, st0_tag)) {
733 push();
734 FPU_copy_to_reg0(&arg, st0_tag);
735 f_cos(&st(0), st0_tag);
736 } else {
737 /* An error, so restore st(0) */
738 FPU_copy_to_reg0(&arg, st0_tag);
739 }
740 }
741
742 /*---------------------------------------------------------------------------*/
743 /* The following all require two arguments: st(0) and st(1) */
744
745 /* A lean, mean kernel for the fprem instructions. This relies upon
746 the division and rounding to an integer in do_fprem giving an
747 exact result. Because of this, rem_kernel() needs to deal only with
748 the least significant 64 bits, the more significant bits of the
749 result must be zero.
750 */
rem_kernel(unsigned long long st0,unsigned long long * y,unsigned long long st1,unsigned long long q,int n)751 static void rem_kernel(unsigned long long st0, unsigned long long *y,
752 unsigned long long st1, unsigned long long q, int n)
753 {
754 int dummy;
755 unsigned long long x;
756
757 x = st0 << n;
758
759 /* Do the required multiplication and subtraction in the one operation */
760
761 /* lsw x -= lsw st1 * lsw q */
762 asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
763 (((unsigned *)&x)[0]), "=m"(((unsigned *)&x)[1]),
764 "=a"(dummy)
765 :"2"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[0])
766 :"%dx");
767 /* msw x -= msw st1 * lsw q */
768 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
769 "=a"(dummy)
770 :"1"(((unsigned *)&st1)[1]), "m"(((unsigned *)&q)[0])
771 :"%dx");
772 /* msw x -= lsw st1 * msw q */
773 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x)[1]),
774 "=a"(dummy)
775 :"1"(((unsigned *)&st1)[0]), "m"(((unsigned *)&q)[1])
776 :"%dx");
777
778 *y = x;
779 }
780
781 /* Remainder of st(0) / st(1) */
782 /* This routine produces exact results, i.e. there is never any
783 rounding or truncation, etc of the result. */
do_fprem(FPU_REG * st0_ptr,u_char st0_tag,int round)784 static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
785 {
786 FPU_REG *st1_ptr = &st(1);
787 u_char st1_tag = FPU_gettagi(1);
788
789 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
790 FPU_REG tmp, st0, st1;
791 u_char st0_sign, st1_sign;
792 u_char tmptag;
793 int tag;
794 int old_cw;
795 int expdif;
796 long long q;
797 unsigned short saved_status;
798 int cc;
799
800 fprem_valid:
801 /* Convert registers for internal use. */
802 st0_sign = FPU_to_exp16(st0_ptr, &st0);
803 st1_sign = FPU_to_exp16(st1_ptr, &st1);
804 expdif = exponent16(&st0) - exponent16(&st1);
805
806 old_cw = control_word;
807 cc = 0;
808
809 /* We want the status following the denorm tests, but don't want
810 the status changed by the arithmetic operations. */
811 saved_status = partial_status;
812 control_word &= ~CW_RC;
813 control_word |= RC_CHOP;
814
815 if (expdif < 64) {
816 /* This should be the most common case */
817
818 if (expdif > -2) {
819 u_char sign = st0_sign ^ st1_sign;
820 tag = FPU_u_div(&st0, &st1, &tmp,
821 PR_64_BITS | RC_CHOP | 0x3f,
822 sign);
823 setsign(&tmp, sign);
824
825 if (exponent(&tmp) >= 0) {
826 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't
827 overflow to 2^64 */
828 q = significand(&tmp);
829
830 rem_kernel(significand(&st0),
831 &significand(&tmp),
832 significand(&st1),
833 q, expdif);
834
835 setexponent16(&tmp, exponent16(&st1));
836 } else {
837 reg_copy(&st0, &tmp);
838 q = 0;
839 }
840
841 if ((round == RC_RND)
842 && (tmp.sigh & 0xc0000000)) {
843 /* We may need to subtract st(1) once more,
844 to get a result <= 1/2 of st(1). */
845 unsigned long long x;
846 expdif =
847 exponent16(&st1) - exponent16(&tmp);
848 if (expdif <= 1) {
849 if (expdif == 0)
850 x = significand(&st1) -
851 significand(&tmp);
852 else /* expdif is 1 */
853 x = (significand(&st1)
854 << 1) -
855 significand(&tmp);
856 if ((x < significand(&tmp)) ||
857 /* or equi-distant (from 0 & st(1)) and q is odd */
858 ((x == significand(&tmp))
859 && (q & 1))) {
860 st0_sign = !st0_sign;
861 significand(&tmp) = x;
862 q++;
863 }
864 }
865 }
866
867 if (q & 4)
868 cc |= SW_C0;
869 if (q & 2)
870 cc |= SW_C3;
871 if (q & 1)
872 cc |= SW_C1;
873 } else {
874 control_word = old_cw;
875 setcc(0);
876 return;
877 }
878 } else {
879 /* There is a large exponent difference ( >= 64 ) */
880 /* To make much sense, the code in this section should
881 be done at high precision. */
882 int exp_1, N;
883 u_char sign;
884
885 /* prevent overflow here */
886 /* N is 'a number between 32 and 63' (p26-113) */
887 reg_copy(&st0, &tmp);
888 tmptag = st0_tag;
889 N = (expdif & 0x0000001f) + 32; /* This choice gives results
890 identical to an AMD 486 */
891 setexponent16(&tmp, N);
892 exp_1 = exponent16(&st1);
893 setexponent16(&st1, 0);
894 expdif -= N;
895
896 sign = getsign(&tmp) ^ st1_sign;
897 tag =
898 FPU_u_div(&tmp, &st1, &tmp,
899 PR_64_BITS | RC_CHOP | 0x3f, sign);
900 setsign(&tmp, sign);
901
902 FPU_round_to_int(&tmp, tag); /* Fortunately, this can't
903 overflow to 2^64 */
904
905 rem_kernel(significand(&st0),
906 &significand(&tmp),
907 significand(&st1),
908 significand(&tmp), exponent(&tmp)
909 );
910 setexponent16(&tmp, exp_1 + expdif);
911
912 /* It is possible for the operation to be complete here.
913 What does the IEEE standard say? The Intel 80486 manual
914 implies that the operation will never be completed at this
915 point, and the behaviour of a real 80486 confirms this.
916 */
917 if (!(tmp.sigh | tmp.sigl)) {
918 /* The result is zero */
919 control_word = old_cw;
920 partial_status = saved_status;
921 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
922 setsign(&st0, st0_sign);
923 #ifdef PECULIAR_486
924 setcc(SW_C2);
925 #else
926 setcc(0);
927 #endif /* PECULIAR_486 */
928 return;
929 }
930 cc = SW_C2;
931 }
932
933 control_word = old_cw;
934 partial_status = saved_status;
935 tag = FPU_normalize_nuo(&tmp);
936 reg_copy(&tmp, st0_ptr);
937
938 /* The only condition to be looked for is underflow,
939 and it can occur here only if underflow is unmasked. */
940 if ((exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
941 && !(control_word & CW_Underflow)) {
942 setcc(cc);
943 tag = arith_underflow(st0_ptr);
944 setsign(st0_ptr, st0_sign);
945 FPU_settag0(tag);
946 return;
947 } else if ((exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero)) {
948 stdexp(st0_ptr);
949 setsign(st0_ptr, st0_sign);
950 } else {
951 tag =
952 FPU_round(st0_ptr, 0, 0, FULL_PRECISION, st0_sign);
953 }
954 FPU_settag0(tag);
955 setcc(cc);
956
957 return;
958 }
959
960 if (st0_tag == TAG_Special)
961 st0_tag = FPU_Special(st0_ptr);
962 if (st1_tag == TAG_Special)
963 st1_tag = FPU_Special(st1_ptr);
964
965 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
966 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
967 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
968 if (denormal_operand() < 0)
969 return;
970 goto fprem_valid;
971 } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
972 FPU_stack_underflow();
973 return;
974 } else if (st0_tag == TAG_Zero) {
975 if (st1_tag == TAG_Valid) {
976 setcc(0);
977 return;
978 } else if (st1_tag == TW_Denormal) {
979 if (denormal_operand() < 0)
980 return;
981 setcc(0);
982 return;
983 } else if (st1_tag == TAG_Zero) {
984 arith_invalid(0);
985 return;
986 } /* fprem(?,0) always invalid */
987 else if (st1_tag == TW_Infinity) {
988 setcc(0);
989 return;
990 }
991 } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
992 if (st1_tag == TAG_Zero) {
993 arith_invalid(0); /* fprem(Valid,Zero) is invalid */
994 return;
995 } else if (st1_tag != TW_NaN) {
996 if (((st0_tag == TW_Denormal)
997 || (st1_tag == TW_Denormal))
998 && (denormal_operand() < 0))
999 return;
1000
1001 if (st1_tag == TW_Infinity) {
1002 /* fprem(Valid,Infinity) is o.k. */
1003 setcc(0);
1004 return;
1005 }
1006 }
1007 } else if (st0_tag == TW_Infinity) {
1008 if (st1_tag != TW_NaN) {
1009 arith_invalid(0); /* fprem(Infinity,?) is invalid */
1010 return;
1011 }
1012 }
1013
1014 /* One of the registers must contain a NaN if we got here. */
1015
1016 #ifdef PARANOID
1017 if ((st0_tag != TW_NaN) && (st1_tag != TW_NaN))
1018 EXCEPTION(EX_INTERNAL | 0x118);
1019 #endif /* PARANOID */
1020
1021 real_2op_NaN(st1_ptr, st1_tag, 0, st1_ptr);
1022
1023 }
1024
1025 /* ST(1) <- ST(1) * log ST; pop ST */
fyl2x(FPU_REG * st0_ptr,u_char st0_tag)1026 static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
1027 {
1028 FPU_REG *st1_ptr = &st(1), exponent;
1029 u_char st1_tag = FPU_gettagi(1);
1030 u_char sign;
1031 int e, tag;
1032
1033 clear_C1();
1034
1035 if ((st0_tag == TAG_Valid) && (st1_tag == TAG_Valid)) {
1036 both_valid:
1037 /* Both regs are Valid or Denormal */
1038 if (signpositive(st0_ptr)) {
1039 if (st0_tag == TW_Denormal)
1040 FPU_to_exp16(st0_ptr, st0_ptr);
1041 else
1042 /* Convert st(0) for internal use. */
1043 setexponent16(st0_ptr, exponent(st0_ptr));
1044
1045 if ((st0_ptr->sigh == 0x80000000)
1046 && (st0_ptr->sigl == 0)) {
1047 /* Special case. The result can be precise. */
1048 u_char esign;
1049 e = exponent16(st0_ptr);
1050 if (e >= 0) {
1051 exponent.sigh = e;
1052 esign = SIGN_POS;
1053 } else {
1054 exponent.sigh = -e;
1055 esign = SIGN_NEG;
1056 }
1057 exponent.sigl = 0;
1058 setexponent16(&exponent, 31);
1059 tag = FPU_normalize_nuo(&exponent);
1060 stdexp(&exponent);
1061 setsign(&exponent, esign);
1062 tag =
1063 FPU_mul(&exponent, tag, 1, FULL_PRECISION);
1064 if (tag >= 0)
1065 FPU_settagi(1, tag);
1066 } else {
1067 /* The usual case */
1068 sign = getsign(st1_ptr);
1069 if (st1_tag == TW_Denormal)
1070 FPU_to_exp16(st1_ptr, st1_ptr);
1071 else
1072 /* Convert st(1) for internal use. */
1073 setexponent16(st1_ptr,
1074 exponent(st1_ptr));
1075 poly_l2(st0_ptr, st1_ptr, sign);
1076 }
1077 } else {
1078 /* negative */
1079 if (arith_invalid(1) < 0)
1080 return;
1081 }
1082
1083 FPU_pop();
1084
1085 return;
1086 }
1087
1088 if (st0_tag == TAG_Special)
1089 st0_tag = FPU_Special(st0_ptr);
1090 if (st1_tag == TAG_Special)
1091 st1_tag = FPU_Special(st1_ptr);
1092
1093 if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1094 FPU_stack_underflow_pop(1);
1095 return;
1096 } else if ((st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal)) {
1097 if (st0_tag == TAG_Zero) {
1098 if (st1_tag == TAG_Zero) {
1099 /* Both args zero is invalid */
1100 if (arith_invalid(1) < 0)
1101 return;
1102 } else {
1103 u_char sign;
1104 sign = getsign(st1_ptr) ^ SIGN_NEG;
1105 if (FPU_divide_by_zero(1, sign) < 0)
1106 return;
1107
1108 setsign(st1_ptr, sign);
1109 }
1110 } else if (st1_tag == TAG_Zero) {
1111 /* st(1) contains zero, st(0) valid <> 0 */
1112 /* Zero is the valid answer */
1113 sign = getsign(st1_ptr);
1114
1115 if (signnegative(st0_ptr)) {
1116 /* log(negative) */
1117 if (arith_invalid(1) < 0)
1118 return;
1119 } else if ((st0_tag == TW_Denormal)
1120 && (denormal_operand() < 0))
1121 return;
1122 else {
1123 if (exponent(st0_ptr) < 0)
1124 sign ^= SIGN_NEG;
1125
1126 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1127 setsign(st1_ptr, sign);
1128 }
1129 } else {
1130 /* One or both operands are denormals. */
1131 if (denormal_operand() < 0)
1132 return;
1133 goto both_valid;
1134 }
1135 } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1136 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1137 return;
1138 }
1139 /* One or both arg must be an infinity */
1140 else if (st0_tag == TW_Infinity) {
1141 if ((signnegative(st0_ptr)) || (st1_tag == TAG_Zero)) {
1142 /* log(-infinity) or 0*log(infinity) */
1143 if (arith_invalid(1) < 0)
1144 return;
1145 } else {
1146 u_char sign = getsign(st1_ptr);
1147
1148 if ((st1_tag == TW_Denormal)
1149 && (denormal_operand() < 0))
1150 return;
1151
1152 FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1153 setsign(st1_ptr, sign);
1154 }
1155 }
1156 /* st(1) must be infinity here */
1157 else if (((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal))
1158 && (signpositive(st0_ptr))) {
1159 if (exponent(st0_ptr) >= 0) {
1160 if ((exponent(st0_ptr) == 0) &&
1161 (st0_ptr->sigh == 0x80000000) &&
1162 (st0_ptr->sigl == 0)) {
1163 /* st(0) holds 1.0 */
1164 /* infinity*log(1) */
1165 if (arith_invalid(1) < 0)
1166 return;
1167 }
1168 /* else st(0) is positive and > 1.0 */
1169 } else {
1170 /* st(0) is positive and < 1.0 */
1171
1172 if ((st0_tag == TW_Denormal)
1173 && (denormal_operand() < 0))
1174 return;
1175
1176 changesign(st1_ptr);
1177 }
1178 } else {
1179 /* st(0) must be zero or negative */
1180 if (st0_tag == TAG_Zero) {
1181 /* This should be invalid, but a real 80486 is happy with it. */
1182
1183 #ifndef PECULIAR_486
1184 sign = getsign(st1_ptr);
1185 if (FPU_divide_by_zero(1, sign) < 0)
1186 return;
1187 #endif /* PECULIAR_486 */
1188
1189 changesign(st1_ptr);
1190 } else if (arith_invalid(1) < 0) /* log(negative) */
1191 return;
1192 }
1193
1194 FPU_pop();
1195 }
1196
fpatan(FPU_REG * st0_ptr,u_char st0_tag)1197 static void fpatan(FPU_REG *st0_ptr, u_char st0_tag)
1198 {
1199 FPU_REG *st1_ptr = &st(1);
1200 u_char st1_tag = FPU_gettagi(1);
1201 int tag;
1202
1203 clear_C1();
1204 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1205 valid_atan:
1206
1207 poly_atan(st0_ptr, st0_tag, st1_ptr, st1_tag);
1208
1209 FPU_pop();
1210
1211 return;
1212 }
1213
1214 if (st0_tag == TAG_Special)
1215 st0_tag = FPU_Special(st0_ptr);
1216 if (st1_tag == TAG_Special)
1217 st1_tag = FPU_Special(st1_ptr);
1218
1219 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1220 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1221 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1222 if (denormal_operand() < 0)
1223 return;
1224
1225 goto valid_atan;
1226 } else if ((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty)) {
1227 FPU_stack_underflow_pop(1);
1228 return;
1229 } else if ((st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
1230 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) >= 0)
1231 FPU_pop();
1232 return;
1233 } else if ((st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
1234 u_char sign = getsign(st1_ptr);
1235 if (st0_tag == TW_Infinity) {
1236 if (st1_tag == TW_Infinity) {
1237 if (signpositive(st0_ptr)) {
1238 FPU_copy_to_reg1(&CONST_PI4, TAG_Valid);
1239 } else {
1240 setpositive(st1_ptr);
1241 tag =
1242 FPU_u_add(&CONST_PI4, &CONST_PI2,
1243 st1_ptr, FULL_PRECISION,
1244 SIGN_POS,
1245 exponent(&CONST_PI4),
1246 exponent(&CONST_PI2));
1247 if (tag >= 0)
1248 FPU_settagi(1, tag);
1249 }
1250 } else {
1251 if ((st1_tag == TW_Denormal)
1252 && (denormal_operand() < 0))
1253 return;
1254
1255 if (signpositive(st0_ptr)) {
1256 FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
1257 setsign(st1_ptr, sign); /* An 80486 preserves the sign */
1258 FPU_pop();
1259 return;
1260 } else {
1261 FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1262 }
1263 }
1264 } else {
1265 /* st(1) is infinity, st(0) not infinity */
1266 if ((st0_tag == TW_Denormal)
1267 && (denormal_operand() < 0))
1268 return;
1269
1270 FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1271 }
1272 setsign(st1_ptr, sign);
1273 } else if (st1_tag == TAG_Zero) {
1274 /* st(0) must be valid or zero */
1275 u_char sign = getsign(st1_ptr);
1276
1277 if ((st0_tag == TW_Denormal) && (denormal_operand() < 0))
1278 return;
1279
1280 if (signpositive(st0_ptr)) {
1281 /* An 80486 preserves the sign */
1282 FPU_pop();
1283 return;
1284 }
1285
1286 FPU_copy_to_reg1(&CONST_PI, TAG_Valid);
1287 setsign(st1_ptr, sign);
1288 } else if (st0_tag == TAG_Zero) {
1289 /* st(1) must be TAG_Valid here */
1290 u_char sign = getsign(st1_ptr);
1291
1292 if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1293 return;
1294
1295 FPU_copy_to_reg1(&CONST_PI2, TAG_Valid);
1296 setsign(st1_ptr, sign);
1297 }
1298 #ifdef PARANOID
1299 else
1300 EXCEPTION(EX_INTERNAL | 0x125);
1301 #endif /* PARANOID */
1302
1303 FPU_pop();
1304 set_precision_flag_up(); /* We do not really know if up or down */
1305 }
1306
fprem(FPU_REG * st0_ptr,u_char st0_tag)1307 static void fprem(FPU_REG *st0_ptr, u_char st0_tag)
1308 {
1309 do_fprem(st0_ptr, st0_tag, RC_CHOP);
1310 }
1311
fprem1(FPU_REG * st0_ptr,u_char st0_tag)1312 static void fprem1(FPU_REG *st0_ptr, u_char st0_tag)
1313 {
1314 do_fprem(st0_ptr, st0_tag, RC_RND);
1315 }
1316
fyl2xp1(FPU_REG * st0_ptr,u_char st0_tag)1317 static void fyl2xp1(FPU_REG *st0_ptr, u_char st0_tag)
1318 {
1319 u_char sign, sign1;
1320 FPU_REG *st1_ptr = &st(1), a, b;
1321 u_char st1_tag = FPU_gettagi(1);
1322
1323 clear_C1();
1324 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1325 valid_yl2xp1:
1326
1327 sign = getsign(st0_ptr);
1328 sign1 = getsign(st1_ptr);
1329
1330 FPU_to_exp16(st0_ptr, &a);
1331 FPU_to_exp16(st1_ptr, &b);
1332
1333 if (poly_l2p1(sign, sign1, &a, &b, st1_ptr))
1334 return;
1335
1336 FPU_pop();
1337 return;
1338 }
1339
1340 if (st0_tag == TAG_Special)
1341 st0_tag = FPU_Special(st0_ptr);
1342 if (st1_tag == TAG_Special)
1343 st1_tag = FPU_Special(st1_ptr);
1344
1345 if (((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
1346 || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
1347 || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal))) {
1348 if (denormal_operand() < 0)
1349 return;
1350
1351 goto valid_yl2xp1;
1352 } else if ((st0_tag == TAG_Empty) | (st1_tag == TAG_Empty)) {
1353 FPU_stack_underflow_pop(1);
1354 return;
1355 } else if (st0_tag == TAG_Zero) {
1356 switch (st1_tag) {
1357 case TW_Denormal:
1358 if (denormal_operand() < 0)
1359 return;
1360 fallthrough;
1361 case TAG_Zero:
1362 case TAG_Valid:
1363 setsign(st0_ptr, getsign(st0_ptr) ^ getsign(st1_ptr));
1364 FPU_copy_to_reg1(st0_ptr, st0_tag);
1365 break;
1366
1367 case TW_Infinity:
1368 /* Infinity*log(1) */
1369 if (arith_invalid(1) < 0)
1370 return;
1371 break;
1372
1373 case TW_NaN:
1374 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1375 return;
1376 break;
1377
1378 default:
1379 #ifdef PARANOID
1380 EXCEPTION(EX_INTERNAL | 0x116);
1381 return;
1382 #endif /* PARANOID */
1383 break;
1384 }
1385 } else if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1386 switch (st1_tag) {
1387 case TAG_Zero:
1388 if (signnegative(st0_ptr)) {
1389 if (exponent(st0_ptr) >= 0) {
1390 /* st(0) holds <= -1.0 */
1391 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1392 changesign(st1_ptr);
1393 #else
1394 if (arith_invalid(1) < 0)
1395 return;
1396 #endif /* PECULIAR_486 */
1397 } else if ((st0_tag == TW_Denormal)
1398 && (denormal_operand() < 0))
1399 return;
1400 else
1401 changesign(st1_ptr);
1402 } else if ((st0_tag == TW_Denormal)
1403 && (denormal_operand() < 0))
1404 return;
1405 break;
1406
1407 case TW_Infinity:
1408 if (signnegative(st0_ptr)) {
1409 if ((exponent(st0_ptr) >= 0) &&
1410 !((st0_ptr->sigh == 0x80000000) &&
1411 (st0_ptr->sigl == 0))) {
1412 /* st(0) holds < -1.0 */
1413 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1414 changesign(st1_ptr);
1415 #else
1416 if (arith_invalid(1) < 0)
1417 return;
1418 #endif /* PECULIAR_486 */
1419 } else if ((st0_tag == TW_Denormal)
1420 && (denormal_operand() < 0))
1421 return;
1422 else
1423 changesign(st1_ptr);
1424 } else if ((st0_tag == TW_Denormal)
1425 && (denormal_operand() < 0))
1426 return;
1427 break;
1428
1429 case TW_NaN:
1430 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1431 return;
1432 }
1433
1434 } else if (st0_tag == TW_NaN) {
1435 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1436 return;
1437 } else if (st0_tag == TW_Infinity) {
1438 if (st1_tag == TW_NaN) {
1439 if (real_2op_NaN(st0_ptr, st0_tag, 1, st0_ptr) < 0)
1440 return;
1441 } else if (signnegative(st0_ptr)) {
1442 #ifndef PECULIAR_486
1443 /* This should have higher priority than denormals, but... */
1444 if (arith_invalid(1) < 0) /* log(-infinity) */
1445 return;
1446 #endif /* PECULIAR_486 */
1447 if ((st1_tag == TW_Denormal)
1448 && (denormal_operand() < 0))
1449 return;
1450 #ifdef PECULIAR_486
1451 /* Denormal operands actually get higher priority */
1452 if (arith_invalid(1) < 0) /* log(-infinity) */
1453 return;
1454 #endif /* PECULIAR_486 */
1455 } else if (st1_tag == TAG_Zero) {
1456 /* log(infinity) */
1457 if (arith_invalid(1) < 0)
1458 return;
1459 }
1460
1461 /* st(1) must be valid here. */
1462
1463 else if ((st1_tag == TW_Denormal) && (denormal_operand() < 0))
1464 return;
1465
1466 /* The Manual says that log(Infinity) is invalid, but a real
1467 80486 sensibly says that it is o.k. */
1468 else {
1469 u_char sign = getsign(st1_ptr);
1470 FPU_copy_to_reg1(&CONST_INF, TAG_Special);
1471 setsign(st1_ptr, sign);
1472 }
1473 }
1474 #ifdef PARANOID
1475 else {
1476 EXCEPTION(EX_INTERNAL | 0x117);
1477 return;
1478 }
1479 #endif /* PARANOID */
1480
1481 FPU_pop();
1482 return;
1483
1484 }
1485
fscale(FPU_REG * st0_ptr,u_char st0_tag)1486 static void fscale(FPU_REG *st0_ptr, u_char st0_tag)
1487 {
1488 FPU_REG *st1_ptr = &st(1);
1489 u_char st1_tag = FPU_gettagi(1);
1490 int old_cw = control_word;
1491 u_char sign = getsign(st0_ptr);
1492
1493 clear_C1();
1494 if (!((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid))) {
1495 long scale;
1496 FPU_REG tmp;
1497
1498 /* Convert register for internal use. */
1499 setexponent16(st0_ptr, exponent(st0_ptr));
1500
1501 valid_scale:
1502
1503 if (exponent(st1_ptr) > 30) {
1504 /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1505
1506 if (signpositive(st1_ptr)) {
1507 EXCEPTION(EX_Overflow);
1508 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1509 } else {
1510 EXCEPTION(EX_Underflow);
1511 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1512 }
1513 setsign(st0_ptr, sign);
1514 return;
1515 }
1516
1517 control_word &= ~CW_RC;
1518 control_word |= RC_CHOP;
1519 reg_copy(st1_ptr, &tmp);
1520 FPU_round_to_int(&tmp, st1_tag); /* This can never overflow here */
1521 control_word = old_cw;
1522 scale = signnegative(st1_ptr) ? -tmp.sigl : tmp.sigl;
1523 scale += exponent16(st0_ptr);
1524
1525 setexponent16(st0_ptr, scale);
1526
1527 /* Use FPU_round() to properly detect under/overflow etc */
1528 FPU_round(st0_ptr, 0, 0, control_word, sign);
1529
1530 return;
1531 }
1532
1533 if (st0_tag == TAG_Special)
1534 st0_tag = FPU_Special(st0_ptr);
1535 if (st1_tag == TAG_Special)
1536 st1_tag = FPU_Special(st1_ptr);
1537
1538 if ((st0_tag == TAG_Valid) || (st0_tag == TW_Denormal)) {
1539 switch (st1_tag) {
1540 case TAG_Valid:
1541 /* st(0) must be a denormal */
1542 if ((st0_tag == TW_Denormal)
1543 && (denormal_operand() < 0))
1544 return;
1545
1546 FPU_to_exp16(st0_ptr, st0_ptr); /* Will not be left on stack */
1547 goto valid_scale;
1548
1549 case TAG_Zero:
1550 if (st0_tag == TW_Denormal)
1551 denormal_operand();
1552 return;
1553
1554 case TW_Denormal:
1555 denormal_operand();
1556 return;
1557
1558 case TW_Infinity:
1559 if ((st0_tag == TW_Denormal)
1560 && (denormal_operand() < 0))
1561 return;
1562
1563 if (signpositive(st1_ptr))
1564 FPU_copy_to_reg0(&CONST_INF, TAG_Special);
1565 else
1566 FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
1567 setsign(st0_ptr, sign);
1568 return;
1569
1570 case TW_NaN:
1571 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1572 return;
1573 }
1574 } else if (st0_tag == TAG_Zero) {
1575 switch (st1_tag) {
1576 case TAG_Valid:
1577 case TAG_Zero:
1578 return;
1579
1580 case TW_Denormal:
1581 denormal_operand();
1582 return;
1583
1584 case TW_Infinity:
1585 if (signpositive(st1_ptr))
1586 arith_invalid(0); /* Zero scaled by +Infinity */
1587 return;
1588
1589 case TW_NaN:
1590 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1591 return;
1592 }
1593 } else if (st0_tag == TW_Infinity) {
1594 switch (st1_tag) {
1595 case TAG_Valid:
1596 case TAG_Zero:
1597 return;
1598
1599 case TW_Denormal:
1600 denormal_operand();
1601 return;
1602
1603 case TW_Infinity:
1604 if (signnegative(st1_ptr))
1605 arith_invalid(0); /* Infinity scaled by -Infinity */
1606 return;
1607
1608 case TW_NaN:
1609 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1610 return;
1611 }
1612 } else if (st0_tag == TW_NaN) {
1613 if (st1_tag != TAG_Empty) {
1614 real_2op_NaN(st1_ptr, st1_tag, 0, st0_ptr);
1615 return;
1616 }
1617 }
1618 #ifdef PARANOID
1619 if (!((st0_tag == TAG_Empty) || (st1_tag == TAG_Empty))) {
1620 EXCEPTION(EX_INTERNAL | 0x115);
1621 return;
1622 }
1623 #endif
1624
1625 /* At least one of st(0), st(1) must be empty */
1626 FPU_stack_underflow();
1627
1628 }
1629
1630 /*---------------------------------------------------------------------------*/
1631
1632 static FUNC_ST0 const trig_table_a[] = {
1633 f2xm1, fyl2x, fptan, fpatan,
1634 fxtract, fprem1, (FUNC_ST0) fdecstp, (FUNC_ST0) fincstp
1635 };
1636
FPU_triga(void)1637 void FPU_triga(void)
1638 {
1639 (trig_table_a[FPU_rm]) (&st(0), FPU_gettag0());
1640 }
1641
1642 static FUNC_ST0 const trig_table_b[] = {
1643 fprem, fyl2xp1, fsqrt_, fsincos, frndint_, fscale, fsin, fcos
1644 };
1645
FPU_trigb(void)1646 void FPU_trigb(void)
1647 {
1648 (trig_table_b[FPU_rm]) (&st(0), FPU_gettag0());
1649 }
1650