1 /*---------------------------------------------------------------------------+
4 | Implementation of the FPU "transcendental" functions. |
6 | Copyright (C) 1992,1993,1994,1997,1999 |
7 | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
8 | Australia. E-mail billm@melbpc.org.au |
11 +---------------------------------------------------------------------------*/
13 #include "fpu_system.h"
14 #include "exception.h"
17 #include "control_w.h"
18 #include "reg_constant.h"
20 static void rem_kernel(unsigned long long st0
, unsigned long long *y
,
21 unsigned long long st1
, unsigned long long q
, int n
);
23 #define BETTER_THAN_486
27 /* Used only by fptan, fsin, fcos, and fsincos. */
28 /* This routine produces very accurate results, similar to
29 using a value of pi with more than 128 bits precision. */
30 /* Limited measurements show no results worse than 64 bit precision
31 except for the results for arguments close to 2^63, where the
32 precision of the result sometimes degrades to about 63.9 bits */
33 static int trig_arg(FPU_REG
*st0_ptr
, int even
)
38 int old_cw
= control_word
, saved_status
= partial_status
;
39 int tag
, st0_tag
= TAG_Valid
;
41 if (exponent(st0_ptr
) >= 63) {
42 partial_status
|= SW_C2
; /* Reduction incomplete. */
46 control_word
&= ~CW_RC
;
47 control_word
|= RC_CHOP
;
50 tag
= FPU_u_div(st0_ptr
, &CONST_PI2
, &tmp
, PR_64_BITS
| RC_CHOP
| 0x3f,
53 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't overflow
55 q
= significand(&tmp
);
57 rem_kernel(significand(st0_ptr
),
59 significand(&CONST_PI2
),
60 q
, exponent(st0_ptr
) - exponent(&CONST_PI2
));
61 setexponent16(&tmp
, exponent(&CONST_PI2
));
62 st0_tag
= FPU_normalize(&tmp
);
63 FPU_copy_to_reg0(&tmp
, st0_tag
);
66 if ((even
&& !(q
& 1)) || (!even
&& (q
& 1))) {
68 FPU_sub(REV
| LOADED
| TAG_Valid
, (int)&CONST_PI2
,
71 #ifdef BETTER_THAN_486
72 /* So far, the results are exact but based upon a 64 bit
73 precision approximation to pi/2. The technique used
74 now is equivalent to using an approximation to pi/2 which
75 is accurate to about 128 bits. */
76 if ((exponent(st0_ptr
) <= exponent(&CONST_PI2extra
) + 64)
78 /* This code gives the effect of having pi/2 to better than
79 128 bits precision. */
81 significand(&tmp
) = q
+ 1;
82 setexponent16(&tmp
, 63);
85 FPU_u_mul(&CONST_PI2extra
, &tmp
, &tmp
,
86 FULL_PRECISION
, SIGN_POS
,
87 exponent(&CONST_PI2extra
) +
89 setsign(&tmp
, getsign(&CONST_PI2extra
));
90 st0_tag
= FPU_add(&tmp
, tmptag
, 0, FULL_PRECISION
);
91 if (signnegative(st0_ptr
)) {
92 /* CONST_PI2extra is negative, so the result of the addition
93 can be negative. This means that the argument is actually
94 in a different quadrant. The correction is always < pi/2,
95 so it can't overflow into yet another quadrant. */
100 #endif /* BETTER_THAN_486 */
102 #ifdef BETTER_THAN_486
104 /* So far, the results are exact but based upon a 64 bit
105 precision approximation to pi/2. The technique used
106 now is equivalent to using an approximation to pi/2 which
107 is accurate to about 128 bits. */
109 && (exponent(st0_ptr
) <= exponent(&CONST_PI2extra
) + 64))
111 /* This code gives the effect of having p/2 to better than
112 128 bits precision. */
114 significand(&tmp
) = q
;
115 setexponent16(&tmp
, 63);
116 FPU_normalize(&tmp
); /* This must return TAG_Valid */
118 FPU_u_mul(&CONST_PI2extra
, &tmp
, &tmp
,
119 FULL_PRECISION
, SIGN_POS
,
120 exponent(&CONST_PI2extra
) +
122 setsign(&tmp
, getsign(&CONST_PI2extra
));
123 st0_tag
= FPU_sub(LOADED
| (tmptag
& 0x0f), (int)&tmp
,
125 if ((exponent(st0_ptr
) == exponent(&CONST_PI2
)) &&
126 ((st0_ptr
->sigh
> CONST_PI2
.sigh
)
127 || ((st0_ptr
->sigh
== CONST_PI2
.sigh
)
128 && (st0_ptr
->sigl
> CONST_PI2
.sigl
)))) {
129 /* CONST_PI2extra is negative, so the result of the
130 subtraction can be larger than pi/2. This means
131 that the argument is actually in a different quadrant.
132 The correction is always < pi/2, so it can't overflow
133 into yet another quadrant. */
135 FPU_sub(REV
| LOADED
| TAG_Valid
,
136 (int)&CONST_PI2
, FULL_PRECISION
);
141 #endif /* BETTER_THAN_486 */
143 FPU_settag0(st0_tag
);
144 control_word
= old_cw
;
145 partial_status
= saved_status
& ~SW_C2
; /* Reduction complete. */
147 return (q
& 3) | even
;
150 /* Convert a long to register */
151 static void convert_l2reg(long const *arg
, int deststnr
)
156 FPU_REG
*dest
= &st(deststnr
);
159 FPU_copy_to_regi(&CONST_Z
, TAG_Zero
, deststnr
);
172 setexponent16(dest
, 31);
173 tag
= FPU_normalize(dest
);
174 FPU_settagi(deststnr
, tag
);
179 static void single_arg_error(FPU_REG
*st0_ptr
, u_char st0_tag
)
181 if (st0_tag
== TAG_Empty
)
182 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
183 else if (st0_tag
== TW_NaN
)
184 real_1op_NaN(st0_ptr
); /* return with a NaN in st(0) */
187 EXCEPTION(EX_INTERNAL
| 0x0112);
188 #endif /* PARANOID */
191 static void single_arg_2_error(FPU_REG
*st0_ptr
, u_char st0_tag
)
197 isNaN
= (exponent(st0_ptr
) == EXP_OVER
)
198 && (st0_ptr
->sigh
& 0x80000000);
199 if (isNaN
&& !(st0_ptr
->sigh
& 0x40000000)) { /* Signaling ? */
200 EXCEPTION(EX_Invalid
);
201 if (control_word
& CW_Invalid
) {
202 /* The masked response */
203 /* Convert to a QNaN */
204 st0_ptr
->sigh
|= 0x40000000;
206 FPU_copy_to_reg0(st0_ptr
, TAG_Special
);
211 FPU_copy_to_reg0(st0_ptr
, TAG_Special
);
213 /* pseudoNaN or other unsupported */
214 EXCEPTION(EX_Invalid
);
215 if (control_word
& CW_Invalid
) {
216 /* The masked response */
217 FPU_copy_to_reg0(&CONST_QNaN
, TAG_Special
);
219 FPU_copy_to_reg0(&CONST_QNaN
, TAG_Special
);
222 break; /* return with a NaN in st(0) */
225 EXCEPTION(EX_INTERNAL
| 0x0112);
226 #endif /* PARANOID */
230 /*---------------------------------------------------------------------------*/
232 static void f2xm1(FPU_REG
*st0_ptr
, u_char tag
)
238 if (tag
== TAG_Valid
) {
239 /* For an 80486 FPU, the result is undefined if the arg is >= 1.0 */
240 if (exponent(st0_ptr
) < 0) {
243 FPU_to_exp16(st0_ptr
, &a
);
245 /* poly_2xm1(x) requires 0 < st(0) < 1. */
246 poly_2xm1(getsign(st0_ptr
), &a
, st0_ptr
);
248 set_precision_flag_up(); /* 80486 appears to always do this */
255 if (tag
== TAG_Special
)
256 tag
= FPU_Special(st0_ptr
);
260 if (denormal_operand() < 0)
264 if (signnegative(st0_ptr
)) {
265 /* -infinity gives -1 (p16-10) */
266 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
267 setnegative(st0_ptr
);
271 single_arg_error(st0_ptr
, tag
);
275 static void fptan(FPU_REG
*st0_ptr
, u_char st0_tag
)
279 u_char arg_sign
= getsign(st0_ptr
);
281 /* Stack underflow has higher priority */
282 if (st0_tag
== TAG_Empty
) {
283 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
284 if (control_word
& CW_Invalid
) {
285 st_new_ptr
= &st(-1);
287 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
292 if (STACK_OVERFLOW
) {
293 FPU_stack_overflow();
297 if (st0_tag
== TAG_Valid
) {
298 if (exponent(st0_ptr
) > -40) {
299 if ((q
= trig_arg(st0_ptr
, 0)) == -1) {
300 /* Operand is out of range */
305 setsign(st0_ptr
, (q
& 1) ^ (arg_sign
!= 0));
306 set_precision_flag_up(); /* We do not really know if up or down */
308 /* For a small arg, the result == the argument */
309 /* Underflow may happen */
313 FPU_to_exp16(st0_ptr
, st0_ptr
);
316 FPU_round(st0_ptr
, 1, 0, FULL_PRECISION
, arg_sign
);
317 FPU_settag0(st0_tag
);
320 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
324 if (st0_tag
== TAG_Zero
) {
326 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
331 if (st0_tag
== TAG_Special
)
332 st0_tag
= FPU_Special(st0_ptr
);
334 if (st0_tag
== TW_Denormal
) {
335 if (denormal_operand() < 0)
341 if (st0_tag
== TW_Infinity
) {
342 /* The 80486 treats infinity as an invalid operand */
343 if (arith_invalid(0) >= 0) {
344 st_new_ptr
= &st(-1);
351 single_arg_2_error(st0_ptr
, st0_tag
);
354 static void fxtract(FPU_REG
*st0_ptr
, u_char st0_tag
)
358 register FPU_REG
*st1_ptr
= st0_ptr
; /* anticipate */
360 if (STACK_OVERFLOW
) {
361 FPU_stack_overflow();
367 if (st0_tag
== TAG_Valid
) {
371 sign
= getsign(st1_ptr
);
372 reg_copy(st1_ptr
, st_new_ptr
);
373 setexponent16(st_new_ptr
, exponent(st_new_ptr
));
377 e
= exponent16(st_new_ptr
);
378 convert_l2reg(&e
, 1);
379 setexponentpos(st_new_ptr
, 0);
380 setsign(st_new_ptr
, sign
);
381 FPU_settag0(TAG_Valid
); /* Needed if arg was a denormal */
383 } else if (st0_tag
== TAG_Zero
) {
384 sign
= getsign(st0_ptr
);
386 if (FPU_divide_by_zero(0, SIGN_NEG
) < 0)
390 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
391 setsign(st_new_ptr
, sign
);
395 if (st0_tag
== TAG_Special
)
396 st0_tag
= FPU_Special(st0_ptr
);
398 if (st0_tag
== TW_Denormal
) {
399 if (denormal_operand() < 0)
403 sign
= getsign(st1_ptr
);
404 FPU_to_exp16(st1_ptr
, st_new_ptr
);
406 } else if (st0_tag
== TW_Infinity
) {
407 sign
= getsign(st0_ptr
);
408 setpositive(st0_ptr
);
410 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
411 setsign(st_new_ptr
, sign
);
413 } else if (st0_tag
== TW_NaN
) {
414 if (real_1op_NaN(st0_ptr
) < 0)
418 FPU_copy_to_reg0(st0_ptr
, TAG_Special
);
420 } else if (st0_tag
== TAG_Empty
) {
421 /* Is this the correct behaviour? */
422 if (control_word
& EX_Invalid
) {
423 FPU_stack_underflow();
425 FPU_stack_underflow();
427 EXCEPTION(EX_StackUnder
);
431 EXCEPTION(EX_INTERNAL
| 0x119);
432 #endif /* PARANOID */
435 static void fdecstp(void)
441 static void fincstp(void)
447 static void fsqrt_(FPU_REG
*st0_ptr
, u_char st0_tag
)
453 if (st0_tag
== TAG_Valid
) {
456 if (signnegative(st0_ptr
)) {
457 arith_invalid(0); /* sqrt(negative) is invalid */
461 /* make st(0) in [1.0 .. 4.0) */
462 expon
= exponent(st0_ptr
);
466 setexponent16(st0_ptr
, (expon
& 1));
468 /* Do the computation, the sign of the result will be positive. */
469 tag
= wm_sqrt(st0_ptr
, 0, 0, control_word
, SIGN_POS
);
470 addexponent(st0_ptr
, expon
>> 1);
475 if (st0_tag
== TAG_Zero
)
478 if (st0_tag
== TAG_Special
)
479 st0_tag
= FPU_Special(st0_ptr
);
481 if (st0_tag
== TW_Infinity
) {
482 if (signnegative(st0_ptr
))
483 arith_invalid(0); /* sqrt(-Infinity) is invalid */
485 } else if (st0_tag
== TW_Denormal
) {
486 if (signnegative(st0_ptr
)) {
487 arith_invalid(0); /* sqrt(negative) is invalid */
491 if (denormal_operand() < 0)
494 FPU_to_exp16(st0_ptr
, st0_ptr
);
496 expon
= exponent16(st0_ptr
);
501 single_arg_error(st0_ptr
, st0_tag
);
505 static void frndint_(FPU_REG
*st0_ptr
, u_char st0_tag
)
509 if (st0_tag
== TAG_Valid
) {
514 sign
= getsign(st0_ptr
);
516 if (exponent(st0_ptr
) > 63)
519 if (st0_tag
== TW_Denormal
) {
520 if (denormal_operand() < 0)
524 /* Fortunately, this can't overflow to 2^64 */
525 if ((flags
= FPU_round_to_int(st0_ptr
, st0_tag
)))
526 set_precision_flag(flags
);
528 setexponent16(st0_ptr
, 63);
529 tag
= FPU_normalize(st0_ptr
);
530 setsign(st0_ptr
, sign
);
535 if (st0_tag
== TAG_Zero
)
538 if (st0_tag
== TAG_Special
)
539 st0_tag
= FPU_Special(st0_ptr
);
541 if (st0_tag
== TW_Denormal
)
543 else if (st0_tag
== TW_Infinity
)
546 single_arg_error(st0_ptr
, st0_tag
);
549 static int fsin(FPU_REG
*st0_ptr
, u_char tag
)
551 u_char arg_sign
= getsign(st0_ptr
);
553 if (tag
== TAG_Valid
) {
556 if (exponent(st0_ptr
) > -40) {
557 if ((q
= trig_arg(st0_ptr
, 0)) == -1) {
558 /* Operand is out of range */
567 setsign(st0_ptr
, getsign(st0_ptr
) ^ arg_sign
);
569 /* We do not really know if up or down */
570 set_precision_flag_up();
573 /* For a small arg, the result == the argument */
574 set_precision_flag_up(); /* Must be up. */
579 if (tag
== TAG_Zero
) {
584 if (tag
== TAG_Special
)
585 tag
= FPU_Special(st0_ptr
);
587 if (tag
== TW_Denormal
) {
588 if (denormal_operand() < 0)
591 /* For a small arg, the result == the argument */
592 /* Underflow may happen */
593 FPU_to_exp16(st0_ptr
, st0_ptr
);
595 tag
= FPU_round(st0_ptr
, 1, 0, FULL_PRECISION
, arg_sign
);
600 } else if (tag
== TW_Infinity
) {
601 /* The 80486 treats infinity as an invalid operand */
605 single_arg_error(st0_ptr
, tag
);
610 static int f_cos(FPU_REG
*st0_ptr
, u_char tag
)
614 st0_sign
= getsign(st0_ptr
);
616 if (tag
== TAG_Valid
) {
619 if (exponent(st0_ptr
) > -40) {
620 if ((exponent(st0_ptr
) < 0)
621 || ((exponent(st0_ptr
) == 0)
622 && (significand(st0_ptr
) <=
623 0xc90fdaa22168c234LL
))) {
626 /* We do not really know if up or down */
627 set_precision_flag_down();
630 } else if ((q
= trig_arg(st0_ptr
, FCOS
)) != -1) {
636 /* We do not really know if up or down */
637 set_precision_flag_down();
641 /* Operand is out of range */
648 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
650 set_precision_flag_down(); /* 80486 appears to do this. */
652 set_precision_flag_up(); /* Must be up. */
653 #endif /* PECULIAR_486 */
656 } else if (tag
== TAG_Zero
) {
657 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
662 if (tag
== TAG_Special
)
663 tag
= FPU_Special(st0_ptr
);
665 if (tag
== TW_Denormal
) {
666 if (denormal_operand() < 0)
670 } else if (tag
== TW_Infinity
) {
671 /* The 80486 treats infinity as an invalid operand */
675 single_arg_error(st0_ptr
, tag
); /* requires st0_ptr == &st(0) */
680 static void fcos(FPU_REG
*st0_ptr
, u_char st0_tag
)
682 f_cos(st0_ptr
, st0_tag
);
685 static void fsincos(FPU_REG
*st0_ptr
, u_char st0_tag
)
691 /* Stack underflow has higher priority */
692 if (st0_tag
== TAG_Empty
) {
693 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
694 if (control_word
& CW_Invalid
) {
695 st_new_ptr
= &st(-1);
697 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
702 if (STACK_OVERFLOW
) {
703 FPU_stack_overflow();
707 if (st0_tag
== TAG_Special
)
708 tag
= FPU_Special(st0_ptr
);
713 single_arg_2_error(st0_ptr
, TW_NaN
);
715 } else if (tag
== TW_Infinity
) {
716 /* The 80486 treats infinity as an invalid operand */
717 if (arith_invalid(0) >= 0) {
718 /* Masked response */
725 reg_copy(st0_ptr
, &arg
);
726 if (!fsin(st0_ptr
, st0_tag
)) {
728 FPU_copy_to_reg0(&arg
, st0_tag
);
729 f_cos(&st(0), st0_tag
);
731 /* An error, so restore st(0) */
732 FPU_copy_to_reg0(&arg
, st0_tag
);
736 /*---------------------------------------------------------------------------*/
737 /* The following all require two arguments: st(0) and st(1) */
739 /* A lean, mean kernel for the fprem instructions. This relies upon
740 the division and rounding to an integer in do_fprem giving an
741 exact result. Because of this, rem_kernel() needs to deal only with
742 the least significant 64 bits, the more significant bits of the
745 static void rem_kernel(unsigned long long st0
, unsigned long long *y
,
746 unsigned long long st1
, unsigned long long q
, int n
)
749 unsigned long long x
;
753 /* Do the required multiplication and subtraction in the one operation */
755 /* lsw x -= lsw st1 * lsw q */
756 asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
757 (((unsigned *)&x
)[0]), "=m"(((unsigned *)&x
)[1]),
759 :"2"(((unsigned *)&st1
)[0]), "m"(((unsigned *)&q
)[0])
761 /* msw x -= msw st1 * lsw q */
762 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x
)[1]),
764 :"1"(((unsigned *)&st1
)[1]), "m"(((unsigned *)&q
)[0])
766 /* msw x -= lsw st1 * msw q */
767 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x
)[1]),
769 :"1"(((unsigned *)&st1
)[0]), "m"(((unsigned *)&q
)[1])
775 /* Remainder of st(0) / st(1) */
776 /* This routine produces exact results, i.e. there is never any
777 rounding or truncation, etc of the result. */
778 static void do_fprem(FPU_REG
*st0_ptr
, u_char st0_tag
, int round
)
780 FPU_REG
*st1_ptr
= &st(1);
781 u_char st1_tag
= FPU_gettagi(1);
783 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
784 FPU_REG tmp
, st0
, st1
;
785 u_char st0_sign
, st1_sign
;
791 unsigned short saved_status
;
795 /* Convert registers for internal use. */
796 st0_sign
= FPU_to_exp16(st0_ptr
, &st0
);
797 st1_sign
= FPU_to_exp16(st1_ptr
, &st1
);
798 expdif
= exponent16(&st0
) - exponent16(&st1
);
800 old_cw
= control_word
;
803 /* We want the status following the denorm tests, but don't want
804 the status changed by the arithmetic operations. */
805 saved_status
= partial_status
;
806 control_word
&= ~CW_RC
;
807 control_word
|= RC_CHOP
;
810 /* This should be the most common case */
813 u_char sign
= st0_sign
^ st1_sign
;
814 tag
= FPU_u_div(&st0
, &st1
, &tmp
,
815 PR_64_BITS
| RC_CHOP
| 0x3f,
819 if (exponent(&tmp
) >= 0) {
820 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't
822 q
= significand(&tmp
);
824 rem_kernel(significand(&st0
),
829 setexponent16(&tmp
, exponent16(&st1
));
831 reg_copy(&st0
, &tmp
);
835 if ((round
== RC_RND
)
836 && (tmp
.sigh
& 0xc0000000)) {
837 /* We may need to subtract st(1) once more,
838 to get a result <= 1/2 of st(1). */
839 unsigned long long x
;
841 exponent16(&st1
) - exponent16(&tmp
);
844 x
= significand(&st1
) -
846 else /* expdif is 1 */
847 x
= (significand(&st1
)
850 if ((x
< significand(&tmp
)) ||
851 /* or equi-distant (from 0 & st(1)) and q is odd */
852 ((x
== significand(&tmp
))
854 st0_sign
= !st0_sign
;
855 significand(&tmp
) = x
;
868 control_word
= old_cw
;
873 /* There is a large exponent difference ( >= 64 ) */
874 /* To make much sense, the code in this section should
875 be done at high precision. */
879 /* prevent overflow here */
880 /* N is 'a number between 32 and 63' (p26-113) */
881 reg_copy(&st0
, &tmp
);
883 N
= (expdif
& 0x0000001f) + 32; /* This choice gives results
884 identical to an AMD 486 */
885 setexponent16(&tmp
, N
);
886 exp_1
= exponent16(&st1
);
887 setexponent16(&st1
, 0);
890 sign
= getsign(&tmp
) ^ st1_sign
;
892 FPU_u_div(&tmp
, &st1
, &tmp
,
893 PR_64_BITS
| RC_CHOP
| 0x3f, sign
);
896 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't
899 rem_kernel(significand(&st0
),
902 significand(&tmp
), exponent(&tmp
)
904 setexponent16(&tmp
, exp_1
+ expdif
);
906 /* It is possible for the operation to be complete here.
907 What does the IEEE standard say? The Intel 80486 manual
908 implies that the operation will never be completed at this
909 point, and the behaviour of a real 80486 confirms this.
911 if (!(tmp
.sigh
| tmp
.sigl
)) {
912 /* The result is zero */
913 control_word
= old_cw
;
914 partial_status
= saved_status
;
915 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
916 setsign(&st0
, st0_sign
);
921 #endif /* PECULIAR_486 */
927 control_word
= old_cw
;
928 partial_status
= saved_status
;
929 tag
= FPU_normalize_nuo(&tmp
);
930 reg_copy(&tmp
, st0_ptr
);
932 /* The only condition to be looked for is underflow,
933 and it can occur here only if underflow is unmasked. */
934 if ((exponent16(&tmp
) <= EXP_UNDER
) && (tag
!= TAG_Zero
)
935 && !(control_word
& CW_Underflow
)) {
937 tag
= arith_underflow(st0_ptr
);
938 setsign(st0_ptr
, st0_sign
);
941 } else if ((exponent16(&tmp
) > EXP_UNDER
) || (tag
== TAG_Zero
)) {
943 setsign(st0_ptr
, st0_sign
);
946 FPU_round(st0_ptr
, 0, 0, FULL_PRECISION
, st0_sign
);
954 if (st0_tag
== TAG_Special
)
955 st0_tag
= FPU_Special(st0_ptr
);
956 if (st1_tag
== TAG_Special
)
957 st1_tag
= FPU_Special(st1_ptr
);
959 if (((st0_tag
== TAG_Valid
) && (st1_tag
== TW_Denormal
))
960 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TAG_Valid
))
961 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TW_Denormal
))) {
962 if (denormal_operand() < 0)
965 } else if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
966 FPU_stack_underflow();
968 } else if (st0_tag
== TAG_Zero
) {
969 if (st1_tag
== TAG_Valid
) {
972 } else if (st1_tag
== TW_Denormal
) {
973 if (denormal_operand() < 0)
977 } else if (st1_tag
== TAG_Zero
) {
980 } /* fprem(?,0) always invalid */
981 else if (st1_tag
== TW_Infinity
) {
985 } else if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
986 if (st1_tag
== TAG_Zero
) {
987 arith_invalid(0); /* fprem(Valid,Zero) is invalid */
989 } else if (st1_tag
!= TW_NaN
) {
990 if (((st0_tag
== TW_Denormal
)
991 || (st1_tag
== TW_Denormal
))
992 && (denormal_operand() < 0))
995 if (st1_tag
== TW_Infinity
) {
996 /* fprem(Valid,Infinity) is o.k. */
1001 } else if (st0_tag
== TW_Infinity
) {
1002 if (st1_tag
!= TW_NaN
) {
1003 arith_invalid(0); /* fprem(Infinity,?) is invalid */
1008 /* One of the registers must contain a NaN if we got here. */
1011 if ((st0_tag
!= TW_NaN
) && (st1_tag
!= TW_NaN
))
1012 EXCEPTION(EX_INTERNAL
| 0x118);
1013 #endif /* PARANOID */
1015 real_2op_NaN(st1_ptr
, st1_tag
, 0, st1_ptr
);
1019 /* ST(1) <- ST(1) * log ST; pop ST */
1020 static void fyl2x(FPU_REG
*st0_ptr
, u_char st0_tag
)
1022 FPU_REG
*st1_ptr
= &st(1), exponent
;
1023 u_char st1_tag
= FPU_gettagi(1);
1029 if ((st0_tag
== TAG_Valid
) && (st1_tag
== TAG_Valid
)) {
1031 /* Both regs are Valid or Denormal */
1032 if (signpositive(st0_ptr
)) {
1033 if (st0_tag
== TW_Denormal
)
1034 FPU_to_exp16(st0_ptr
, st0_ptr
);
1036 /* Convert st(0) for internal use. */
1037 setexponent16(st0_ptr
, exponent(st0_ptr
));
1039 if ((st0_ptr
->sigh
== 0x80000000)
1040 && (st0_ptr
->sigl
== 0)) {
1041 /* Special case. The result can be precise. */
1043 e
= exponent16(st0_ptr
);
1052 setexponent16(&exponent
, 31);
1053 tag
= FPU_normalize_nuo(&exponent
);
1055 setsign(&exponent
, esign
);
1057 FPU_mul(&exponent
, tag
, 1, FULL_PRECISION
);
1059 FPU_settagi(1, tag
);
1061 /* The usual case */
1062 sign
= getsign(st1_ptr
);
1063 if (st1_tag
== TW_Denormal
)
1064 FPU_to_exp16(st1_ptr
, st1_ptr
);
1066 /* Convert st(1) for internal use. */
1067 setexponent16(st1_ptr
,
1069 poly_l2(st0_ptr
, st1_ptr
, sign
);
1073 if (arith_invalid(1) < 0)
1082 if (st0_tag
== TAG_Special
)
1083 st0_tag
= FPU_Special(st0_ptr
);
1084 if (st1_tag
== TAG_Special
)
1085 st1_tag
= FPU_Special(st1_ptr
);
1087 if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
1088 FPU_stack_underflow_pop(1);
1090 } else if ((st0_tag
<= TW_Denormal
) && (st1_tag
<= TW_Denormal
)) {
1091 if (st0_tag
== TAG_Zero
) {
1092 if (st1_tag
== TAG_Zero
) {
1093 /* Both args zero is invalid */
1094 if (arith_invalid(1) < 0)
1098 sign
= getsign(st1_ptr
) ^ SIGN_NEG
;
1099 if (FPU_divide_by_zero(1, sign
) < 0)
1102 setsign(st1_ptr
, sign
);
1104 } else if (st1_tag
== TAG_Zero
) {
1105 /* st(1) contains zero, st(0) valid <> 0 */
1106 /* Zero is the valid answer */
1107 sign
= getsign(st1_ptr
);
1109 if (signnegative(st0_ptr
)) {
1111 if (arith_invalid(1) < 0)
1113 } else if ((st0_tag
== TW_Denormal
)
1114 && (denormal_operand() < 0))
1117 if (exponent(st0_ptr
) < 0)
1120 FPU_copy_to_reg1(&CONST_Z
, TAG_Zero
);
1121 setsign(st1_ptr
, sign
);
1124 /* One or both operands are denormals. */
1125 if (denormal_operand() < 0)
1129 } else if ((st0_tag
== TW_NaN
) || (st1_tag
== TW_NaN
)) {
1130 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1133 /* One or both arg must be an infinity */
1134 else if (st0_tag
== TW_Infinity
) {
1135 if ((signnegative(st0_ptr
)) || (st1_tag
== TAG_Zero
)) {
1136 /* log(-infinity) or 0*log(infinity) */
1137 if (arith_invalid(1) < 0)
1140 u_char sign
= getsign(st1_ptr
);
1142 if ((st1_tag
== TW_Denormal
)
1143 && (denormal_operand() < 0))
1146 FPU_copy_to_reg1(&CONST_INF
, TAG_Special
);
1147 setsign(st1_ptr
, sign
);
1150 /* st(1) must be infinity here */
1151 else if (((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
))
1152 && (signpositive(st0_ptr
))) {
1153 if (exponent(st0_ptr
) >= 0) {
1154 if ((exponent(st0_ptr
) == 0) &&
1155 (st0_ptr
->sigh
== 0x80000000) &&
1156 (st0_ptr
->sigl
== 0)) {
1157 /* st(0) holds 1.0 */
1158 /* infinity*log(1) */
1159 if (arith_invalid(1) < 0)
1162 /* else st(0) is positive and > 1.0 */
1164 /* st(0) is positive and < 1.0 */
1166 if ((st0_tag
== TW_Denormal
)
1167 && (denormal_operand() < 0))
1170 changesign(st1_ptr
);
1173 /* st(0) must be zero or negative */
1174 if (st0_tag
== TAG_Zero
) {
1175 /* This should be invalid, but a real 80486 is happy with it. */
1177 #ifndef PECULIAR_486
1178 sign
= getsign(st1_ptr
);
1179 if (FPU_divide_by_zero(1, sign
) < 0)
1181 #endif /* PECULIAR_486 */
1183 changesign(st1_ptr
);
1184 } else if (arith_invalid(1) < 0) /* log(negative) */
1191 static void fpatan(FPU_REG
*st0_ptr
, u_char st0_tag
)
1193 FPU_REG
*st1_ptr
= &st(1);
1194 u_char st1_tag
= FPU_gettagi(1);
1198 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1201 poly_atan(st0_ptr
, st0_tag
, st1_ptr
, st1_tag
);
1208 if (st0_tag
== TAG_Special
)
1209 st0_tag
= FPU_Special(st0_ptr
);
1210 if (st1_tag
== TAG_Special
)
1211 st1_tag
= FPU_Special(st1_ptr
);
1213 if (((st0_tag
== TAG_Valid
) && (st1_tag
== TW_Denormal
))
1214 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TAG_Valid
))
1215 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TW_Denormal
))) {
1216 if (denormal_operand() < 0)
1220 } else if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
1221 FPU_stack_underflow_pop(1);
1223 } else if ((st0_tag
== TW_NaN
) || (st1_tag
== TW_NaN
)) {
1224 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) >= 0)
1227 } else if ((st0_tag
== TW_Infinity
) || (st1_tag
== TW_Infinity
)) {
1228 u_char sign
= getsign(st1_ptr
);
1229 if (st0_tag
== TW_Infinity
) {
1230 if (st1_tag
== TW_Infinity
) {
1231 if (signpositive(st0_ptr
)) {
1232 FPU_copy_to_reg1(&CONST_PI4
, TAG_Valid
);
1234 setpositive(st1_ptr
);
1236 FPU_u_add(&CONST_PI4
, &CONST_PI2
,
1237 st1_ptr
, FULL_PRECISION
,
1239 exponent(&CONST_PI4
),
1240 exponent(&CONST_PI2
));
1242 FPU_settagi(1, tag
);
1245 if ((st1_tag
== TW_Denormal
)
1246 && (denormal_operand() < 0))
1249 if (signpositive(st0_ptr
)) {
1250 FPU_copy_to_reg1(&CONST_Z
, TAG_Zero
);
1251 setsign(st1_ptr
, sign
); /* An 80486 preserves the sign */
1255 FPU_copy_to_reg1(&CONST_PI
, TAG_Valid
);
1259 /* st(1) is infinity, st(0) not infinity */
1260 if ((st0_tag
== TW_Denormal
)
1261 && (denormal_operand() < 0))
1264 FPU_copy_to_reg1(&CONST_PI2
, TAG_Valid
);
1266 setsign(st1_ptr
, sign
);
1267 } else if (st1_tag
== TAG_Zero
) {
1268 /* st(0) must be valid or zero */
1269 u_char sign
= getsign(st1_ptr
);
1271 if ((st0_tag
== TW_Denormal
) && (denormal_operand() < 0))
1274 if (signpositive(st0_ptr
)) {
1275 /* An 80486 preserves the sign */
1280 FPU_copy_to_reg1(&CONST_PI
, TAG_Valid
);
1281 setsign(st1_ptr
, sign
);
1282 } else if (st0_tag
== TAG_Zero
) {
1283 /* st(1) must be TAG_Valid here */
1284 u_char sign
= getsign(st1_ptr
);
1286 if ((st1_tag
== TW_Denormal
) && (denormal_operand() < 0))
1289 FPU_copy_to_reg1(&CONST_PI2
, TAG_Valid
);
1290 setsign(st1_ptr
, sign
);
1294 EXCEPTION(EX_INTERNAL
| 0x125);
1295 #endif /* PARANOID */
1298 set_precision_flag_up(); /* We do not really know if up or down */
1301 static void fprem(FPU_REG
*st0_ptr
, u_char st0_tag
)
1303 do_fprem(st0_ptr
, st0_tag
, RC_CHOP
);
1306 static void fprem1(FPU_REG
*st0_ptr
, u_char st0_tag
)
1308 do_fprem(st0_ptr
, st0_tag
, RC_RND
);
1311 static void fyl2xp1(FPU_REG
*st0_ptr
, u_char st0_tag
)
1314 FPU_REG
*st1_ptr
= &st(1), a
, b
;
1315 u_char st1_tag
= FPU_gettagi(1);
1318 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1321 sign
= getsign(st0_ptr
);
1322 sign1
= getsign(st1_ptr
);
1324 FPU_to_exp16(st0_ptr
, &a
);
1325 FPU_to_exp16(st1_ptr
, &b
);
1327 if (poly_l2p1(sign
, sign1
, &a
, &b
, st1_ptr
))
1334 if (st0_tag
== TAG_Special
)
1335 st0_tag
= FPU_Special(st0_ptr
);
1336 if (st1_tag
== TAG_Special
)
1337 st1_tag
= FPU_Special(st1_ptr
);
1339 if (((st0_tag
== TAG_Valid
) && (st1_tag
== TW_Denormal
))
1340 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TAG_Valid
))
1341 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TW_Denormal
))) {
1342 if (denormal_operand() < 0)
1346 } else if ((st0_tag
== TAG_Empty
) | (st1_tag
== TAG_Empty
)) {
1347 FPU_stack_underflow_pop(1);
1349 } else if (st0_tag
== TAG_Zero
) {
1352 if (denormal_operand() < 0)
1357 setsign(st0_ptr
, getsign(st0_ptr
) ^ getsign(st1_ptr
));
1358 FPU_copy_to_reg1(st0_ptr
, st0_tag
);
1362 /* Infinity*log(1) */
1363 if (arith_invalid(1) < 0)
1368 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1374 EXCEPTION(EX_INTERNAL
| 0x116);
1376 #endif /* PARANOID */
1379 } else if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
1382 if (signnegative(st0_ptr
)) {
1383 if (exponent(st0_ptr
) >= 0) {
1384 /* st(0) holds <= -1.0 */
1385 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1386 changesign(st1_ptr
);
1388 if (arith_invalid(1) < 0)
1390 #endif /* PECULIAR_486 */
1391 } else if ((st0_tag
== TW_Denormal
)
1392 && (denormal_operand() < 0))
1395 changesign(st1_ptr
);
1396 } else if ((st0_tag
== TW_Denormal
)
1397 && (denormal_operand() < 0))
1402 if (signnegative(st0_ptr
)) {
1403 if ((exponent(st0_ptr
) >= 0) &&
1404 !((st0_ptr
->sigh
== 0x80000000) &&
1405 (st0_ptr
->sigl
== 0))) {
1406 /* st(0) holds < -1.0 */
1407 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1408 changesign(st1_ptr
);
1410 if (arith_invalid(1) < 0)
1412 #endif /* PECULIAR_486 */
1413 } else if ((st0_tag
== TW_Denormal
)
1414 && (denormal_operand() < 0))
1417 changesign(st1_ptr
);
1418 } else if ((st0_tag
== TW_Denormal
)
1419 && (denormal_operand() < 0))
1424 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1428 } else if (st0_tag
== TW_NaN
) {
1429 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1431 } else if (st0_tag
== TW_Infinity
) {
1432 if (st1_tag
== TW_NaN
) {
1433 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1435 } else if (signnegative(st0_ptr
)) {
1436 #ifndef PECULIAR_486
1437 /* This should have higher priority than denormals, but... */
1438 if (arith_invalid(1) < 0) /* log(-infinity) */
1440 #endif /* PECULIAR_486 */
1441 if ((st1_tag
== TW_Denormal
)
1442 && (denormal_operand() < 0))
1445 /* Denormal operands actually get higher priority */
1446 if (arith_invalid(1) < 0) /* log(-infinity) */
1448 #endif /* PECULIAR_486 */
1449 } else if (st1_tag
== TAG_Zero
) {
1451 if (arith_invalid(1) < 0)
1455 /* st(1) must be valid here. */
1457 else if ((st1_tag
== TW_Denormal
) && (denormal_operand() < 0))
1460 /* The Manual says that log(Infinity) is invalid, but a real
1461 80486 sensibly says that it is o.k. */
1463 u_char sign
= getsign(st1_ptr
);
1464 FPU_copy_to_reg1(&CONST_INF
, TAG_Special
);
1465 setsign(st1_ptr
, sign
);
1470 EXCEPTION(EX_INTERNAL
| 0x117);
1473 #endif /* PARANOID */
1480 static void fscale(FPU_REG
*st0_ptr
, u_char st0_tag
)
1482 FPU_REG
*st1_ptr
= &st(1);
1483 u_char st1_tag
= FPU_gettagi(1);
1484 int old_cw
= control_word
;
1485 u_char sign
= getsign(st0_ptr
);
1488 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1492 /* Convert register for internal use. */
1493 setexponent16(st0_ptr
, exponent(st0_ptr
));
1497 if (exponent(st1_ptr
) > 30) {
1498 /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1500 if (signpositive(st1_ptr
)) {
1501 EXCEPTION(EX_Overflow
);
1502 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
1504 EXCEPTION(EX_Underflow
);
1505 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
1507 setsign(st0_ptr
, sign
);
1511 control_word
&= ~CW_RC
;
1512 control_word
|= RC_CHOP
;
1513 reg_copy(st1_ptr
, &tmp
);
1514 FPU_round_to_int(&tmp
, st1_tag
); /* This can never overflow here */
1515 control_word
= old_cw
;
1516 scale
= signnegative(st1_ptr
) ? -tmp
.sigl
: tmp
.sigl
;
1517 scale
+= exponent16(st0_ptr
);
1519 setexponent16(st0_ptr
, scale
);
1521 /* Use FPU_round() to properly detect under/overflow etc */
1522 FPU_round(st0_ptr
, 0, 0, control_word
, sign
);
1527 if (st0_tag
== TAG_Special
)
1528 st0_tag
= FPU_Special(st0_ptr
);
1529 if (st1_tag
== TAG_Special
)
1530 st1_tag
= FPU_Special(st1_ptr
);
1532 if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
1535 /* st(0) must be a denormal */
1536 if ((st0_tag
== TW_Denormal
)
1537 && (denormal_operand() < 0))
1540 FPU_to_exp16(st0_ptr
, st0_ptr
); /* Will not be left on stack */
1544 if (st0_tag
== TW_Denormal
)
1553 if ((st0_tag
== TW_Denormal
)
1554 && (denormal_operand() < 0))
1557 if (signpositive(st1_ptr
))
1558 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
1560 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
1561 setsign(st0_ptr
, sign
);
1565 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1568 } else if (st0_tag
== TAG_Zero
) {
1579 if (signpositive(st1_ptr
))
1580 arith_invalid(0); /* Zero scaled by +Infinity */
1584 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1587 } else if (st0_tag
== TW_Infinity
) {
1598 if (signnegative(st1_ptr
))
1599 arith_invalid(0); /* Infinity scaled by -Infinity */
1603 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1606 } else if (st0_tag
== TW_NaN
) {
1607 if (st1_tag
!= TAG_Empty
) {
1608 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1613 if (!((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
))) {
1614 EXCEPTION(EX_INTERNAL
| 0x115);
1619 /* At least one of st(0), st(1) must be empty */
1620 FPU_stack_underflow();
1624 /*---------------------------------------------------------------------------*/
1626 static FUNC_ST0
const trig_table_a
[] = {
1627 f2xm1
, fyl2x
, fptan
, fpatan
,
1628 fxtract
, fprem1
, (FUNC_ST0
) fdecstp
, (FUNC_ST0
) fincstp
1631 void FPU_triga(void)
1633 (trig_table_a
[FPU_rm
]) (&st(0), FPU_gettag0());
1636 static FUNC_ST0
const trig_table_b
[] = {
1637 fprem
, fyl2xp1
, fsqrt_
, fsincos
, frndint_
, fscale
, (FUNC_ST0
) fsin
, fcos
1640 void FPU_trigb(void)
1642 (trig_table_b
[FPU_rm
]) (&st(0), FPU_gettag0());