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(FPU_REG
*st0_ptr
, u_char st0_tag
)
442 static void fincstp(FPU_REG
*st0_ptr
, u_char st0_tag
)
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 f_sin(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 void fsin(FPU_REG
*st0_ptr
, u_char tag
)
616 static int f_cos(FPU_REG
*st0_ptr
, u_char tag
)
620 st0_sign
= getsign(st0_ptr
);
622 if (tag
== TAG_Valid
) {
625 if (exponent(st0_ptr
) > -40) {
626 if ((exponent(st0_ptr
) < 0)
627 || ((exponent(st0_ptr
) == 0)
628 && (significand(st0_ptr
) <=
629 0xc90fdaa22168c234LL
))) {
632 /* We do not really know if up or down */
633 set_precision_flag_down();
636 } else if ((q
= trig_arg(st0_ptr
, FCOS
)) != -1) {
642 /* We do not really know if up or down */
643 set_precision_flag_down();
647 /* Operand is out of range */
654 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
656 set_precision_flag_down(); /* 80486 appears to do this. */
658 set_precision_flag_up(); /* Must be up. */
659 #endif /* PECULIAR_486 */
662 } else if (tag
== TAG_Zero
) {
663 FPU_copy_to_reg0(&CONST_1
, TAG_Valid
);
668 if (tag
== TAG_Special
)
669 tag
= FPU_Special(st0_ptr
);
671 if (tag
== TW_Denormal
) {
672 if (denormal_operand() < 0)
676 } else if (tag
== TW_Infinity
) {
677 /* The 80486 treats infinity as an invalid operand */
681 single_arg_error(st0_ptr
, tag
); /* requires st0_ptr == &st(0) */
686 static void fcos(FPU_REG
*st0_ptr
, u_char st0_tag
)
688 f_cos(st0_ptr
, st0_tag
);
691 static void fsincos(FPU_REG
*st0_ptr
, u_char st0_tag
)
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);
703 FPU_stack_underflow(); /* Puts a QNaN in the new st(0) */
708 if (STACK_OVERFLOW
) {
709 FPU_stack_overflow();
713 if (st0_tag
== TAG_Special
)
714 tag
= FPU_Special(st0_ptr
);
719 single_arg_2_error(st0_ptr
, TW_NaN
);
721 } else if (tag
== TW_Infinity
) {
722 /* The 80486 treats infinity as an invalid operand */
723 if (arith_invalid(0) >= 0) {
724 /* Masked response */
731 reg_copy(st0_ptr
, &arg
);
732 if (!f_sin(st0_ptr
, st0_tag
)) {
734 FPU_copy_to_reg0(&arg
, st0_tag
);
735 f_cos(&st(0), st0_tag
);
737 /* An error, so restore st(0) */
738 FPU_copy_to_reg0(&arg
, st0_tag
);
742 /*---------------------------------------------------------------------------*/
743 /* The following all require two arguments: st(0) and st(1) */
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
751 static void rem_kernel(unsigned long long st0
, unsigned long long *y
,
752 unsigned long long st1
, unsigned long long q
, int n
)
755 unsigned long long x
;
759 /* Do the required multiplication and subtraction in the one operation */
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]),
765 :"2"(((unsigned *)&st1
)[0]), "m"(((unsigned *)&q
)[0])
767 /* msw x -= msw st1 * lsw q */
768 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x
)[1]),
770 :"1"(((unsigned *)&st1
)[1]), "m"(((unsigned *)&q
)[0])
772 /* msw x -= lsw st1 * msw q */
773 asm volatile ("mull %3; subl %%eax,%0":"=m" (((unsigned *)&x
)[1]),
775 :"1"(((unsigned *)&st1
)[0]), "m"(((unsigned *)&q
)[1])
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. */
784 static void do_fprem(FPU_REG
*st0_ptr
, u_char st0_tag
, int round
)
786 FPU_REG
*st1_ptr
= &st(1);
787 u_char st1_tag
= FPU_gettagi(1);
789 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
790 FPU_REG tmp
, st0
, st1
;
791 u_char st0_sign
, st1_sign
;
797 unsigned short saved_status
;
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
);
806 old_cw
= control_word
;
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
;
816 /* This should be the most common case */
819 u_char sign
= st0_sign
^ st1_sign
;
820 tag
= FPU_u_div(&st0
, &st1
, &tmp
,
821 PR_64_BITS
| RC_CHOP
| 0x3f,
825 if (exponent(&tmp
) >= 0) {
826 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't
828 q
= significand(&tmp
);
830 rem_kernel(significand(&st0
),
835 setexponent16(&tmp
, exponent16(&st1
));
837 reg_copy(&st0
, &tmp
);
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
;
847 exponent16(&st1
) - exponent16(&tmp
);
850 x
= significand(&st1
) -
852 else /* expdif is 1 */
853 x
= (significand(&st1
)
856 if ((x
< significand(&tmp
)) ||
857 /* or equi-distant (from 0 & st(1)) and q is odd */
858 ((x
== significand(&tmp
))
860 st0_sign
= !st0_sign
;
861 significand(&tmp
) = x
;
874 control_word
= old_cw
;
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. */
885 /* prevent overflow here */
886 /* N is 'a number between 32 and 63' (p26-113) */
887 reg_copy(&st0
, &tmp
);
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);
896 sign
= getsign(&tmp
) ^ st1_sign
;
898 FPU_u_div(&tmp
, &st1
, &tmp
,
899 PR_64_BITS
| RC_CHOP
| 0x3f, sign
);
902 FPU_round_to_int(&tmp
, tag
); /* Fortunately, this can't
905 rem_kernel(significand(&st0
),
908 significand(&tmp
), exponent(&tmp
)
910 setexponent16(&tmp
, exp_1
+ expdif
);
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.
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
);
927 #endif /* PECULIAR_486 */
933 control_word
= old_cw
;
934 partial_status
= saved_status
;
935 tag
= FPU_normalize_nuo(&tmp
);
936 reg_copy(&tmp
, st0_ptr
);
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
)) {
943 tag
= arith_underflow(st0_ptr
);
944 setsign(st0_ptr
, st0_sign
);
947 } else if ((exponent16(&tmp
) > EXP_UNDER
) || (tag
== TAG_Zero
)) {
949 setsign(st0_ptr
, st0_sign
);
952 FPU_round(st0_ptr
, 0, 0, FULL_PRECISION
, st0_sign
);
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
);
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)
971 } else if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
972 FPU_stack_underflow();
974 } else if (st0_tag
== TAG_Zero
) {
975 if (st1_tag
== TAG_Valid
) {
978 } else if (st1_tag
== TW_Denormal
) {
979 if (denormal_operand() < 0)
983 } else if (st1_tag
== TAG_Zero
) {
986 } /* fprem(?,0) always invalid */
987 else if (st1_tag
== TW_Infinity
) {
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 */
995 } else if (st1_tag
!= TW_NaN
) {
996 if (((st0_tag
== TW_Denormal
)
997 || (st1_tag
== TW_Denormal
))
998 && (denormal_operand() < 0))
1001 if (st1_tag
== TW_Infinity
) {
1002 /* fprem(Valid,Infinity) is o.k. */
1007 } else if (st0_tag
== TW_Infinity
) {
1008 if (st1_tag
!= TW_NaN
) {
1009 arith_invalid(0); /* fprem(Infinity,?) is invalid */
1014 /* One of the registers must contain a NaN if we got here. */
1017 if ((st0_tag
!= TW_NaN
) && (st1_tag
!= TW_NaN
))
1018 EXCEPTION(EX_INTERNAL
| 0x118);
1019 #endif /* PARANOID */
1021 real_2op_NaN(st1_ptr
, st1_tag
, 0, st1_ptr
);
1025 /* ST(1) <- ST(1) * log ST; pop ST */
1026 static void fyl2x(FPU_REG
*st0_ptr
, u_char st0_tag
)
1028 FPU_REG
*st1_ptr
= &st(1), exponent
;
1029 u_char st1_tag
= FPU_gettagi(1);
1035 if ((st0_tag
== TAG_Valid
) && (st1_tag
== TAG_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
);
1042 /* Convert st(0) for internal use. */
1043 setexponent16(st0_ptr
, exponent(st0_ptr
));
1045 if ((st0_ptr
->sigh
== 0x80000000)
1046 && (st0_ptr
->sigl
== 0)) {
1047 /* Special case. The result can be precise. */
1049 e
= exponent16(st0_ptr
);
1058 setexponent16(&exponent
, 31);
1059 tag
= FPU_normalize_nuo(&exponent
);
1061 setsign(&exponent
, esign
);
1063 FPU_mul(&exponent
, tag
, 1, FULL_PRECISION
);
1065 FPU_settagi(1, tag
);
1067 /* The usual case */
1068 sign
= getsign(st1_ptr
);
1069 if (st1_tag
== TW_Denormal
)
1070 FPU_to_exp16(st1_ptr
, st1_ptr
);
1072 /* Convert st(1) for internal use. */
1073 setexponent16(st1_ptr
,
1075 poly_l2(st0_ptr
, st1_ptr
, sign
);
1079 if (arith_invalid(1) < 0)
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
);
1093 if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
1094 FPU_stack_underflow_pop(1);
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)
1104 sign
= getsign(st1_ptr
) ^ SIGN_NEG
;
1105 if (FPU_divide_by_zero(1, sign
) < 0)
1108 setsign(st1_ptr
, sign
);
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
);
1115 if (signnegative(st0_ptr
)) {
1117 if (arith_invalid(1) < 0)
1119 } else if ((st0_tag
== TW_Denormal
)
1120 && (denormal_operand() < 0))
1123 if (exponent(st0_ptr
) < 0)
1126 FPU_copy_to_reg1(&CONST_Z
, TAG_Zero
);
1127 setsign(st1_ptr
, sign
);
1130 /* One or both operands are denormals. */
1131 if (denormal_operand() < 0)
1135 } else if ((st0_tag
== TW_NaN
) || (st1_tag
== TW_NaN
)) {
1136 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
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)
1146 u_char sign
= getsign(st1_ptr
);
1148 if ((st1_tag
== TW_Denormal
)
1149 && (denormal_operand() < 0))
1152 FPU_copy_to_reg1(&CONST_INF
, TAG_Special
);
1153 setsign(st1_ptr
, sign
);
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)
1168 /* else st(0) is positive and > 1.0 */
1170 /* st(0) is positive and < 1.0 */
1172 if ((st0_tag
== TW_Denormal
)
1173 && (denormal_operand() < 0))
1176 changesign(st1_ptr
);
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. */
1183 #ifndef PECULIAR_486
1184 sign
= getsign(st1_ptr
);
1185 if (FPU_divide_by_zero(1, sign
) < 0)
1187 #endif /* PECULIAR_486 */
1189 changesign(st1_ptr
);
1190 } else if (arith_invalid(1) < 0) /* log(negative) */
1197 static void fpatan(FPU_REG
*st0_ptr
, u_char st0_tag
)
1199 FPU_REG
*st1_ptr
= &st(1);
1200 u_char st1_tag
= FPU_gettagi(1);
1204 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1207 poly_atan(st0_ptr
, st0_tag
, st1_ptr
, st1_tag
);
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
);
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)
1226 } else if ((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
)) {
1227 FPU_stack_underflow_pop(1);
1229 } else if ((st0_tag
== TW_NaN
) || (st1_tag
== TW_NaN
)) {
1230 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) >= 0)
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
);
1240 setpositive(st1_ptr
);
1242 FPU_u_add(&CONST_PI4
, &CONST_PI2
,
1243 st1_ptr
, FULL_PRECISION
,
1245 exponent(&CONST_PI4
),
1246 exponent(&CONST_PI2
));
1248 FPU_settagi(1, tag
);
1251 if ((st1_tag
== TW_Denormal
)
1252 && (denormal_operand() < 0))
1255 if (signpositive(st0_ptr
)) {
1256 FPU_copy_to_reg1(&CONST_Z
, TAG_Zero
);
1257 setsign(st1_ptr
, sign
); /* An 80486 preserves the sign */
1261 FPU_copy_to_reg1(&CONST_PI
, TAG_Valid
);
1265 /* st(1) is infinity, st(0) not infinity */
1266 if ((st0_tag
== TW_Denormal
)
1267 && (denormal_operand() < 0))
1270 FPU_copy_to_reg1(&CONST_PI2
, TAG_Valid
);
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
);
1277 if ((st0_tag
== TW_Denormal
) && (denormal_operand() < 0))
1280 if (signpositive(st0_ptr
)) {
1281 /* An 80486 preserves the sign */
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
);
1292 if ((st1_tag
== TW_Denormal
) && (denormal_operand() < 0))
1295 FPU_copy_to_reg1(&CONST_PI2
, TAG_Valid
);
1296 setsign(st1_ptr
, sign
);
1300 EXCEPTION(EX_INTERNAL
| 0x125);
1301 #endif /* PARANOID */
1304 set_precision_flag_up(); /* We do not really know if up or down */
1307 static void fprem(FPU_REG
*st0_ptr
, u_char st0_tag
)
1309 do_fprem(st0_ptr
, st0_tag
, RC_CHOP
);
1312 static void fprem1(FPU_REG
*st0_ptr
, u_char st0_tag
)
1314 do_fprem(st0_ptr
, st0_tag
, RC_RND
);
1317 static void fyl2xp1(FPU_REG
*st0_ptr
, u_char st0_tag
)
1320 FPU_REG
*st1_ptr
= &st(1), a
, b
;
1321 u_char st1_tag
= FPU_gettagi(1);
1324 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1327 sign
= getsign(st0_ptr
);
1328 sign1
= getsign(st1_ptr
);
1330 FPU_to_exp16(st0_ptr
, &a
);
1331 FPU_to_exp16(st1_ptr
, &b
);
1333 if (poly_l2p1(sign
, sign1
, &a
, &b
, st1_ptr
))
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
);
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)
1352 } else if ((st0_tag
== TAG_Empty
) | (st1_tag
== TAG_Empty
)) {
1353 FPU_stack_underflow_pop(1);
1355 } else if (st0_tag
== TAG_Zero
) {
1358 if (denormal_operand() < 0)
1363 setsign(st0_ptr
, getsign(st0_ptr
) ^ getsign(st1_ptr
));
1364 FPU_copy_to_reg1(st0_ptr
, st0_tag
);
1368 /* Infinity*log(1) */
1369 if (arith_invalid(1) < 0)
1374 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1380 EXCEPTION(EX_INTERNAL
| 0x116);
1382 #endif /* PARANOID */
1385 } else if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
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
);
1394 if (arith_invalid(1) < 0)
1396 #endif /* PECULIAR_486 */
1397 } else if ((st0_tag
== TW_Denormal
)
1398 && (denormal_operand() < 0))
1401 changesign(st1_ptr
);
1402 } else if ((st0_tag
== TW_Denormal
)
1403 && (denormal_operand() < 0))
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
);
1416 if (arith_invalid(1) < 0)
1418 #endif /* PECULIAR_486 */
1419 } else if ((st0_tag
== TW_Denormal
)
1420 && (denormal_operand() < 0))
1423 changesign(st1_ptr
);
1424 } else if ((st0_tag
== TW_Denormal
)
1425 && (denormal_operand() < 0))
1430 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
1434 } else if (st0_tag
== TW_NaN
) {
1435 if (real_2op_NaN(st0_ptr
, st0_tag
, 1, st0_ptr
) < 0)
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)
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) */
1446 #endif /* PECULIAR_486 */
1447 if ((st1_tag
== TW_Denormal
)
1448 && (denormal_operand() < 0))
1451 /* Denormal operands actually get higher priority */
1452 if (arith_invalid(1) < 0) /* log(-infinity) */
1454 #endif /* PECULIAR_486 */
1455 } else if (st1_tag
== TAG_Zero
) {
1457 if (arith_invalid(1) < 0)
1461 /* st(1) must be valid here. */
1463 else if ((st1_tag
== TW_Denormal
) && (denormal_operand() < 0))
1466 /* The Manual says that log(Infinity) is invalid, but a real
1467 80486 sensibly says that it is o.k. */
1469 u_char sign
= getsign(st1_ptr
);
1470 FPU_copy_to_reg1(&CONST_INF
, TAG_Special
);
1471 setsign(st1_ptr
, sign
);
1476 EXCEPTION(EX_INTERNAL
| 0x117);
1479 #endif /* PARANOID */
1486 static void fscale(FPU_REG
*st0_ptr
, u_char st0_tag
)
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
);
1494 if (!((st0_tag
^ TAG_Valid
) | (st1_tag
^ TAG_Valid
))) {
1498 /* Convert register for internal use. */
1499 setexponent16(st0_ptr
, exponent(st0_ptr
));
1503 if (exponent(st1_ptr
) > 30) {
1504 /* 2^31 is far too large, would require 2^(2^30) or 2^(-2^30) */
1506 if (signpositive(st1_ptr
)) {
1507 EXCEPTION(EX_Overflow
);
1508 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
1510 EXCEPTION(EX_Underflow
);
1511 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
1513 setsign(st0_ptr
, sign
);
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
);
1525 setexponent16(st0_ptr
, scale
);
1527 /* Use FPU_round() to properly detect under/overflow etc */
1528 FPU_round(st0_ptr
, 0, 0, control_word
, sign
);
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
);
1538 if ((st0_tag
== TAG_Valid
) || (st0_tag
== TW_Denormal
)) {
1541 /* st(0) must be a denormal */
1542 if ((st0_tag
== TW_Denormal
)
1543 && (denormal_operand() < 0))
1546 FPU_to_exp16(st0_ptr
, st0_ptr
); /* Will not be left on stack */
1550 if (st0_tag
== TW_Denormal
)
1559 if ((st0_tag
== TW_Denormal
)
1560 && (denormal_operand() < 0))
1563 if (signpositive(st1_ptr
))
1564 FPU_copy_to_reg0(&CONST_INF
, TAG_Special
);
1566 FPU_copy_to_reg0(&CONST_Z
, TAG_Zero
);
1567 setsign(st0_ptr
, sign
);
1571 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1574 } else if (st0_tag
== TAG_Zero
) {
1585 if (signpositive(st1_ptr
))
1586 arith_invalid(0); /* Zero scaled by +Infinity */
1590 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1593 } else if (st0_tag
== TW_Infinity
) {
1604 if (signnegative(st1_ptr
))
1605 arith_invalid(0); /* Infinity scaled by -Infinity */
1609 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1612 } else if (st0_tag
== TW_NaN
) {
1613 if (st1_tag
!= TAG_Empty
) {
1614 real_2op_NaN(st1_ptr
, st1_tag
, 0, st0_ptr
);
1619 if (!((st0_tag
== TAG_Empty
) || (st1_tag
== TAG_Empty
))) {
1620 EXCEPTION(EX_INTERNAL
| 0x115);
1625 /* At least one of st(0), st(1) must be empty */
1626 FPU_stack_underflow();
1630 /*---------------------------------------------------------------------------*/
1632 static FUNC_ST0
const trig_table_a
[] = {
1633 f2xm1
, fyl2x
, fptan
, fpatan
,
1634 fxtract
, fprem1
, fdecstp
, fincstp
,
1637 void FPU_triga(void)
1639 (trig_table_a
[FPU_rm
]) (&st(0), FPU_gettag0());
1642 static FUNC_ST0
const trig_table_b
[] = {
1643 fprem
, fyl2xp1
, fsqrt_
, fsincos
, frndint_
, fscale
, fsin
, fcos
1646 void FPU_trigb(void)
1648 (trig_table_b
[FPU_rm
]) (&st(0), FPU_gettag0());