4 * Implementation of the FPU "transcendental" functions.
7 * Copyright (C) 1992,1993,1994
8 * W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
9 * Australia. E-mail billm@vaxc.cc.monash.edu.au
10 * All rights reserved.
12 * This copyright notice covers the redistribution and use of the
13 * FPU emulator developed by W. Metzenthen. It covers only its use
14 * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
15 * use is not permitted under this copyright.
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
20 * 1. Redistributions of source code must retain the above copyright
21 * notice, this list of conditions and the following disclaimer.
22 * 2. Redistributions in binary form must include information specifying
23 * that source code for the emulator is freely available and include
25 * a) an offer to provide the source code for a nominal distribution
27 * b) list at least two alternative methods whereby the source
28 * can be obtained, e.g. a publically accessible bulletin board
29 * and an anonymous ftp site from which the software can be
31 * 3. All advertising materials specifically mentioning features or use of
32 * this emulator must acknowledge that it was developed by W. Metzenthen.
33 * 4. The name of W. Metzenthen may not be used to endorse or promote
34 * products derived from this software without specific prior written
37 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
38 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
39 * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
40 * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
41 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
42 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
43 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
44 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
45 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
46 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
49 * The purpose of this copyright, based upon the Berkeley copyright, is to
50 * ensure that the covered software remains freely available to everyone.
52 * The software (with necessary differences) is also available, but under
53 * the terms of the GNU copyleft, for the Linux operating system and for
54 * the djgpp ms-dos extender.
56 * W. Metzenthen June 1994.
59 * $FreeBSD: src/sys/gnu/i386/fpemul/fpu_trig.c,v 1.10 1999/08/28 00:42:52 peter Exp $
60 * $DragonFly: src/sys/platform/pc32/gnu/fpemul/Attic/fpu_trig.c,v 1.4 2005/01/06 08:33:11 asmodai Exp $
65 #include <sys/param.h>
67 #include <sys/systm.h> /* for printf() in EXCEPTION() */
70 #include <machine/pcb.h>
73 #include "fpu_system.h"
74 #include "exception.h"
76 #include "reg_constant.h"
77 #include "control_w.h"
79 static void convert_l2reg(long *arg, FPU_REG * dest);
87 int old_cw = control_word;
89 control_word &= ~CW_RC;
90 control_word |= RC_CHOP;
93 reg_div(", &CONST_PI2, ", FULL_PRECISION);
95 reg_move(", &tmp);
97 if (tmp.sigh & 0x80000000)
98 return -1; /* |Arg| is >= 2^63 */
99 tmp.exp = EXP_BIAS + 63;
100 q = *(long long *) &(tmp.sigl);
103 reg_sub(", &tmp, X, FULL_PRECISION);
106 control_word = old_cw;
111 /* Convert a long to register */
113 convert_l2reg(long *arg, FPU_REG * dest)
118 reg_move(&CONST_Z, dest);
122 dest->sign = SIGN_POS;
125 dest->sign = SIGN_NEG;
130 dest->exp = EXP_BIAS + 31;
131 dest->tag = TW_Valid;
137 single_arg_error(void)
139 switch (FPU_st0_tag) {
141 if (!(FPU_st0_ptr->sigh & 0x40000000)) { /* Signaling ? */
142 EXCEPTION(EX_Invalid);
143 /* Convert to a QNaN */
144 FPU_st0_ptr->sigh |= 0x40000000;
146 break; /* return with a NaN in st(0) */
148 stack_underflow(); /* Puts a QNaN in st(0) */
152 EXCEPTION(EX_INTERNAL | 0x0112);
153 #endif /* PARANOID */
158 /*---------------------------------------------------------------------------*/
163 switch (FPU_st0_tag) {
168 #ifdef DENORM_OPERAND
169 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
171 #endif /* DENORM_OPERAND */
173 if (FPU_st0_ptr->sign == SIGN_POS) {
174 /* poly_2xm1(x) requires 0 < x < 1. */
175 if (poly_2xm1(FPU_st0_ptr, &rv))
177 reg_mul(&rv, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
179 /* **** Should change poly_2xm1() to at least handle numbers near 0 */
180 /* poly_2xm1(x) doesn't handle negative
182 /* So we compute (poly_2xm1(x+1)-1)/2, for -1
184 reg_add(FPU_st0_ptr, &CONST_1, &tmp, FULL_PRECISION);
185 poly_2xm1(&tmp, &rv);
186 reg_mul(&rv, &tmp, &tmp, FULL_PRECISION);
187 reg_sub(&tmp, &CONST_1, FPU_st0_ptr, FULL_PRECISION);
189 if (FPU_st0_ptr->exp <= EXP_UNDER)
190 arith_underflow(FPU_st0_ptr);
197 if (FPU_st0_ptr->sign == SIGN_NEG) {
198 /* -infinity gives -1 (p16-10) */
199 reg_move(&CONST_1, FPU_st0_ptr);
200 FPU_st0_ptr->sign = SIGN_NEG;
213 char arg_sign = FPU_st0_ptr->sign;
215 if (STACK_OVERFLOW) {
219 switch (FPU_st0_tag) {
222 #ifdef DENORM_OPERAND
223 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
225 #endif /* DENORM_OPERAND */
227 FPU_st0_ptr->sign = SIGN_POS;
228 if ((q = trig_arg(FPU_st0_ptr)) != -1) {
230 reg_sub(&CONST_1, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
232 poly_tan(FPU_st0_ptr, FPU_st0_ptr);
234 FPU_st0_ptr->sign = (q & 1) ^ arg_sign;
236 if (FPU_st0_ptr->exp <= EXP_UNDER)
237 arith_underflow(FPU_st0_ptr);
240 reg_move(&CONST_1, FPU_st0_ptr);
243 /* Operand is out of range */
245 FPU_st0_ptr->sign = arg_sign; /* restore st(0) */
250 /* Operand is out of range */
252 FPU_st0_ptr->sign = arg_sign; /* restore st(0) */
256 reg_move(&CONST_1, FPU_st0_ptr);
270 FPU_REG *st1_ptr = FPU_st0_ptr; /* anticipate */
272 if (STACK_OVERFLOW) {
276 if (!(FPU_st0_tag ^ TW_Valid)) {
279 #ifdef DENORM_OPERAND
280 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
282 #endif /* DENORM_OPERAND */
285 reg_move(st1_ptr, FPU_st0_ptr);
286 FPU_st0_ptr->exp = EXP_BIAS;
287 e = st1_ptr->exp - EXP_BIAS;
288 convert_l2reg(&e, st1_ptr);
291 if (FPU_st0_tag == TW_Zero) {
292 char sign = FPU_st0_ptr->sign;
293 divide_by_zero(SIGN_NEG, FPU_st0_ptr);
295 reg_move(&CONST_Z, FPU_st0_ptr);
296 FPU_st0_ptr->sign = sign;
299 if (FPU_st0_tag == TW_Infinity) {
300 char sign = FPU_st0_ptr->sign;
301 FPU_st0_ptr->sign = SIGN_POS;
303 reg_move(&CONST_INF, FPU_st0_ptr);
304 FPU_st0_ptr->sign = sign;
307 if (FPU_st0_tag == TW_NaN) {
308 if (!(FPU_st0_ptr->sigh & 0x40000000)) { /* Signaling ? */
309 EXCEPTION(EX_Invalid);
310 /* Convert to a QNaN */
311 FPU_st0_ptr->sigh |= 0x40000000;
314 reg_move(st1_ptr, FPU_st0_ptr);
317 if (FPU_st0_tag == TW_Empty) {
318 /* Is this the correct
320 if (control_word & EX_Invalid) {
325 EXCEPTION(EX_StackUnder);
329 EXCEPTION(EX_INTERNAL | 0x119);
330 #endif /* PARANOID */
337 top--; /* FPU_st0_ptr will be fixed in math_emulate()
338 * before the next instr */
344 top++; /* FPU_st0_ptr will be fixed in math_emulate()
345 * before the next instr */
352 if (!(FPU_st0_tag ^ TW_Valid)) {
355 if (FPU_st0_ptr->sign == SIGN_NEG) {
356 arith_invalid(FPU_st0_ptr); /* sqrt(negative) is
360 #ifdef DENORM_OPERAND
361 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
363 #endif /* DENORM_OPERAND */
365 expon = FPU_st0_ptr->exp - EXP_BIAS;
366 FPU_st0_ptr->exp = EXP_BIAS + (expon & 1); /* make st(0) in [1.0
369 wm_sqrt(FPU_st0_ptr, control_word); /* Do the computation */
371 FPU_st0_ptr->exp += expon >> 1;
372 FPU_st0_ptr->sign = SIGN_POS;
374 if (FPU_st0_tag == TW_Zero)
377 if (FPU_st0_tag == TW_Infinity) {
378 if (FPU_st0_ptr->sign == SIGN_NEG)
379 arith_invalid(FPU_st0_ptr); /* sqrt(-Infinity) is
393 if (!(FPU_st0_tag ^ TW_Valid)) {
394 if (FPU_st0_ptr->exp > EXP_BIAS + 63)
397 #ifdef DENORM_OPERAND
398 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
400 #endif /* DENORM_OPERAND */
402 round_to_int(FPU_st0_ptr); /* Fortunately, this can't
403 * overflow to 2^64 */
404 FPU_st0_ptr->exp = EXP_BIAS + 63;
405 normalize(FPU_st0_ptr);
408 if ((FPU_st0_tag == TW_Zero) || (FPU_st0_tag == TW_Infinity))
418 char arg_sign = FPU_st0_ptr->sign;
420 if (FPU_st0_tag == TW_Valid) {
422 FPU_st0_ptr->sign = SIGN_POS;
424 #ifdef DENORM_OPERAND
425 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
427 #endif /* DENORM_OPERAND */
429 if ((q = trig_arg(FPU_st0_ptr)) != -1) {
433 reg_sub(&CONST_1, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
435 poly_sine(FPU_st0_ptr, &rv);
439 rv.sign ^= SIGN_POS ^ SIGN_NEG;
441 reg_move(&rv, FPU_st0_ptr);
443 if (FPU_st0_ptr->exp <= EXP_UNDER)
444 arith_underflow(FPU_st0_ptr);
446 set_precision_flag_up(); /* We do not really know
451 /* Operand is out of range */
453 FPU_st0_ptr->sign = arg_sign; /* restore st(0) */
457 if (FPU_st0_tag == TW_Zero) {
461 if (FPU_st0_tag == TW_Infinity) {
462 /* Operand is out of range */
464 FPU_st0_ptr->sign = arg_sign; /* restore st(0) */
474 char arg_sign = arg->sign;
476 if (arg->tag == TW_Valid) {
479 #ifdef DENORM_OPERAND
480 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
482 #endif /* DENORM_OPERAND */
484 arg->sign = SIGN_POS;
485 if ((q = trig_arg(arg)) != -1) {
489 reg_sub(&CONST_1, arg, arg, FULL_PRECISION);
495 rv.sign ^= SIGN_POS ^ SIGN_NEG;
498 set_precision_flag_up(); /* We do not really know
503 /* Operand is out of range */
505 arg->sign = arg_sign; /* restore st(0) */
509 if (arg->tag == TW_Zero) {
510 reg_move(&CONST_1, arg);
514 if (FPU_st0_tag == TW_Infinity) {
515 /* Operand is out of range */
517 arg->sign = arg_sign; /* restore st(0) */
520 single_arg_error(); /* requires arg ==
540 if (STACK_OVERFLOW) {
544 reg_move(FPU_st0_ptr, &arg);
548 reg_move(&arg, FPU_st0_ptr);
553 /*---------------------------------------------------------------------------*/
554 /* The following all require two arguments: st(0) and st(1) */
556 /* remainder of st(0) / st(1) */
557 /* Assumes that st(0) and st(1) are both TW_Valid */
559 fprem_kernel(int round)
561 FPU_REG *st1_ptr = &st(1);
562 char st1_tag = st1_ptr->tag;
564 if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
566 int old_cw = control_word;
567 int expdif = FPU_st0_ptr->exp - (st1_ptr)->exp;
569 #ifdef DENORM_OPERAND
570 if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
571 (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
573 #endif /* DENORM_OPERAND */
575 control_word &= ~CW_RC;
576 control_word |= round;
579 /* This should be the most common case */
583 reg_div(FPU_st0_ptr, st1_ptr, &tmp, FULL_PRECISION);
585 round_to_int(&tmp); /* Fortunately, this can't
586 * overflow to 2^64 */
587 tmp.exp = EXP_BIAS + 63;
588 q = *(long long *) &(tmp.sigl);
591 reg_mul(st1_ptr, &tmp, &tmp, FULL_PRECISION);
592 reg_sub(FPU_st0_ptr, &tmp, FPU_st0_ptr, FULL_PRECISION);
603 /* There is a large exponent difference ( >= 64 ) */
606 reg_div(FPU_st0_ptr, st1_ptr, &tmp, FULL_PRECISION);
607 /* N is 'a number between 32 and 63' (p26-113) */
608 N_exp = (tmp.exp & 31) + 32;
609 tmp.exp = EXP_BIAS + N_exp;
611 round_to_int(&tmp); /* Fortunately, this can't
612 * overflow to 2^64 */
613 tmp.exp = EXP_BIAS + 63;
616 tmp.exp = EXP_BIAS + expdif - N_exp;
618 reg_mul(st1_ptr, &tmp, &tmp, FULL_PRECISION);
619 reg_sub(FPU_st0_ptr, &tmp, FPU_st0_ptr, FULL_PRECISION);
623 control_word = old_cw;
625 if (FPU_st0_ptr->exp <= EXP_UNDER)
626 arith_underflow(FPU_st0_ptr);
629 if ((FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty)) {
633 if (FPU_st0_tag == TW_Zero) {
634 if (st1_tag == TW_Valid) {
636 #ifdef DENORM_OPERAND
637 if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
639 #endif /* DENORM_OPERAND */
644 if (st1_tag == TW_Zero) {
645 arith_invalid(FPU_st0_ptr);
648 /* fprem(?,0) always invalid */
650 if (st1_tag == TW_Infinity) {
655 if (FPU_st0_tag == TW_Valid) {
656 if (st1_tag == TW_Zero) {
657 arith_invalid(FPU_st0_ptr); /* fprem(Valid,Zero) is
661 if (st1_tag != TW_NaN) {
662 #ifdef DENORM_OPERAND
663 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
665 #endif /* DENORM_OPERAND */
667 if (st1_tag == TW_Infinity) {
676 if (FPU_st0_tag == TW_Infinity) {
677 if (st1_tag != TW_NaN) {
678 arith_invalid(FPU_st0_ptr); /* fprem(Infinity,?) is
683 /* One of the registers must contain a NaN is we got here. */
686 if ((FPU_st0_tag != TW_NaN) && (st1_tag != TW_NaN))
687 EXCEPTION(EX_INTERNAL | 0x118);
688 #endif /* PARANOID */
690 real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
695 /* ST(1) <- ST(1) * log ST; pop ST */
699 FPU_REG *st1_ptr = &st(1);
700 char st1_tag = st1_ptr->tag;
702 if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
703 if (FPU_st0_ptr->sign == SIGN_POS) {
704 int saved_control, saved_status;
706 #ifdef DENORM_OPERAND
707 if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
708 (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
710 #endif /* DENORM_OPERAND */
712 /* We use the general purpose arithmetic, so we need
714 saved_status = status_word;
715 saved_control = control_word;
716 control_word = FULL_PRECISION;
718 poly_l2(FPU_st0_ptr, FPU_st0_ptr);
720 /* Enough of the basic arithmetic is done now */
721 control_word = saved_control;
722 status_word = saved_status;
724 /* Let the multiply set the flags */
725 reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
728 FPU_st0_ptr = &st(0);
732 FPU_st0_ptr = &st(0);
733 arith_invalid(FPU_st0_ptr); /* st(0) cannot be
738 if ((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty)) {
739 stack_underflow_pop(1);
742 if ((FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
743 real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
747 if ((FPU_st0_tag <= TW_Zero) && (st1_tag <= TW_Zero)) {
748 /* one of the args is zero, the other
749 * valid, or both zero */
750 if (FPU_st0_tag == TW_Zero) {
752 FPU_st0_ptr = &st(0);
753 if (FPU_st0_ptr->tag == TW_Zero)
754 arith_invalid(FPU_st0_ptr); /* Both args zero is
758 * specifically covered in the
759 * manual, but divide-by-zero
760 * would seem to be the best
761 * response. However, a real
762 * 80486 does it this way... */
764 if (FPU_st0_ptr->tag == TW_Infinity) {
765 reg_move(&CONST_INF, FPU_st0_ptr);
768 #endif /* PECULIAR_486 */
770 divide_by_zero(st1_ptr->sign ^ SIGN_NEG ^ SIGN_POS, FPU_st0_ptr);
773 /* st(1) contains zero, st(0)
775 /* Zero is the valid answer */
776 char sign = st1_ptr->sign;
778 #ifdef DENORM_OPERAND
779 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
781 #endif /* DENORM_OPERAND */
782 if (FPU_st0_ptr->sign == SIGN_NEG) {
784 FPU_st0_ptr = &st(0);
785 arith_invalid(FPU_st0_ptr); /* log(negative) */
788 if (FPU_st0_ptr->exp < EXP_BIAS)
789 sign ^= SIGN_NEG ^ SIGN_POS;
791 FPU_st0_ptr = &st(0);
792 reg_move(&CONST_Z, FPU_st0_ptr);
793 FPU_st0_ptr->sign = sign;
797 /* One or both arg must be an infinity */
799 if (FPU_st0_tag == TW_Infinity) {
800 if ((FPU_st0_ptr->sign == SIGN_NEG) || (st1_tag == TW_Zero)) {
802 FPU_st0_ptr = &st(0);
803 arith_invalid(FPU_st0_ptr); /* log(-infinity) or
807 char sign = st1_ptr->sign;
809 #ifdef DENORM_OPERAND
810 if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
812 #endif /* DENORM_OPERAND */
815 FPU_st0_ptr = &st(0);
816 reg_move(&CONST_INF, FPU_st0_ptr);
817 FPU_st0_ptr->sign = sign;
821 /* st(1) must be infinity here */
823 if ((FPU_st0_tag == TW_Valid) && (FPU_st0_ptr->sign == SIGN_POS)) {
824 if (FPU_st0_ptr->exp >= EXP_BIAS) {
825 if ((FPU_st0_ptr->exp == EXP_BIAS) &&
826 (FPU_st0_ptr->sigh == 0x80000000) &&
827 (FPU_st0_ptr->sigl == 0)) {
834 FPU_st0_ptr = &st(0);
835 arith_invalid(FPU_st0_ptr); /* infinity*log(1) */
847 #ifdef DENORM_OPERAND
848 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
850 #endif /* DENORM_OPERAND */
852 st1_ptr->sign ^= SIGN_NEG;
857 /* st(0) must be zero
859 if (FPU_st0_ptr->tag == TW_Zero) {
861 FPU_st0_ptr = st1_ptr;
862 st1_ptr->sign ^= SIGN_NEG ^ SIGN_POS;
870 divide_by_zero(st1_ptr->sign, FPU_st0_ptr);
871 #endif /* PECULIAR_486 */
874 FPU_st0_ptr = st1_ptr;
875 arith_invalid(FPU_st0_ptr); /* log(negative) */
885 FPU_REG *st1_ptr = &st(1);
886 char st1_tag = st1_ptr->tag;
888 if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
889 int saved_control, saved_status;
891 int quadrant = st1_ptr->sign | ((FPU_st0_ptr->sign) << 1);
893 #ifdef DENORM_OPERAND
894 if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
895 (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
897 #endif /* DENORM_OPERAND */
899 /* We use the general purpose arithmetic so we need to save
901 saved_status = status_word;
902 saved_control = control_word;
903 control_word = FULL_PRECISION;
905 st1_ptr->sign = FPU_st0_ptr->sign = SIGN_POS;
906 if (compare(st1_ptr) == COMP_A_lt_B) {
908 reg_div(FPU_st0_ptr, st1_ptr, &sum, FULL_PRECISION);
910 reg_div(st1_ptr, FPU_st0_ptr, &sum, FULL_PRECISION);
915 reg_sub(&CONST_PI2, &sum, &sum, FULL_PRECISION);
918 reg_sub(&CONST_PI, &sum, &sum, FULL_PRECISION);
921 sum.sign ^= SIGN_POS ^ SIGN_NEG;
923 /* All of the basic arithmetic is done now */
924 control_word = saved_control;
925 status_word = saved_status;
927 reg_move(&sum, st1_ptr);
929 if ((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty)) {
930 stack_underflow_pop(1);
933 if ((FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
934 real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
938 if ((FPU_st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
939 char sign = st1_ptr->sign;
940 if (FPU_st0_tag == TW_Infinity) {
941 if (st1_tag == TW_Infinity) {
942 if (FPU_st0_ptr->sign == SIGN_POS) {
943 reg_move(&CONST_PI4, st1_ptr);
945 reg_add(&CONST_PI4, &CONST_PI2, st1_ptr, FULL_PRECISION);
948 #ifdef DENORM_OPERAND
949 if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
951 #endif /* DENORM_OPERAND */
953 if (FPU_st0_ptr->sign == SIGN_POS) {
954 reg_move(&CONST_Z, st1_ptr);
958 reg_move(&CONST_PI, st1_ptr);
961 /* st(1) is infinity, st(0)
963 #ifdef DENORM_OPERAND
964 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
966 #endif /* DENORM_OPERAND */
968 reg_move(&CONST_PI2, st1_ptr);
970 st1_ptr->sign = sign;
972 if (st1_tag == TW_Zero) {
973 /* st(0) must be valid or zero */
974 char sign = st1_ptr->sign;
976 #ifdef DENORM_OPERAND
977 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
979 #endif /* DENORM_OPERAND */
981 if (FPU_st0_ptr->sign == SIGN_POS) {
982 reg_move(&CONST_Z, st1_ptr);
986 reg_move(&CONST_PI, st1_ptr);
987 st1_ptr->sign = sign;
989 if (FPU_st0_tag == TW_Zero) {
992 char sign = st1_ptr->sign;
994 #ifdef DENORM_OPERAND
995 if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
997 #endif /* DENORM_OPERAND */
999 reg_move(&CONST_PI2, st1_ptr);
1000 st1_ptr->sign = sign;
1004 EXCEPTION(EX_INTERNAL | 0x220);
1005 #endif /* PARANOID */
1008 set_precision_flag_up();/* We do not really know if up or down */
1015 fprem_kernel(RC_CHOP);
1022 fprem_kernel(RC_RND);
1029 FPU_REG *st1_ptr = &st(1);
1030 char st1_tag = st1_ptr->tag;
1032 if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
1033 int saved_control, saved_status;
1035 #ifdef DENORM_OPERAND
1036 if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
1037 (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
1039 #endif /* DENORM_OPERAND */
1041 /* We use the general purpose arithmetic so we need to save
1043 saved_status = status_word;
1044 saved_control = control_word;
1045 control_word = FULL_PRECISION;
1047 if (poly_l2p1(FPU_st0_ptr, FPU_st0_ptr)) {
1048 arith_invalid(st1_ptr); /* poly_l2p1() returned
1053 /* Enough of the basic arithmetic is done now */
1054 control_word = saved_control;
1055 status_word = saved_status;
1057 /* Let the multiply set the flags */
1058 reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
1062 if ((FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty)) {
1063 stack_underflow_pop(1);
1066 if (FPU_st0_tag == TW_Zero) {
1067 if (st1_tag <= TW_Zero) {
1069 #ifdef DENORM_OPERAND
1070 if ((st1_tag == TW_Valid) && (st1_ptr->exp <= EXP_UNDER) &&
1071 (denormal_operand()))
1073 #endif /* DENORM_OPERAND */
1075 st1_ptr->sign ^= FPU_st0_ptr->sign;
1076 reg_move(FPU_st0_ptr, st1_ptr);
1078 if (st1_tag == TW_Infinity) {
1079 arith_invalid(st1_ptr); /* Infinity*log(1) */
1083 if (st1_tag == TW_NaN) {
1084 real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1090 EXCEPTION(EX_INTERNAL | 0x116);
1093 #endif /* PARANOID */
1097 if (FPU_st0_tag == TW_Valid) {
1098 if (st1_tag == TW_Zero) {
1099 if (FPU_st0_ptr->sign == SIGN_NEG) {
1100 if (FPU_st0_ptr->exp >= EXP_BIAS) {
1103 arith_invalid(st1_ptr); /* infinity*log(1) */
1107 #ifdef DENORM_OPERAND
1108 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1110 #endif /* DENORM_OPERAND */
1111 st1_ptr->sign ^= SIGN_POS ^ SIGN_NEG;
1115 #ifdef DENORM_OPERAND
1116 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1118 #endif /* DENORM_OPERAND */
1122 if (st1_tag == TW_Infinity) {
1123 if (FPU_st0_ptr->sign == SIGN_NEG) {
1124 if ((FPU_st0_ptr->exp >= EXP_BIAS) &&
1125 !((FPU_st0_ptr->sigh == 0x80000000) &&
1126 (FPU_st0_ptr->sigl == 0))) {
1129 arith_invalid(st1_ptr);
1133 #ifdef DENORM_OPERAND
1134 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1136 #endif /* DENORM_OPERAND */
1137 st1_ptr->sign ^= SIGN_POS ^ SIGN_NEG;
1141 #ifdef DENORM_OPERAND
1142 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1144 #endif /* DENORM_OPERAND */
1148 if (st1_tag == TW_NaN) {
1149 real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1154 if (FPU_st0_tag == TW_NaN) {
1155 real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1159 if (FPU_st0_tag == TW_Infinity) {
1160 if (st1_tag == TW_NaN) {
1161 real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
1165 if ((FPU_st0_ptr->sign == SIGN_NEG) ||
1166 (st1_tag == TW_Zero)) {
1167 arith_invalid(st1_ptr); /* log(infinity) */
1171 /* st(1) must be valid
1174 #ifdef DENORM_OPERAND
1175 if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1177 #endif /* DENORM_OPERAND */
1180 * that log(Infinity)
1182 * real 80486 sensibly
1186 char sign = st1_ptr->sign;
1187 reg_move(&CONST_INF, st1_ptr);
1188 st1_ptr->sign = sign;
1195 EXCEPTION(EX_INTERNAL | 0x117);
1197 #endif /* PARANOID */
1204 FPU_REG *st1_ptr = &st(1);
1205 char st1_tag = st1_ptr->tag;
1206 int old_cw = control_word;
1208 if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
1212 #ifdef DENORM_OPERAND
1213 if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
1214 (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
1216 #endif /* DENORM_OPERAND */
1218 if (st1_ptr->exp > EXP_BIAS + 30) {
1219 /* 2^31 is far too large, would require 2^(2^30) or
1223 if (st1_ptr->sign == SIGN_POS) {
1224 EXCEPTION(EX_Overflow);
1225 sign = FPU_st0_ptr->sign;
1226 reg_move(&CONST_INF, FPU_st0_ptr);
1227 FPU_st0_ptr->sign = sign;
1229 EXCEPTION(EX_Underflow);
1230 sign = FPU_st0_ptr->sign;
1231 reg_move(&CONST_Z, FPU_st0_ptr);
1232 FPU_st0_ptr->sign = sign;
1236 control_word &= ~CW_RC;
1237 control_word |= RC_CHOP;
1238 reg_move(st1_ptr, &tmp);
1239 round_to_int(&tmp); /* This can never overflow here */
1240 control_word = old_cw;
1241 scale = st1_ptr->sign ? -tmp.sigl : tmp.sigl;
1242 scale += FPU_st0_ptr->exp;
1243 FPU_st0_ptr->exp = scale;
1245 /* Use round_reg() to properly detect under/overflow etc */
1246 round_reg(FPU_st0_ptr, 0, control_word);
1250 if (FPU_st0_tag == TW_Valid) {
1251 if (st1_tag == TW_Zero) {
1253 #ifdef DENORM_OPERAND
1254 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1256 #endif /* DENORM_OPERAND */
1260 if (st1_tag == TW_Infinity) {
1261 char sign = st1_ptr->sign;
1263 #ifdef DENORM_OPERAND
1264 if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1266 #endif /* DENORM_OPERAND */
1268 if (sign == SIGN_POS) {
1269 reg_move(&CONST_INF, FPU_st0_ptr);
1271 reg_move(&CONST_Z, FPU_st0_ptr);
1272 FPU_st0_ptr->sign = sign;
1275 if (st1_tag == TW_NaN) {
1276 real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
1280 if (FPU_st0_tag == TW_Zero) {
1281 if (st1_tag == TW_Valid) {
1283 #ifdef DENORM_OPERAND
1284 if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1286 #endif /* DENORM_OPERAND */
1290 if (st1_tag == TW_Zero) {
1293 if (st1_tag == TW_Infinity) {
1294 if (st1_ptr->sign == SIGN_NEG)
1297 arith_invalid(FPU_st0_ptr); /* Zero scaled by
1302 if (st1_tag == TW_NaN) {
1303 real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
1307 if (FPU_st0_tag == TW_Infinity) {
1308 if (st1_tag == TW_Valid) {
1310 #ifdef DENORM_OPERAND
1311 if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
1313 #endif /* DENORM_OPERAND */
1317 if (((st1_tag == TW_Infinity) && (st1_ptr->sign == SIGN_POS))
1318 || (st1_tag == TW_Zero))
1321 if (st1_tag == TW_Infinity) {
1322 arith_invalid(FPU_st0_ptr); /* Infinity scaled by
1326 if (st1_tag == TW_NaN) {
1327 real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
1331 if (FPU_st0_tag == TW_NaN) {
1332 if (st1_tag != TW_Empty) {
1333 real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
1338 if (!((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty))) {
1339 EXCEPTION(EX_INTERNAL | 0x115);
1344 /* At least one of st(0), st(1) must be empty */
1350 /*---------------------------------------------------------------------------*/
1352 static FUNC trig_table_a[] = {
1353 f2xm1, fyl2x, fptan, fpatan, fxtract, fprem1, fdecstp, fincstp
1359 (trig_table_a[FPU_rm]) ();
1363 static FUNC trig_table_b[] =
1365 fprem, fyl2xp1, fsqrt_, fsincos, frndint_, emu_fscale, fsin, fcos
1371 (trig_table_b[FPU_rm]) ();