1 // SPDX-License-Identifier: GPL-2.0
2 /*---------------------------------------------------------------------------+
5 | Implementation of the FPU "transcendental" functions. |
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 |
12 +---------------------------------------------------------------------------*/
14 #include "fpu_system.h"
15 #include "exception.h"
18 #include "control_w.h"
19 #include "reg_constant.h"
21 static void rem_kernel(unsigned long long st0
, unsigned long long *y
,
22 unsigned long long st1
, unsigned long long q
, int n
);
24 #define BETTER_THAN_486
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 */
34 static int trig_arg(FPU_REG
*st0_ptr
, int even
)
39 int old_cw
= control_word
, saved_status
= partial_status
;
40 int tag
, st0_tag
= TAG_Valid
;
42 if (exponent(st0_ptr
) >= 63) {
43 partial_status
|= SW_C2
; /* Reduction incomplete. */
47 control_word
&= ~CW_RC
;
48 control_word
|= RC_CHOP
;
51 tag
= FPU_u_div(st0_ptr
, &CONST_PI2
, &tmp
, PR_64_BITS
| RC_CHOP
| 0x3f,
54 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't overflow
56 q
= significand(&tmp
);
58 rem_kernel(significand(st0_ptr
),
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
);
67 if ((even
&& !(q
& 1)) || (!even
&& (q
& 1))) {
69 FPU_sub(REV
| LOADED
| TAG_Valid
, (int)&CONST_PI2
,
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)
79 /* This code gives the effect of having pi/2 to better than
80 128 bits precision. */
82 significand(&tmp
) = q
+ 1;
83 setexponent16(&tmp
, 63);
86 FPU_u_mul(&CONST_PI2extra
, &tmp
, &tmp
,
87 FULL_PRECISION
, SIGN_POS
,
88 exponent(&CONST_PI2extra
) +
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. */
101 #endif /* BETTER_THAN_486 */
103 #ifdef BETTER_THAN_486
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. */
110 && (exponent(st0_ptr
) <= exponent(&CONST_PI2extra
) + 64))
112 /* This code gives the effect of having p/2 to better than
113 128 bits precision. */
115 significand(&tmp
) = q
;
116 setexponent16(&tmp
, 63);
117 FPU_normalize(&tmp
); /* This must return TAG_Valid */
119 FPU_u_mul(&CONST_PI2extra
, &tmp
, &tmp
,
120 FULL_PRECISION
, SIGN_POS
,
121 exponent(&CONST_PI2extra
) +
123 setsign(&tmp
, getsign(&CONST_PI2extra
));
124 st0_tag
= FPU_sub(LOADED
| (tmptag
& 0x0f), (int)&tmp
,
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. */
136 FPU_sub(REV
| LOADED
| TAG_Valid
,
137 (int)&CONST_PI2
, FULL_PRECISION
);
142 #endif /* BETTER_THAN_486 */
144 FPU_settag0(st0_tag
);
145 control_word
= old_cw
;
146 partial_status
= saved_status
& ~SW_C2
; /* Reduction complete. */
148 return (q
& 3) | even
;
151 /* Convert a long to register */
152 static void convert_l2reg(long const *arg
, int deststnr
)
157 FPU_REG
*dest
= &st(deststnr
);
160 FPU_copy_to_regi(&CONST_Z
, TAG_Zero
, deststnr
);
173 setexponent16(dest
, 31);
174 tag
= FPU_normalize(dest
);
175 FPU_settagi(deststnr
, tag
);
180 static void single_arg_error(FPU_REG
*st0_ptr
, u_char st0_tag
)
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) */
188 EXCEPTION(EX_INTERNAL
| 0x0112);
189 #endif /* PARANOID */
192 static void single_arg_2_error(FPU_REG
*st0_ptr
, u_char st0_tag
)
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;
207 FPU_copy_to_reg0(st0_ptr
, TAG_Special
);
212 FPU_copy_to_reg0(st0_ptr
, TAG_Special
);
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
);
220 FPU_copy_to_reg0(&CONST_QNaN
, TAG_Special
);
223 break; /* return with a NaN in st(0) */
226 EXCEPTION(EX_INTERNAL
| 0x0112);
227 #endif /* PARANOID */
231 /*---------------------------------------------------------------------------*/
233 static void f2xm1(FPU_REG
*st0_ptr
, u_char tag
)
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) {
244 FPU_to_exp16(st0_ptr
, &a
);
246 /* poly_2xm1(x) requires 0 < st(0) < 1. */
247 poly_2xm1(getsign(st0_ptr
), &a
, st0_ptr
);
249 set_precision_flag_up(); /* 80486 appears to always do this */
256 if (tag
== TAG_Special
)
257 tag
= FPU_Special(st0_ptr
);
261 if (denormal_operand() < 0)
265 if (signnegative(st0_ptr
)) {
266 /* -infinity gives -1 (p16-10) */
267 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
268 setnegative(st0_ptr
);
272 single_arg_error(st0_ptr
, tag
);
276 static void fptan(FPU_REG
*st0_ptr
, u_char st0_tag
)
280 u_char arg_sign
= getsign(st0_ptr
);
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);
288 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
293 if (STACK_OVERFLOW
) {
294 FPU_stack_overflow();
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 */
306 setsign(st0_ptr
, (q
& 1) ^ (arg_sign
!= 0));
307 set_precision_flag_up(); /* We do not really know if up or down */
309 /* For a small arg, the result == the argument */
310 /* Underflow may happen */
314 FPU_to_exp16(st0_ptr
, st0_ptr
);
317 FPU_round(st0_ptr
, 1, 0, FULL_PRECISION
, arg_sign
);
318 FPU_settag0(st0_tag
);
321 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
325 if (st0_tag
== TAG_Zero
) {
327 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
332 if (st0_tag
== TAG_Special
)
333 st0_tag
= FPU_Special(st0_ptr
);
335 if (st0_tag
== TW_Denormal
) {
336 if (denormal_operand() < 0)
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);
352 single_arg_2_error(st0_ptr
, st0_tag
);
355 static void fxtract(FPU_REG
*st0_ptr
, u_char st0_tag
)
359 register FPU_REG
*st1_ptr
= st0_ptr
; /* anticipate */
361 if (STACK_OVERFLOW
) {
362 FPU_stack_overflow();
368 if (st0_tag
== TAG_Valid
) {
372 sign
= getsign(st1_ptr
);
373 reg_copy(st1_ptr
, st_new_ptr
);
374 setexponent16(st_new_ptr
, exponent(st_new_ptr
));
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 */
384 } else if (st0_tag
== TAG_Zero
) {
385 sign
= getsign(st0_ptr
);
387 if (FPU_divide_by_zero(0, SIGN_NEG
) < 0)
391 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
392 setsign(st_new_ptr
, sign
);
396 if (st0_tag
== TAG_Special
)
397 st0_tag
= FPU_Special(st0_ptr
);
399 if (st0_tag
== TW_Denormal
) {
400 if (denormal_operand() < 0)
404 sign
= getsign(st1_ptr
);
405 FPU_to_exp16(st1_ptr
, st_new_ptr
);
407 } else if (st0_tag
== TW_Infinity
) {
408 sign
= getsign(st0_ptr
);
409 setpositive(st0_ptr
);
411 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
412 setsign(st_new_ptr
, sign
);
414 } else if (st0_tag
== TW_NaN
) {
415 if (real_1op_NaN(st0_ptr
) < 0)
419 FPU_copy_to_reg0(st0_ptr
, TAG_Special
);
421 } else if (st0_tag
== TAG_Empty
) {
422 /* Is this the correct behaviour? */
423 if (control_word
& EX_Invalid
) {
424 FPU_stack_underflow();
426 FPU_stack_underflow();
428 EXCEPTION(EX_StackUnder
);
432 EXCEPTION(EX_INTERNAL
| 0x119);
433 #endif /* PARANOID */
436 static void fdecstp(void)
442 static void fincstp(void)
448 static void fsqrt_(FPU_REG
*st0_ptr
, u_char st0_tag
)
454 if (st0_tag
== TAG_Valid
) {
457 if (signnegative(st0_ptr
)) {
458 arith_invalid(0); /* sqrt(negative) is invalid */
462 /* make st(0) in [1.0 .. 4.0) */
463 expon
= exponent(st0_ptr
);
467 setexponent16(st0_ptr
, (expon
& 1));
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);
476 if (st0_tag
== TAG_Zero
)
479 if (st0_tag
== TAG_Special
)
480 st0_tag
= FPU_Special(st0_ptr
);
482 if (st0_tag
== TW_Infinity
) {
483 if (signnegative(st0_ptr
))
484 arith_invalid(0); /* sqrt(-Infinity) is invalid */
486 } else if (st0_tag
== TW_Denormal
) {
487 if (signnegative(st0_ptr
)) {
488 arith_invalid(0); /* sqrt(negative) is invalid */
492 if (denormal_operand() < 0)
495 FPU_to_exp16(st0_ptr
, st0_ptr
);
497 expon
= exponent16(st0_ptr
);
502 single_arg_error(st0_ptr
, st0_tag
);
506 static void frndint_(FPU_REG
*st0_ptr
, u_char st0_tag
)
510 if (st0_tag
== TAG_Valid
) {
515 sign
= getsign(st0_ptr
);
517 if (exponent(st0_ptr
) > 63)
520 if (st0_tag
== TW_Denormal
) {
521 if (denormal_operand() < 0)
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
);
529 setexponent16(st0_ptr
, 63);
530 tag
= FPU_normalize(st0_ptr
);
531 setsign(st0_ptr
, sign
);
536 if (st0_tag
== TAG_Zero
)
539 if (st0_tag
== TAG_Special
)
540 st0_tag
= FPU_Special(st0_ptr
);
542 if (st0_tag
== TW_Denormal
)
544 else if (st0_tag
== TW_Infinity
)
547 single_arg_error(st0_ptr
, st0_tag
);
550 static int fsin(FPU_REG
*st0_ptr
, u_char tag
)
552 u_char arg_sign
= getsign(st0_ptr
);
554 if (tag
== TAG_Valid
) {
557 if (exponent(st0_ptr
) > -40) {
558 if ((q
= trig_arg(st0_ptr
, 0)) == -1) {
559 /* Operand is out of range */
568 setsign(st0_ptr
, getsign(st0_ptr
) ^ arg_sign
);
570 /* We do not really know if up or down */
571 set_precision_flag_up();
574 /* For a small arg, the result == the argument */
575 set_precision_flag_up(); /* Must be up. */
580 if (tag
== TAG_Zero
) {
585 if (tag
== TAG_Special
)
586 tag
= FPU_Special(st0_ptr
);
588 if (tag
== TW_Denormal
) {
589 if (denormal_operand() < 0)
592 /* For a small arg, the result == the argument */
593 /* Underflow may happen */
594 FPU_to_exp16(st0_ptr
, st0_ptr
);
596 tag
= FPU_round(st0_ptr
, 1, 0, FULL_PRECISION
, arg_sign
);
601 } else if (tag
== TW_Infinity
) {
602 /* The 80486 treats infinity as an invalid operand */
606 single_arg_error(st0_ptr
, tag
);
611 static int f_cos(FPU_REG
*st0_ptr
, u_char tag
)
615 st0_sign
= getsign(st0_ptr
);
617 if (tag
== TAG_Valid
) {
620 if (exponent(st0_ptr
) > -40) {
621 if ((exponent(st0_ptr
) < 0)
622 || ((exponent(st0_ptr
) == 0)
623 && (significand(st0_ptr
) <=
624 0xc90fdaa22168c234LL
))) {
627 /* We do not really know if up or down */
628 set_precision_flag_down();
631 } else if ((q
= trig_arg(st0_ptr
, FCOS
)) != -1) {
637 /* We do not really know if up or down */
638 set_precision_flag_down();
642 /* Operand is out of range */
649 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
651 set_precision_flag_down(); /* 80486 appears to do this. */
653 set_precision_flag_up(); /* Must be up. */
654 #endif /* PECULIAR_486 */
657 } else if (tag
== TAG_Zero
) {
658 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
663 if (tag
== TAG_Special
)
664 tag
= FPU_Special(st0_ptr
);
666 if (tag
== TW_Denormal
) {
667 if (denormal_operand() < 0)
671 } else if (tag
== TW_Infinity
) {
672 /* The 80486 treats infinity as an invalid operand */
676 single_arg_error(st0_ptr
, tag
); /* requires st0_ptr == &st(0) */
681 static void fcos(FPU_REG
*st0_ptr
, u_char st0_tag
)
683 f_cos(st0_ptr
, st0_tag
);
686 static void fsincos(FPU_REG
*st0_ptr
, u_char st0_tag
)
692 /* Stack underflow has higher priority */
693 if (st0_tag
== TAG_Empty
) {
694 FPU_stack_underflow(); /* Puts a QNaN in st(0) */
695 if (control_word
& CW_Invalid
) {
696 st_new_ptr
= &st(-1);
698 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
703 if (STACK_OVERFLOW
) {
704 FPU_stack_overflow();
708 if (st0_tag
== TAG_Special
)
709 tag
= FPU_Special(st0_ptr
);
714 single_arg_2_error(st0_ptr
, TW_NaN
);
716 } else if (tag
== TW_Infinity
) {
717 /* The 80486 treats infinity as an invalid operand */
718 if (arith_invalid(0) >= 0) {
719 /* Masked response */
726 reg_copy(st0_ptr
, &arg
);
727 if (!fsin(st0_ptr
, st0_tag
)) {
729 FPU_copy_to_reg0(&arg
, st0_tag
);
730 f_cos(&st(0), st0_tag
);
732 /* An error, so restore st(0) */
733 FPU_copy_to_reg0(&arg
, st0_tag
);
737 /*---------------------------------------------------------------------------*/
738 /* The following all require two arguments: st(0) and st(1) */
740 /* A lean, mean kernel for the fprem instructions. This relies upon
741 the division and rounding to an integer in do_fprem giving an
742 exact result. Because of this, rem_kernel() needs to deal only with
743 the least significant 64 bits, the more significant bits of the
746 static void rem_kernel(unsigned long long st0
, unsigned long long *y
,
747 unsigned long long st1
, unsigned long long q
, int n
)
750 unsigned long long x
;
754 /* Do the required multiplication and subtraction in the one operation */
756 /* lsw x -= lsw st1 * lsw q */
757 asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1":"=m"
758 (((unsigned *)&x
)[0]), "=m"(((unsigned *)&x
)[1]),
760 :"2"(((unsigned *)&st1
)[0]), "m"(((unsigned *)&q
)[0])
762 /* msw x -= msw st1 * lsw q */
763 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x
)[1]),
765 :"1"(((unsigned *)&st1
)[1]), "m"(((unsigned *)&q
)[0])
767 /* msw x -= lsw st1 * msw q */
768 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x
)[1]),
770 :"1"(((unsigned *)&st1
)[0]), "m"(((unsigned *)&q
)[1])
776 /* Remainder of st(0) / st(1) */
777 /* This routine produces exact results, i.e. there is never any
778 rounding or truncation, etc of the result. */
779 static void do_fprem(FPU_REG
*st0_ptr
, u_char st0_tag
, int round
)
781 FPU_REG
*st1_ptr
= &st(1);
782 u_char st1_tag
= FPU_gettagi(1);
784 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
785 FPU_REG tmp
, st0
, st1
;
786 u_char st0_sign
, st1_sign
;
792 unsigned short saved_status
;
796 /* Convert registers for internal use. */
797 st0_sign
= FPU_to_exp16(st0_ptr
, &st0
);
798 st1_sign
= FPU_to_exp16(st1_ptr
, &st1
);
799 expdif
= exponent16(&st0
) - exponent16(&st1
);
801 old_cw
= control_word
;
804 /* We want the status following the denorm tests, but don't want
805 the status changed by the arithmetic operations. */
806 saved_status
= partial_status
;
807 control_word
&= ~CW_RC
;
808 control_word
|= RC_CHOP
;
811 /* This should be the most common case */
814 u_char sign
= st0_sign
^ st1_sign
;
815 tag
= FPU_u_div(&st0
, &st1
, &tmp
,
816 PR_64_BITS
| RC_CHOP
| 0x3f,
820 if (exponent(&tmp
) >= 0) {
821 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't
823 q
= significand(&tmp
);
825 rem_kernel(significand(&st0
),
830 setexponent16(&tmp
, exponent16(&st1
));
832 reg_copy(&st0
, &tmp
);
836 if ((round
== RC_RND
)
837 && (tmp
.sigh
& 0xc0000000)) {
838 /* We may need to subtract st(1) once more,
839 to get a result <= 1/2 of st(1). */
840 unsigned long long x
;
842 exponent16(&st1
) - exponent16(&tmp
);
845 x
= significand(&st1
) -
847 else /* expdif is 1 */
848 x
= (significand(&st1
)
851 if ((x
< significand(&tmp
)) ||
852 /* or equi-distant (from 0 & st(1)) and q is odd */
853 ((x
== significand(&tmp
))
855 st0_sign
= !st0_sign
;
856 significand(&tmp
) = x
;
869 control_word
= old_cw
;
874 /* There is a large exponent difference ( >= 64 ) */
875 /* To make much sense, the code in this section should
876 be done at high precision. */
880 /* prevent overflow here */
881 /* N is 'a number between 32 and 63' (p26-113) */
882 reg_copy(&st0
, &tmp
);
884 N
= (expdif
& 0x0000001f) + 32; /* This choice gives results
885 identical to an AMD 486 */
886 setexponent16(&tmp
, N
);
887 exp_1
= exponent16(&st1
);
888 setexponent16(&st1
, 0);
891 sign
= getsign(&tmp
) ^ st1_sign
;
893 FPU_u_div(&tmp
, &st1
, &tmp
,
894 PR_64_BITS
| RC_CHOP
| 0x3f, sign
);
897 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't
900 rem_kernel(significand(&st0
),
903 significand(&tmp
), exponent(&tmp
)
905 setexponent16(&tmp
, exp_1
+ expdif
);
907 /* It is possible for the operation to be complete here.
908 What does the IEEE standard say? The Intel 80486 manual
909 implies that the operation will never be completed at this
910 point, and the behaviour of a real 80486 confirms this.
912 if (!(tmp
.sigh
| tmp
.sigl
)) {
913 /* The result is zero */
914 control_word
= old_cw
;
915 partial_status
= saved_status
;
916 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
917 setsign(&st0
, st0_sign
);
922 #endif /* PECULIAR_486 */
928 control_word
= old_cw
;
929 partial_status
= saved_status
;
930 tag
= FPU_normalize_nuo(&tmp
);
931 reg_copy(&tmp
, st0_ptr
);
933 /* The only condition to be looked for is underflow,
934 and it can occur here only if underflow is unmasked. */
935 if ((exponent16(&tmp
) <= EXP_UNDER
) && (tag
!= TAG_Zero
)
936 && !(control_word
& CW_Underflow
)) {
938 tag
= arith_underflow(st0_ptr
);
939 setsign(st0_ptr
, st0_sign
);
942 } else if ((exponent16(&tmp
) > EXP_UNDER
) || (tag
== TAG_Zero
)) {
944 setsign(st0_ptr
, st0_sign
);
947 FPU_round(st0_ptr
, 0, 0, FULL_PRECISION
, st0_sign
);
955 if (st0_tag
== TAG_Special
)
956 st0_tag
= FPU_Special(st0_ptr
);
957 if (st1_tag
== TAG_Special
)
958 st1_tag
= FPU_Special(st1_ptr
);
960 if (((st0_tag
== TAG_Valid
) && (st1_tag
== TW_Denormal
))
961 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TAG_Valid
))
962 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TW_Denormal
))) {
963 if (denormal_operand() < 0)
966 } else if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
967 FPU_stack_underflow();
969 } else if (st0_tag
== TAG_Zero
) {
970 if (st1_tag
== TAG_Valid
) {
973 } else if (st1_tag
== TW_Denormal
) {
974 if (denormal_operand() < 0)
978 } else if (st1_tag
== TAG_Zero
) {
981 } /* fprem(?,0) always invalid */
982 else if (st1_tag
== TW_Infinity
) {
986 } else if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
987 if (st1_tag
== TAG_Zero
) {
988 arith_invalid(0); /* fprem(Valid,Zero) is invalid */
990 } else if (st1_tag
!= TW_NaN
) {
991 if (((st0_tag
== TW_Denormal
)
992 || (st1_tag
== TW_Denormal
))
993 && (denormal_operand() < 0))
996 if (st1_tag
== TW_Infinity
) {
997 /* fprem(Valid,Infinity) is o.k. */
1002 } else if (st0_tag
== TW_Infinity
) {
1003 if (st1_tag
!= TW_NaN
) {
1004 arith_invalid(0); /* fprem(Infinity,?) is invalid */
1009 /* One of the registers must contain a NaN if we got here. */
1012 if ((st0_tag
!= TW_NaN
) && (st1_tag
!= TW_NaN
))
1013 EXCEPTION(EX_INTERNAL
| 0x118);
1014 #endif /* PARANOID */
1016 real_2op_NaN(st1_ptr
, st1_tag
, 0, st1_ptr
);
1020 /* ST(1) <- ST(1) * log ST; pop ST */
1021 static void fyl2x(FPU_REG
*st0_ptr
, u_char st0_tag
)
1023 FPU_REG
*st1_ptr
= &st(1), exponent
;
1024 u_char st1_tag
= FPU_gettagi(1);
1030 if ((st0_tag
== TAG_Valid
) && (st1_tag
== TAG_Valid
)) {
1032 /* Both regs are Valid or Denormal */
1033 if (signpositive(st0_ptr
)) {
1034 if (st0_tag
== TW_Denormal
)
1035 FPU_to_exp16(st0_ptr
, st0_ptr
);
1037 /* Convert st(0) for internal use. */
1038 setexponent16(st0_ptr
, exponent(st0_ptr
));
1040 if ((st0_ptr
->sigh
== 0x80000000)
1041 && (st0_ptr
->sigl
== 0)) {
1042 /* Special case. The result can be precise. */
1044 e
= exponent16(st0_ptr
);
1053 setexponent16(&exponent
, 31);
1054 tag
= FPU_normalize_nuo(&exponent
);
1056 setsign(&exponent
, esign
);
1058 FPU_mul(&exponent
, tag
, 1, FULL_PRECISION
);
1060 FPU_settagi(1, tag
);
1062 /* The usual case */
1063 sign
= getsign(st1_ptr
);
1064 if (st1_tag
== TW_Denormal
)
1065 FPU_to_exp16(st1_ptr
, st1_ptr
);
1067 /* Convert st(1) for internal use. */
1068 setexponent16(st1_ptr
,
1070 poly_l2(st0_ptr
, st1_ptr
, sign
);
1074 if (arith_invalid(1) < 0)
1083 if (st0_tag
== TAG_Special
)
1084 st0_tag
= FPU_Special(st0_ptr
);
1085 if (st1_tag
== TAG_Special
)
1086 st1_tag
= FPU_Special(st1_ptr
);
1088 if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
1089 FPU_stack_underflow_pop(1);
1091 } else if ((st0_tag
<= TW_Denormal
) && (st1_tag
<= TW_Denormal
)) {
1092 if (st0_tag
== TAG_Zero
) {
1093 if (st1_tag
== TAG_Zero
) {
1094 /* Both args zero is invalid */
1095 if (arith_invalid(1) < 0)
1099 sign
= getsign(st1_ptr
) ^ SIGN_NEG
;
1100 if (FPU_divide_by_zero(1, sign
) < 0)
1103 setsign(st1_ptr
, sign
);
1105 } else if (st1_tag
== TAG_Zero
) {
1106 /* st(1) contains zero, st(0) valid <> 0 */
1107 /* Zero is the valid answer */
1108 sign
= getsign(st1_ptr
);
1110 if (signnegative(st0_ptr
)) {
1112 if (arith_invalid(1) < 0)
1114 } else if ((st0_tag
== TW_Denormal
)
1115 && (denormal_operand() < 0))
1118 if (exponent(st0_ptr
) < 0)
1121 FPU_copy_to_reg1(&CONST_Z
, TAG_Zero
);
1122 setsign(st1_ptr
, sign
);
1125 /* One or both operands are denormals. */
1126 if (denormal_operand() < 0)
1130 } else if ((st0_tag
== TW_NaN
) || (st1_tag
== TW_NaN
)) {
1131 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1134 /* One or both arg must be an infinity */
1135 else if (st0_tag
== TW_Infinity
) {
1136 if ((signnegative(st0_ptr
)) || (st1_tag
== TAG_Zero
)) {
1137 /* log(-infinity) or 0*log(infinity) */
1138 if (arith_invalid(1) < 0)
1141 u_char sign
= getsign(st1_ptr
);
1143 if ((st1_tag
== TW_Denormal
)
1144 && (denormal_operand() < 0))
1147 FPU_copy_to_reg1(&CONST_INF
, TAG_Special
);
1148 setsign(st1_ptr
, sign
);
1151 /* st(1) must be infinity here */
1152 else if (((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
))
1153 && (signpositive(st0_ptr
))) {
1154 if (exponent(st0_ptr
) >= 0) {
1155 if ((exponent(st0_ptr
) == 0) &&
1156 (st0_ptr
->sigh
== 0x80000000) &&
1157 (st0_ptr
->sigl
== 0)) {
1158 /* st(0) holds 1.0 */
1159 /* infinity*log(1) */
1160 if (arith_invalid(1) < 0)
1163 /* else st(0) is positive and > 1.0 */
1165 /* st(0) is positive and < 1.0 */
1167 if ((st0_tag
== TW_Denormal
)
1168 && (denormal_operand() < 0))
1171 changesign(st1_ptr
);
1174 /* st(0) must be zero or negative */
1175 if (st0_tag
== TAG_Zero
) {
1176 /* This should be invalid, but a real 80486 is happy with it. */
1178 #ifndef PECULIAR_486
1179 sign
= getsign(st1_ptr
);
1180 if (FPU_divide_by_zero(1, sign
) < 0)
1182 #endif /* PECULIAR_486 */
1184 changesign(st1_ptr
);
1185 } else if (arith_invalid(1) < 0) /* log(negative) */
1192 static void fpatan(FPU_REG
*st0_ptr
, u_char st0_tag
)
1194 FPU_REG
*st1_ptr
= &st(1);
1195 u_char st1_tag
= FPU_gettagi(1);
1199 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1202 poly_atan(st0_ptr
, st0_tag
, st1_ptr
, st1_tag
);
1209 if (st0_tag
== TAG_Special
)
1210 st0_tag
= FPU_Special(st0_ptr
);
1211 if (st1_tag
== TAG_Special
)
1212 st1_tag
= FPU_Special(st1_ptr
);
1214 if (((st0_tag
== TAG_Valid
) && (st1_tag
== TW_Denormal
))
1215 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TAG_Valid
))
1216 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TW_Denormal
))) {
1217 if (denormal_operand() < 0)
1221 } else if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
1222 FPU_stack_underflow_pop(1);
1224 } else if ((st0_tag
== TW_NaN
) || (st1_tag
== TW_NaN
)) {
1225 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) >= 0)
1228 } else if ((st0_tag
== TW_Infinity
) || (st1_tag
== TW_Infinity
)) {
1229 u_char sign
= getsign(st1_ptr
);
1230 if (st0_tag
== TW_Infinity
) {
1231 if (st1_tag
== TW_Infinity
) {
1232 if (signpositive(st0_ptr
)) {
1233 FPU_copy_to_reg1(&CONST_PI4
, TAG_Valid
);
1235 setpositive(st1_ptr
);
1237 FPU_u_add(&CONST_PI4
, &CONST_PI2
,
1238 st1_ptr
, FULL_PRECISION
,
1240 exponent(&CONST_PI4
),
1241 exponent(&CONST_PI2
));
1243 FPU_settagi(1, tag
);
1246 if ((st1_tag
== TW_Denormal
)
1247 && (denormal_operand() < 0))
1250 if (signpositive(st0_ptr
)) {
1251 FPU_copy_to_reg1(&CONST_Z
, TAG_Zero
);
1252 setsign(st1_ptr
, sign
); /* An 80486 preserves the sign */
1256 FPU_copy_to_reg1(&CONST_PI
, TAG_Valid
);
1260 /* st(1) is infinity, st(0) not infinity */
1261 if ((st0_tag
== TW_Denormal
)
1262 && (denormal_operand() < 0))
1265 FPU_copy_to_reg1(&CONST_PI2
, TAG_Valid
);
1267 setsign(st1_ptr
, sign
);
1268 } else if (st1_tag
== TAG_Zero
) {
1269 /* st(0) must be valid or zero */
1270 u_char sign
= getsign(st1_ptr
);
1272 if ((st0_tag
== TW_Denormal
) && (denormal_operand() < 0))
1275 if (signpositive(st0_ptr
)) {
1276 /* An 80486 preserves the sign */
1281 FPU_copy_to_reg1(&CONST_PI
, TAG_Valid
);
1282 setsign(st1_ptr
, sign
);
1283 } else if (st0_tag
== TAG_Zero
) {
1284 /* st(1) must be TAG_Valid here */
1285 u_char sign
= getsign(st1_ptr
);
1287 if ((st1_tag
== TW_Denormal
) && (denormal_operand() < 0))
1290 FPU_copy_to_reg1(&CONST_PI2
, TAG_Valid
);
1291 setsign(st1_ptr
, sign
);
1295 EXCEPTION(EX_INTERNAL
| 0x125);
1296 #endif /* PARANOID */
1299 set_precision_flag_up(); /* We do not really know if up or down */
1302 static void fprem(FPU_REG
*st0_ptr
, u_char st0_tag
)
1304 do_fprem(st0_ptr
, st0_tag
, RC_CHOP
);
1307 static void fprem1(FPU_REG
*st0_ptr
, u_char st0_tag
)
1309 do_fprem(st0_ptr
, st0_tag
, RC_RND
);
1312 static void fyl2xp1(FPU_REG
*st0_ptr
, u_char st0_tag
)
1315 FPU_REG
*st1_ptr
= &st(1), a
, b
;
1316 u_char st1_tag
= FPU_gettagi(1);
1319 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1322 sign
= getsign(st0_ptr
);
1323 sign1
= getsign(st1_ptr
);
1325 FPU_to_exp16(st0_ptr
, &a
);
1326 FPU_to_exp16(st1_ptr
, &b
);
1328 if (poly_l2p1(sign
, sign1
, &a
, &b
, st1_ptr
))
1335 if (st0_tag
== TAG_Special
)
1336 st0_tag
= FPU_Special(st0_ptr
);
1337 if (st1_tag
== TAG_Special
)
1338 st1_tag
= FPU_Special(st1_ptr
);
1340 if (((st0_tag
== TAG_Valid
) && (st1_tag
== TW_Denormal
))
1341 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TAG_Valid
))
1342 || ((st0_tag
== TW_Denormal
) && (st1_tag
== TW_Denormal
))) {
1343 if (denormal_operand() < 0)
1347 } else if ((st0_tag
== TAG_Empty
) | (st1_tag
== TAG_Empty
)) {
1348 FPU_stack_underflow_pop(1);
1350 } else if (st0_tag
== TAG_Zero
) {
1353 if (denormal_operand() < 0)
1358 setsign(st0_ptr
, getsign(st0_ptr
) ^ getsign(st1_ptr
));
1359 FPU_copy_to_reg1(st0_ptr
, st0_tag
);
1363 /* Infinity*log(1) */
1364 if (arith_invalid(1) < 0)
1369 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1375 EXCEPTION(EX_INTERNAL
| 0x116);
1377 #endif /* PARANOID */
1380 } else if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
1383 if (signnegative(st0_ptr
)) {
1384 if (exponent(st0_ptr
) >= 0) {
1385 /* st(0) holds <= -1.0 */
1386 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1387 changesign(st1_ptr
);
1389 if (arith_invalid(1) < 0)
1391 #endif /* PECULIAR_486 */
1392 } else if ((st0_tag
== TW_Denormal
)
1393 && (denormal_operand() < 0))
1396 changesign(st1_ptr
);
1397 } else if ((st0_tag
== TW_Denormal
)
1398 && (denormal_operand() < 0))
1403 if (signnegative(st0_ptr
)) {
1404 if ((exponent(st0_ptr
) >= 0) &&
1405 !((st0_ptr
->sigh
== 0x80000000) &&
1406 (st0_ptr
->sigl
== 0))) {
1407 /* st(0) holds < -1.0 */
1408 #ifdef PECULIAR_486 /* Stupid 80486 doesn't worry about log(negative). */
1409 changesign(st1_ptr
);
1411 if (arith_invalid(1) < 0)
1413 #endif /* PECULIAR_486 */
1414 } else if ((st0_tag
== TW_Denormal
)
1415 && (denormal_operand() < 0))
1418 changesign(st1_ptr
);
1419 } else if ((st0_tag
== TW_Denormal
)
1420 && (denormal_operand() < 0))
1425 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1429 } else if (st0_tag
== TW_NaN
) {
1430 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1432 } else if (st0_tag
== TW_Infinity
) {
1433 if (st1_tag
== TW_NaN
) {
1434 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1436 } else if (signnegative(st0_ptr
)) {
1437 #ifndef PECULIAR_486
1438 /* This should have higher priority than denormals, but... */
1439 if (arith_invalid(1) < 0) /* log(-infinity) */
1441 #endif /* PECULIAR_486 */
1442 if ((st1_tag
== TW_Denormal
)
1443 && (denormal_operand() < 0))
1446 /* Denormal operands actually get higher priority */
1447 if (arith_invalid(1) < 0) /* log(-infinity) */
1449 #endif /* PECULIAR_486 */
1450 } else if (st1_tag
== TAG_Zero
) {
1452 if (arith_invalid(1) < 0)
1456 /* st(1) must be valid here. */
1458 else if ((st1_tag
== TW_Denormal
) && (denormal_operand() < 0))
1461 /* The Manual says that log(Infinity) is invalid, but a real
1462 80486 sensibly says that it is o.k. */
1464 u_char sign
= getsign(st1_ptr
);
1465 FPU_copy_to_reg1(&CONST_INF
, TAG_Special
);
1466 setsign(st1_ptr
, sign
);
1471 EXCEPTION(EX_INTERNAL
| 0x117);
1474 #endif /* PARANOID */
1481 static void fscale(FPU_REG
*st0_ptr
, u_char st0_tag
)
1483 FPU_REG
*st1_ptr
= &st(1);
1484 u_char st1_tag
= FPU_gettagi(1);
1485 int old_cw
= control_word
;
1486 u_char sign
= getsign(st0_ptr
);
1489 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1493 /* Convert register for internal use. */
1494 setexponent16(st0_ptr
, exponent(st0_ptr
));
1498 if (exponent(st1_ptr
) > 30) {
1499 /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1501 if (signpositive(st1_ptr
)) {
1502 EXCEPTION(EX_Overflow
);
1503 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
1505 EXCEPTION(EX_Underflow
);
1506 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
1508 setsign(st0_ptr
, sign
);
1512 control_word
&= ~CW_RC
;
1513 control_word
|= RC_CHOP
;
1514 reg_copy(st1_ptr
, &tmp
);
1515 FPU_round_to_int(&tmp
, st1_tag
); /* This can never overflow here */
1516 control_word
= old_cw
;
1517 scale
= signnegative(st1_ptr
) ? -tmp
.sigl
: tmp
.sigl
;
1518 scale
+= exponent16(st0_ptr
);
1520 setexponent16(st0_ptr
, scale
);
1522 /* Use FPU_round() to properly detect under/overflow etc */
1523 FPU_round(st0_ptr
, 0, 0, control_word
, sign
);
1528 if (st0_tag
== TAG_Special
)
1529 st0_tag
= FPU_Special(st0_ptr
);
1530 if (st1_tag
== TAG_Special
)
1531 st1_tag
= FPU_Special(st1_ptr
);
1533 if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
1536 /* st(0) must be a denormal */
1537 if ((st0_tag
== TW_Denormal
)
1538 && (denormal_operand() < 0))
1541 FPU_to_exp16(st0_ptr
, st0_ptr
); /* Will not be left on stack */
1545 if (st0_tag
== TW_Denormal
)
1554 if ((st0_tag
== TW_Denormal
)
1555 && (denormal_operand() < 0))
1558 if (signpositive(st1_ptr
))
1559 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
1561 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
1562 setsign(st0_ptr
, sign
);
1566 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1569 } else if (st0_tag
== TAG_Zero
) {
1580 if (signpositive(st1_ptr
))
1581 arith_invalid(0); /* Zero scaled by +Infinity */
1585 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1588 } else if (st0_tag
== TW_Infinity
) {
1599 if (signnegative(st1_ptr
))
1600 arith_invalid(0); /* Infinity scaled by -Infinity */
1604 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1607 } else if (st0_tag
== TW_NaN
) {
1608 if (st1_tag
!= TAG_Empty
) {
1609 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1614 if (!((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
))) {
1615 EXCEPTION(EX_INTERNAL
| 0x115);
1620 /* At least one of st(0), st(1) must be empty */
1621 FPU_stack_underflow();
1625 /*---------------------------------------------------------------------------*/
1627 static FUNC_ST0
const trig_table_a
[] = {
1628 f2xm1
, fyl2x
, fptan
, fpatan
,
1629 fxtract
, fprem1
, (FUNC_ST0
) fdecstp
, (FUNC_ST0
) fincstp
1632 void FPU_triga(void)
1634 (trig_table_a
[FPU_rm
]) (&st(0), FPU_gettag0());
1637 static FUNC_ST0
const trig_table_b
[] = {
1638 fprem
, fyl2xp1
, fsqrt_
, fsincos
, frndint_
, fscale
, (FUNC_ST0
) fsin
, fcos
1641 void FPU_trigb(void)
1643 (trig_table_b
[FPU_rm
]) (&st(0), FPU_gettag0());