Everyone, say hello to a real floating point emulator!

The fpemul written for Linux by W. Metzenthen: ported to NetBSD and
then to FreeBSD, and now back to OpenBSD.

1 files changed, 1370 insertions, 0 deletions

diff --git a/sys/gnu/arch/i386/fpemul/fpu_trig.c b/sys/gnu/arch/i386/fpemul/fpu_trig.c
new file mode 100644
index 00000000000..8867cb61849
--- /dev/null
+++ b/sys/gnu/arch/i386/fpemul/fpu_trig.c
@@ -0,0 +1,1370 @@
+/*	$OpenBSD: fpu_trig.c,v 1.1 1996/08/27 10:32:50 downsj Exp $	*/
+/*
+ *  fpu_trig.c
+ *
+ * Implementation of the FPU "transcendental" functions.
+ *
+ *
+ * Copyright (C) 1992,1993,1994
+ *                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,
+ *                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au
+ * All rights reserved.
+ *
+ * This copyright notice covers the redistribution and use of the
+ * FPU emulator developed by W. Metzenthen. It covers only its use
+ * in the 386BSD, FreeBSD and NetBSD operating systems. Any other
+ * use is not permitted under this copyright.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must include information specifying
+ *    that source code for the emulator is freely available and include
+ *    either:
+ *      a) an offer to provide the source code for a nominal distribution
+ *         fee, or
+ *      b) list at least two alternative methods whereby the source
+ *         can be obtained, e.g. a publically accessible bulletin board
+ *         and an anonymous ftp site from which the software can be
+ *         downloaded.
+ * 3. All advertising materials specifically mentioning features or use of
+ *    this emulator must acknowledge that it was developed by W. Metzenthen.
+ * 4. The name of W. Metzenthen may not be used to endorse or promote
+ *    products derived from this software without specific prior written
+ *    permission.
+ *
+ * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,
+ * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
+ * AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL
+ * W. METZENTHEN BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ *
+ * The purpose of this copyright, based upon the Berkeley copyright, is to
+ * ensure that the covered software remains freely available to everyone.
+ *
+ * The software (with necessary differences) is also available, but under
+ * the terms of the GNU copyleft, for the Linux operating system and for
+ * the djgpp ms-dos extender.
+ *
+ * W. Metzenthen   June 1994.
+ *
+ *
+ *    $FreeBSD: fpu_trig.c,v 1.5 1995/12/17 21:13:58 phk Exp $
+ *
+ */
+
+
+#include <sys/param.h>
+#include <sys/proc.h>
+#include <machine/cpu.h>
+#include <machine/pcb.h>
+
+#include <gnu/arch/i386/fpemul/fpu_emu.h>
+#include <gnu/arch/i386/fpemul/fpu_system.h>
+#include <gnu/arch/i386/fpemul/exception.h>
+#include <gnu/arch/i386/fpemul/status_w.h>
+#include <gnu/arch/i386/fpemul/reg_constant.h>
+#include <gnu/arch/i386/fpemul/control_w.h>
+
+static void convert_l2reg(long *arg, FPU_REG * dest);
+
+static int
+trig_arg(FPU_REG * X)
+{
+	FPU_REG tmp, quot;
+	int     rv;
+	long long q;
+	int     old_cw = control_word;
+
+	control_word &= ~CW_RC;
+	control_word |= RC_CHOP;
+
+	reg_move(X, &quot);
+	reg_div(&quot, &CONST_PI2, &quot, FULL_PRECISION);
+
+	reg_move(&quot, &tmp);
+	round_to_int(&tmp);
+	if (tmp.sigh & 0x80000000)
+		return -1;	/* |Arg| is >= 2^63 */
+	tmp.exp = EXP_BIAS + 63;
+	q = *(long long *) &(tmp.sigl);
+	normalize(&tmp);
+
+	reg_sub(&quot, &tmp, X, FULL_PRECISION);
+	rv = q & 7;
+
+	control_word = old_cw;
+	return rv;;
+}
+
+
+/* Convert a long to register */
+static void
+convert_l2reg(long *arg, FPU_REG * dest)
+{
+	long    num = *arg;
+
+	if (num == 0) {
+		reg_move(&CONST_Z, dest);
+		return;
+	}
+	if (num > 0)
+		dest->sign = SIGN_POS;
+	else {
+		num = -num;
+		dest->sign = SIGN_NEG;
+	}
+
+	dest->sigh = num;
+	dest->sigl = 0;
+	dest->exp = EXP_BIAS + 31;
+	dest->tag = TW_Valid;
+	normalize(dest);
+}
+
+
+static void
+single_arg_error(void)
+{
+	switch (FPU_st0_tag) {
+		case TW_NaN:
+		if (!(FPU_st0_ptr->sigh & 0x40000000)) {	/* Signaling ? */
+			EXCEPTION(EX_Invalid);
+			/* Convert to a QNaN */
+			FPU_st0_ptr->sigh |= 0x40000000;
+		}
+		break;		/* return with a NaN in st(0) */
+	case TW_Empty:
+		stack_underflow();	/* Puts a QNaN in st(0) */
+		break;
+#ifdef PARANOID
+	default:
+		EXCEPTION(EX_INTERNAL | 0x0112);
+#endif				/* PARANOID */
+	}
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+static void
+f2xm1(void)
+{
+	switch (FPU_st0_tag) {
+		case TW_Valid:
+		{
+			FPU_REG rv, tmp;
+
+#ifdef DENORM_OPERAND
+			if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+				return;
+#endif				/* DENORM_OPERAND */
+
+			if (FPU_st0_ptr->sign == SIGN_POS) {
+				/* poly_2xm1(x) requires 0 < x < 1. */
+				if (poly_2xm1(FPU_st0_ptr, &rv))
+					return;	/* error */
+				reg_mul(&rv, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
+			} else {
+/* **** Should change poly_2xm1() to at least handle numbers near 0 */
+				/* poly_2xm1(x) doesn't handle negative
+				 * numbers. */
+				/* So we compute (poly_2xm1(x+1)-1)/2, for -1
+				 * < x < 0 */
+				reg_add(FPU_st0_ptr, &CONST_1, &tmp, FULL_PRECISION);
+				poly_2xm1(&tmp, &rv);
+				reg_mul(&rv, &tmp, &tmp, FULL_PRECISION);
+				reg_sub(&tmp, &CONST_1, FPU_st0_ptr, FULL_PRECISION);
+				FPU_st0_ptr->exp--;
+				if (FPU_st0_ptr->exp <= EXP_UNDER)
+					arith_underflow(FPU_st0_ptr);
+			}
+			return;
+		}
+	case TW_Zero:
+		return;
+	case TW_Infinity:
+		if (FPU_st0_ptr->sign == SIGN_NEG) {
+			/* -infinity gives -1 (p16-10) */
+			reg_move(&CONST_1, FPU_st0_ptr);
+			FPU_st0_ptr->sign = SIGN_NEG;
+		}
+		return;
+	default:
+		single_arg_error();
+	}
+}
+
+static void
+fptan(void)
+{
+	FPU_REG *st_new_ptr;
+	int     q;
+	char    arg_sign = FPU_st0_ptr->sign;
+
+	if (STACK_OVERFLOW) {
+		stack_overflow();
+		return;
+	}
+	switch (FPU_st0_tag) {
+	case TW_Valid:
+
+#ifdef DENORM_OPERAND
+		if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		FPU_st0_ptr->sign = SIGN_POS;
+		if ((q = trig_arg(FPU_st0_ptr)) != -1) {
+			if (q & 1)
+				reg_sub(&CONST_1, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
+
+			poly_tan(FPU_st0_ptr, FPU_st0_ptr);
+
+			FPU_st0_ptr->sign = (q & 1) ^ arg_sign;
+
+			if (FPU_st0_ptr->exp <= EXP_UNDER)
+				arith_underflow(FPU_st0_ptr);
+
+			push();
+			reg_move(&CONST_1, FPU_st0_ptr);
+			setcc(0);
+		} else {
+			/* Operand is out of range */
+			setcc(SW_C2);
+			FPU_st0_ptr->sign = arg_sign;	/* restore st(0) */
+			return;
+		}
+		break;
+	case TW_Infinity:
+		/* Operand is out of range */
+		setcc(SW_C2);
+		FPU_st0_ptr->sign = arg_sign;	/* restore st(0) */
+		return;
+	case TW_Zero:
+		push();
+		reg_move(&CONST_1, FPU_st0_ptr);
+		setcc(0);
+		break;
+	default:
+		single_arg_error();
+		break;
+	}
+}
+
+
+static void
+fxtract(void)
+{
+	FPU_REG *st_new_ptr;
+	register FPU_REG *st1_ptr = FPU_st0_ptr;	/* anticipate */
+
+	if (STACK_OVERFLOW) {
+		stack_overflow();
+		return;
+	}
+	if (!(FPU_st0_tag ^ TW_Valid)) {
+		long    e;
+
+#ifdef DENORM_OPERAND
+		if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		push();
+		reg_move(st1_ptr, FPU_st0_ptr);
+		FPU_st0_ptr->exp = EXP_BIAS;
+		e = st1_ptr->exp - EXP_BIAS;
+		convert_l2reg(&e, st1_ptr);
+		return;
+	} else
+		if (FPU_st0_tag == TW_Zero) {
+			char    sign = FPU_st0_ptr->sign;
+			divide_by_zero(SIGN_NEG, FPU_st0_ptr);
+			push();
+			reg_move(&CONST_Z, FPU_st0_ptr);
+			FPU_st0_ptr->sign = sign;
+			return;
+		} else
+			if (FPU_st0_tag == TW_Infinity) {
+				char    sign = FPU_st0_ptr->sign;
+				FPU_st0_ptr->sign = SIGN_POS;
+				push();
+				reg_move(&CONST_INF, FPU_st0_ptr);
+				FPU_st0_ptr->sign = sign;
+				return;
+			} else
+				if (FPU_st0_tag == TW_NaN) {
+					if (!(FPU_st0_ptr->sigh & 0x40000000)) {	/* Signaling ? */
+						EXCEPTION(EX_Invalid);
+						/* Convert to a QNaN */
+						FPU_st0_ptr->sigh |= 0x40000000;
+					}
+					push();
+					reg_move(st1_ptr, FPU_st0_ptr);
+					return;
+				} else
+					if (FPU_st0_tag == TW_Empty) {
+						/* Is this the correct
+						 * behaviour? */
+						if (control_word & EX_Invalid) {
+							stack_underflow();
+							push();
+							stack_underflow();
+						} else
+							EXCEPTION(EX_StackUnder);
+					}
+#ifdef PARANOID
+					else
+						EXCEPTION(EX_INTERNAL | 0x119);
+#endif				/* PARANOID */
+}
+
+
+static void
+fdecstp(void)
+{
+	top--;			/* FPU_st0_ptr will be fixed in math_emulate()
+				 * before the next instr */
+}
+
+static void
+fincstp(void)
+{
+	top++;			/* FPU_st0_ptr will be fixed in math_emulate()
+				 * before the next instr */
+}
+
+
+static void
+fsqrt_(void)
+{
+	if (!(FPU_st0_tag ^ TW_Valid)) {
+		int     expon;
+
+		if (FPU_st0_ptr->sign == SIGN_NEG) {
+			arith_invalid(FPU_st0_ptr);	/* sqrt(negative) is
+							 * invalid */
+			return;
+		}
+#ifdef DENORM_OPERAND
+		if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		expon = FPU_st0_ptr->exp - EXP_BIAS;
+		FPU_st0_ptr->exp = EXP_BIAS + (expon & 1);	/* make st(0) in  [1.0
+								 * .. 4.0) */
+
+		wm_sqrt(FPU_st0_ptr, control_word);	/* Do the computation */
+
+		FPU_st0_ptr->exp += expon >> 1;
+		FPU_st0_ptr->sign = SIGN_POS;
+	} else
+		if (FPU_st0_tag == TW_Zero)
+			return;
+		else
+			if (FPU_st0_tag == TW_Infinity) {
+				if (FPU_st0_ptr->sign == SIGN_NEG)
+					arith_invalid(FPU_st0_ptr);	/* sqrt(-Infinity) is
+									 * invalid */
+				return;
+			} else {
+				single_arg_error();
+				return;
+			}
+
+}
+
+
+static void
+frndint_(void)
+{
+	if (!(FPU_st0_tag ^ TW_Valid)) {
+		if (FPU_st0_ptr->exp > EXP_BIAS + 63)
+			return;
+
+#ifdef DENORM_OPERAND
+		if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		round_to_int(FPU_st0_ptr);	/* Fortunately, this can't
+						 * overflow to 2^64 */
+		FPU_st0_ptr->exp = EXP_BIAS + 63;
+		normalize(FPU_st0_ptr);
+		return;
+	} else
+		if ((FPU_st0_tag == TW_Zero) || (FPU_st0_tag == TW_Infinity))
+			return;
+		else
+			single_arg_error();
+}
+
+
+static void
+fsin(void)
+{
+	char    arg_sign = FPU_st0_ptr->sign;
+
+	if (FPU_st0_tag == TW_Valid) {
+		int     q;
+		FPU_st0_ptr->sign = SIGN_POS;
+
+#ifdef DENORM_OPERAND
+		if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		if ((q = trig_arg(FPU_st0_ptr)) != -1) {
+			FPU_REG rv;
+
+			if (q & 1)
+				reg_sub(&CONST_1, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
+
+			poly_sine(FPU_st0_ptr, &rv);
+
+			setcc(0);
+			if (q & 2)
+				rv.sign ^= SIGN_POS ^ SIGN_NEG;
+			rv.sign ^= arg_sign;
+			reg_move(&rv, FPU_st0_ptr);
+
+			if (FPU_st0_ptr->exp <= EXP_UNDER)
+				arith_underflow(FPU_st0_ptr);
+
+			set_precision_flag_up();	/* We do not really know
+							 * if up or down */
+
+			return;
+		} else {
+			/* Operand is out of range */
+			setcc(SW_C2);
+			FPU_st0_ptr->sign = arg_sign;	/* restore st(0) */
+			return;
+		}
+	} else
+		if (FPU_st0_tag == TW_Zero) {
+			setcc(0);
+			return;
+		} else
+			if (FPU_st0_tag == TW_Infinity) {
+				/* Operand is out of range */
+				setcc(SW_C2);
+				FPU_st0_ptr->sign = arg_sign;	/* restore st(0) */
+				return;
+			} else
+				single_arg_error();
+}
+
+
+static int
+f_cos(FPU_REG * arg)
+{
+	char    arg_sign = arg->sign;
+
+	if (arg->tag == TW_Valid) {
+		int     q;
+
+#ifdef DENORM_OPERAND
+		if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+			return 1;
+#endif				/* DENORM_OPERAND */
+
+		arg->sign = SIGN_POS;
+		if ((q = trig_arg(arg)) != -1) {
+			FPU_REG rv;
+
+			if (!(q & 1))
+				reg_sub(&CONST_1, arg, arg, FULL_PRECISION);
+
+			poly_sine(arg, &rv);
+
+			setcc(0);
+			if ((q + 1) & 2)
+				rv.sign ^= SIGN_POS ^ SIGN_NEG;
+			reg_move(&rv, arg);
+
+			set_precision_flag_up();	/* We do not really know
+							 * if up or down */
+
+			return 0;
+		} else {
+			/* Operand is out of range */
+			setcc(SW_C2);
+			arg->sign = arg_sign;	/* restore st(0) */
+			return 1;
+		}
+	} else
+		if (arg->tag == TW_Zero) {
+			reg_move(&CONST_1, arg);
+			setcc(0);
+			return 0;
+		} else
+			if (FPU_st0_tag == TW_Infinity) {
+				/* Operand is out of range */
+				setcc(SW_C2);
+				arg->sign = arg_sign;	/* restore st(0) */
+				return 1;
+			} else {
+				single_arg_error();	/* requires arg ==
+							 * &st(0) */
+				return 1;
+			}
+}
+
+
+static void
+fcos(void)
+{
+	f_cos(FPU_st0_ptr);
+}
+
+
+static void
+fsincos(void)
+{
+	FPU_REG *st_new_ptr;
+	FPU_REG arg;
+
+	if (STACK_OVERFLOW) {
+		stack_overflow();
+		return;
+	}
+	reg_move(FPU_st0_ptr, &arg);
+	if (!f_cos(&arg)) {
+		fsin();
+		push();
+		reg_move(&arg, FPU_st0_ptr);
+	}
+}
+
+
+/*---------------------------------------------------------------------------*/
+/* The following all require two arguments: st(0) and st(1) */
+
+/* remainder of st(0) / st(1) */
+/* Assumes that st(0) and st(1) are both TW_Valid */
+static void
+fprem_kernel(int round)
+{
+	FPU_REG *st1_ptr = &st(1);
+	char    st1_tag = st1_ptr->tag;
+
+	if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
+		FPU_REG tmp;
+		int     old_cw = control_word;
+		int     expdif = FPU_st0_ptr->exp - (st1_ptr)->exp;
+
+#ifdef DENORM_OPERAND
+		if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
+			(st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		control_word &= ~CW_RC;
+		control_word |= round;
+
+		if (expdif < 64) {
+			/* This should be the most common case */
+			long long q;
+			int     c = 0;
+
+			reg_div(FPU_st0_ptr, st1_ptr, &tmp, FULL_PRECISION);
+
+			round_to_int(&tmp);	/* Fortunately, this can't
+						 * overflow to 2^64 */
+			tmp.exp = EXP_BIAS + 63;
+			q = *(long long *) &(tmp.sigl);
+			normalize(&tmp);
+
+			reg_mul(st1_ptr, &tmp, &tmp, FULL_PRECISION);
+			reg_sub(FPU_st0_ptr, &tmp, FPU_st0_ptr, FULL_PRECISION);
+
+			if (q & 4)
+				c |= SW_C3;
+			if (q & 2)
+				c |= SW_C1;
+			if (q & 1)
+				c |= SW_C0;
+
+			setcc(c);
+		} else {
+			/* There is a large exponent difference ( >= 64 ) */
+			int     N_exp;
+
+			reg_div(FPU_st0_ptr, st1_ptr, &tmp, FULL_PRECISION);
+			/* N is 'a number between 32 and 63' (p26-113) */
+			N_exp = (tmp.exp & 31) + 32;
+			tmp.exp = EXP_BIAS + N_exp;
+
+			round_to_int(&tmp);	/* Fortunately, this can't
+						 * overflow to 2^64 */
+			tmp.exp = EXP_BIAS + 63;
+			normalize(&tmp);
+
+			tmp.exp = EXP_BIAS + expdif - N_exp;
+
+			reg_mul(st1_ptr, &tmp, &tmp, FULL_PRECISION);
+			reg_sub(FPU_st0_ptr, &tmp, FPU_st0_ptr, FULL_PRECISION);
+
+			setcc(SW_C2);
+		}
+		control_word = old_cw;
+
+		if (FPU_st0_ptr->exp <= EXP_UNDER)
+			arith_underflow(FPU_st0_ptr);
+		return;
+	} else
+		if ((FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty)) {
+			stack_underflow();
+			return;
+		} else
+			if (FPU_st0_tag == TW_Zero) {
+				if (st1_tag == TW_Valid) {
+
+#ifdef DENORM_OPERAND
+					if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+						return;
+#endif				/* DENORM_OPERAND */
+
+					setcc(0);
+					return;
+				} else
+					if (st1_tag == TW_Zero) {
+						arith_invalid(FPU_st0_ptr);
+						return;
+					}
+				/* fprem(?,0) always invalid */
+					else
+						if (st1_tag == TW_Infinity) {
+							setcc(0);
+							return;
+						}
+			} else
+				if (FPU_st0_tag == TW_Valid) {
+					if (st1_tag == TW_Zero) {
+						arith_invalid(FPU_st0_ptr);	/* fprem(Valid,Zero) is
+										 * invalid */
+						return;
+					} else
+						if (st1_tag != TW_NaN) {
+#ifdef DENORM_OPERAND
+							if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+								return;
+#endif				/* DENORM_OPERAND */
+
+							if (st1_tag == TW_Infinity) {
+								/* fprem(Valid,
+								 * Infinity)
+								 * is o.k. */
+								setcc(0);
+								return;
+							}
+						}
+				} else
+					if (FPU_st0_tag == TW_Infinity) {
+						if (st1_tag != TW_NaN) {
+							arith_invalid(FPU_st0_ptr);	/* fprem(Infinity,?) is
+											 * invalid */
+							return;
+						}
+					}
+	/* One of the registers must contain a NaN is we got here. */
+
+#ifdef PARANOID
+	if ((FPU_st0_tag != TW_NaN) && (st1_tag != TW_NaN))
+		EXCEPTION(EX_INTERNAL | 0x118);
+#endif				/* PARANOID */
+
+	real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
+
+}
+
+
+/* ST(1) <- ST(1) * log ST;  pop ST */
+static void
+fyl2x(void)
+{
+	FPU_REG *st1_ptr = &st(1);
+	char    st1_tag = st1_ptr->tag;
+
+	if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
+		if (FPU_st0_ptr->sign == SIGN_POS) {
+			int     saved_control, saved_status;
+
+#ifdef DENORM_OPERAND
+			if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
+				(st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
+				return;
+#endif				/* DENORM_OPERAND */
+
+			/* We use the general purpose arithmetic, so we need
+			 * to save these. */
+			saved_status = status_word;
+			saved_control = control_word;
+			control_word = FULL_PRECISION;
+
+			poly_l2(FPU_st0_ptr, FPU_st0_ptr);
+
+			/* Enough of the basic arithmetic is done now */
+			control_word = saved_control;
+			status_word = saved_status;
+
+			/* Let the multiply set the flags */
+			reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
+
+			pop();
+			FPU_st0_ptr = &st(0);
+		} else {
+			/* negative	 */
+			pop();
+			FPU_st0_ptr = &st(0);
+			arith_invalid(FPU_st0_ptr);	/* st(0) cannot be
+							 * negative */
+			return;
+		}
+	} else
+		if ((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty)) {
+			stack_underflow_pop(1);
+			return;
+		} else
+			if ((FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
+				real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
+				pop();
+				return;
+			} else
+				if ((FPU_st0_tag <= TW_Zero) && (st1_tag <= TW_Zero)) {
+					/* one of the args is zero, the other
+					 * valid, or both zero */
+					if (FPU_st0_tag == TW_Zero) {
+						pop();
+						FPU_st0_ptr = &st(0);
+						if (FPU_st0_ptr->tag == TW_Zero)
+							arith_invalid(FPU_st0_ptr);	/* Both args zero is
+											 * invalid */
+#ifdef PECULIAR_486
+						/* This case is not
+						 * specifically covered in the
+						 * manual, but divide-by-zero
+						 * would seem to be the best
+						 * response. However, a real
+						 * 80486 does it this way... */
+						else
+							if (FPU_st0_ptr->tag == TW_Infinity) {
+								reg_move(&CONST_INF, FPU_st0_ptr);
+								return;
+							}
+#endif				/* PECULIAR_486 */
+							else
+								divide_by_zero(st1_ptr->sign ^ SIGN_NEG ^ SIGN_POS, FPU_st0_ptr);
+						return;
+					} else {
+						/* st(1) contains zero, st(0)
+						 * valid <> 0 */
+						/* Zero is the valid answer */
+						char    sign = st1_ptr->sign;
+
+#ifdef DENORM_OPERAND
+						if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+							return;
+#endif				/* DENORM_OPERAND */
+						if (FPU_st0_ptr->sign == SIGN_NEG) {
+							pop();
+							FPU_st0_ptr = &st(0);
+							arith_invalid(FPU_st0_ptr);	/* log(negative) */
+							return;
+						}
+						if (FPU_st0_ptr->exp < EXP_BIAS)
+							sign ^= SIGN_NEG ^ SIGN_POS;
+						pop();
+						FPU_st0_ptr = &st(0);
+						reg_move(&CONST_Z, FPU_st0_ptr);
+						FPU_st0_ptr->sign = sign;
+						return;
+					}
+				}
+	/* One or both arg must be an infinity */
+				else
+					if (FPU_st0_tag == TW_Infinity) {
+						if ((FPU_st0_ptr->sign == SIGN_NEG) || (st1_tag == TW_Zero)) {
+							pop();
+							FPU_st0_ptr = &st(0);
+							arith_invalid(FPU_st0_ptr);	/* log(-infinity) or
+											 * 0*log(infinity) */
+							return;
+						} else {
+							char    sign = st1_ptr->sign;
+
+#ifdef DENORM_OPERAND
+							if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+								return;
+#endif				/* DENORM_OPERAND */
+
+							pop();
+							FPU_st0_ptr = &st(0);
+							reg_move(&CONST_INF, FPU_st0_ptr);
+							FPU_st0_ptr->sign = sign;
+							return;
+						}
+					}
+	/* st(1) must be infinity here */
+					else
+						if ((FPU_st0_tag == TW_Valid) && (FPU_st0_ptr->sign == SIGN_POS)) {
+							if (FPU_st0_ptr->exp >= EXP_BIAS) {
+								if ((FPU_st0_ptr->exp == EXP_BIAS) &&
+								    (FPU_st0_ptr->sigh == 0x80000000) &&
+								    (FPU_st0_ptr->sigl == 0)) {
+									/* st(0
+									 * )
+									 * hold
+									 * s
+									 * 1.0 */
+									pop();
+									FPU_st0_ptr = &st(0);
+									arith_invalid(FPU_st0_ptr);	/* infinity*log(1) */
+									return;
+								}
+								/* st(0) is
+								 * positive
+								 * and > 1.0 */
+								pop();
+							} else {
+								/* st(0) is
+								 * positive
+								 * and < 1.0 */
+
+#ifdef DENORM_OPERAND
+								if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+									return;
+#endif				/* DENORM_OPERAND */
+
+								st1_ptr->sign ^= SIGN_NEG;
+								pop();
+							}
+							return;
+						} else {
+							/* st(0) must be zero
+							 * or negative */
+							if (FPU_st0_ptr->tag == TW_Zero) {
+								pop();
+								FPU_st0_ptr = st1_ptr;
+								st1_ptr->sign ^= SIGN_NEG ^ SIGN_POS;
+								/* This should
+								 * be invalid,
+								 * but a real
+								 * 80486 is
+								 * happy with
+								 * it. */
+#ifndef PECULIAR_486
+								divide_by_zero(st1_ptr->sign, FPU_st0_ptr);
+#endif				/* PECULIAR_486 */
+							} else {
+								pop();
+								FPU_st0_ptr = st1_ptr;
+								arith_invalid(FPU_st0_ptr);	/* log(negative) */
+							}
+							return;
+						}
+}
+
+
+static void
+fpatan(void)
+{
+	FPU_REG *st1_ptr = &st(1);
+	char    st1_tag = st1_ptr->tag;
+
+	if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
+		int     saved_control, saved_status;
+		FPU_REG sum;
+		int     quadrant = st1_ptr->sign | ((FPU_st0_ptr->sign) << 1);
+
+#ifdef DENORM_OPERAND
+		if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
+			(st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		/* We use the general purpose arithmetic so we need to save
+		 * these. */
+		saved_status = status_word;
+		saved_control = control_word;
+		control_word = FULL_PRECISION;
+
+		st1_ptr->sign = FPU_st0_ptr->sign = SIGN_POS;
+		if (compare(st1_ptr) == COMP_A_lt_B) {
+			quadrant |= 4;
+			reg_div(FPU_st0_ptr, st1_ptr, &sum, FULL_PRECISION);
+		} else
+			reg_div(st1_ptr, FPU_st0_ptr, &sum, FULL_PRECISION);
+
+		poly_atan(&sum);
+
+		if (quadrant & 4) {
+			reg_sub(&CONST_PI2, &sum, &sum, FULL_PRECISION);
+		}
+		if (quadrant & 2) {
+			reg_sub(&CONST_PI, &sum, &sum, FULL_PRECISION);
+		}
+		if (quadrant & 1)
+			sum.sign ^= SIGN_POS ^ SIGN_NEG;
+
+		/* All of the basic arithmetic is done now */
+		control_word = saved_control;
+		status_word = saved_status;
+
+		reg_move(&sum, st1_ptr);
+	} else
+		if ((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty)) {
+			stack_underflow_pop(1);
+			return;
+		} else
+			if ((FPU_st0_tag == TW_NaN) || (st1_tag == TW_NaN)) {
+				real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
+				pop();
+				return;
+			} else
+				if ((FPU_st0_tag == TW_Infinity) || (st1_tag == TW_Infinity)) {
+					char    sign = st1_ptr->sign;
+					if (FPU_st0_tag == TW_Infinity) {
+						if (st1_tag == TW_Infinity) {
+							if (FPU_st0_ptr->sign == SIGN_POS) {
+								reg_move(&CONST_PI4, st1_ptr);
+							} else
+								reg_add(&CONST_PI4, &CONST_PI2, st1_ptr, FULL_PRECISION);
+						} else {
+
+#ifdef DENORM_OPERAND
+							if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+								return;
+#endif				/* DENORM_OPERAND */
+
+							if (FPU_st0_ptr->sign == SIGN_POS) {
+								reg_move(&CONST_Z, st1_ptr);
+								pop();
+								return;
+							} else
+								reg_move(&CONST_PI, st1_ptr);
+						}
+					} else {
+						/* st(1) is infinity, st(0)
+						 * not infinity */
+#ifdef DENORM_OPERAND
+						if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+							return;
+#endif				/* DENORM_OPERAND */
+
+						reg_move(&CONST_PI2, st1_ptr);
+					}
+					st1_ptr->sign = sign;
+				} else
+					if (st1_tag == TW_Zero) {
+						/* st(0) must be valid or zero */
+						char    sign = st1_ptr->sign;
+
+#ifdef DENORM_OPERAND
+						if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+							return;
+#endif				/* DENORM_OPERAND */
+
+						if (FPU_st0_ptr->sign == SIGN_POS) {
+							reg_move(&CONST_Z, st1_ptr);
+							pop();
+							return;
+						} else
+							reg_move(&CONST_PI, st1_ptr);
+						st1_ptr->sign = sign;
+					} else
+						if (FPU_st0_tag == TW_Zero) {
+							/* st(1) must be
+							 * TW_Valid here */
+							char    sign = st1_ptr->sign;
+
+#ifdef DENORM_OPERAND
+							if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+								return;
+#endif				/* DENORM_OPERAND */
+
+							reg_move(&CONST_PI2, st1_ptr);
+							st1_ptr->sign = sign;
+						}
+#ifdef PARANOID
+						else
+							EXCEPTION(EX_INTERNAL | 0x220);
+#endif				/* PARANOID */
+
+	pop();
+	set_precision_flag_up();/* We do not really know if up or down */
+}
+
+
+static void
+fprem(void)
+{
+	fprem_kernel(RC_CHOP);
+}
+
+
+static void
+fprem1(void)
+{
+	fprem_kernel(RC_RND);
+}
+
+
+static void
+fyl2xp1(void)
+{
+	FPU_REG *st1_ptr = &st(1);
+	char    st1_tag = st1_ptr->tag;
+
+	if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
+		int     saved_control, saved_status;
+
+#ifdef DENORM_OPERAND
+		if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
+			(st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		/* We use the general purpose arithmetic so we need to save
+		 * these. */
+		saved_status = status_word;
+		saved_control = control_word;
+		control_word = FULL_PRECISION;
+
+		if (poly_l2p1(FPU_st0_ptr, FPU_st0_ptr)) {
+			arith_invalid(st1_ptr);	/* poly_l2p1() returned
+						 * invalid */
+			pop();
+			return;
+		}
+		/* Enough of the basic arithmetic is done now */
+		control_word = saved_control;
+		status_word = saved_status;
+
+		/* Let the multiply set the flags */
+		reg_mul(FPU_st0_ptr, st1_ptr, st1_ptr, FULL_PRECISION);
+
+		pop();
+	} else
+		if ((FPU_st0_tag == TW_Empty) | (st1_tag == TW_Empty)) {
+			stack_underflow_pop(1);
+			return;
+		} else
+			if (FPU_st0_tag == TW_Zero) {
+				if (st1_tag <= TW_Zero) {
+
+#ifdef DENORM_OPERAND
+					if ((st1_tag == TW_Valid) && (st1_ptr->exp <= EXP_UNDER) &&
+					    (denormal_operand()))
+						return;
+#endif				/* DENORM_OPERAND */
+
+					st1_ptr->sign ^= FPU_st0_ptr->sign;
+					reg_move(FPU_st0_ptr, st1_ptr);
+				} else
+					if (st1_tag == TW_Infinity) {
+						arith_invalid(st1_ptr);	/* Infinity*log(1) */
+						pop();
+						return;
+					} else
+						if (st1_tag == TW_NaN) {
+							real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
+							pop();
+							return;
+						}
+#ifdef PARANOID
+						else {
+							EXCEPTION(EX_INTERNAL | 0x116);
+							return;
+						}
+#endif				/* PARANOID */
+				pop();
+				return;
+			} else
+				if (FPU_st0_tag == TW_Valid) {
+					if (st1_tag == TW_Zero) {
+						if (FPU_st0_ptr->sign == SIGN_NEG) {
+							if (FPU_st0_ptr->exp >= EXP_BIAS) {
+								/* st(0) holds
+								 * <= -1.0 */
+								arith_invalid(st1_ptr);	/* infinity*log(1) */
+								pop();
+								return;
+							}
+#ifdef DENORM_OPERAND
+							if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+								return;
+#endif				/* DENORM_OPERAND */
+							st1_ptr->sign ^= SIGN_POS ^ SIGN_NEG;
+							pop();
+							return;
+						}
+#ifdef DENORM_OPERAND
+						if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+							return;
+#endif				/* DENORM_OPERAND */
+						pop();
+						return;
+					}
+					if (st1_tag == TW_Infinity) {
+						if (FPU_st0_ptr->sign == SIGN_NEG) {
+							if ((FPU_st0_ptr->exp >= EXP_BIAS) &&
+							    !((FPU_st0_ptr->sigh == 0x80000000) &&
+								(FPU_st0_ptr->sigl == 0))) {
+								/* st(0) holds
+								 * < -1.0 */
+								arith_invalid(st1_ptr);
+								pop();
+								return;
+							}
+#ifdef DENORM_OPERAND
+							if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+								return;
+#endif				/* DENORM_OPERAND */
+							st1_ptr->sign ^= SIGN_POS ^ SIGN_NEG;
+							pop();
+							return;
+						}
+#ifdef DENORM_OPERAND
+						if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+							return;
+#endif				/* DENORM_OPERAND */
+						pop();
+						return;
+					}
+					if (st1_tag == TW_NaN) {
+						real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
+						pop();
+						return;
+					}
+				} else
+					if (FPU_st0_tag == TW_NaN) {
+						real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
+						pop();
+						return;
+					} else
+						if (FPU_st0_tag == TW_Infinity) {
+							if (st1_tag == TW_NaN) {
+								real_2op_NaN(FPU_st0_ptr, st1_ptr, st1_ptr);
+								pop();
+								return;
+							} else
+								if ((FPU_st0_ptr->sign == SIGN_NEG) ||
+								    (st1_tag == TW_Zero)) {
+									arith_invalid(st1_ptr);	/* log(infinity) */
+									pop();
+									return;
+								}
+							/* st(1) must be valid
+							 * here. */
+
+#ifdef DENORM_OPERAND
+							if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+								return;
+#endif				/* DENORM_OPERAND */
+
+							/* The Manual says
+							 * that log(Infinity)
+							 * is invalid, but a
+							 * real 80486 sensibly
+							 * says that it is
+							 * o.k. */
+							{
+								char    sign = st1_ptr->sign;
+								reg_move(&CONST_INF, st1_ptr);
+								st1_ptr->sign = sign;
+							}
+							pop();
+							return;
+						}
+#ifdef PARANOID
+						else {
+							EXCEPTION(EX_INTERNAL | 0x117);
+						}
+#endif				/* PARANOID */
+}
+
+
+static void
+emu_fscale(void)
+{
+	FPU_REG *st1_ptr = &st(1);
+	char    st1_tag = st1_ptr->tag;
+	int     old_cw = control_word;
+
+	if (!((FPU_st0_tag ^ TW_Valid) | (st1_tag ^ TW_Valid))) {
+		long    scale;
+		FPU_REG tmp;
+
+#ifdef DENORM_OPERAND
+		if (((FPU_st0_ptr->exp <= EXP_UNDER) ||
+			(st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()))
+			return;
+#endif				/* DENORM_OPERAND */
+
+		if (st1_ptr->exp > EXP_BIAS + 30) {
+			/* 2^31 is far too large, would require 2^(2^30) or
+			 * 2^(-2^30) */
+			char    sign;
+
+			if (st1_ptr->sign == SIGN_POS) {
+				EXCEPTION(EX_Overflow);
+				sign = FPU_st0_ptr->sign;
+				reg_move(&CONST_INF, FPU_st0_ptr);
+				FPU_st0_ptr->sign = sign;
+			} else {
+				EXCEPTION(EX_Underflow);
+				sign = FPU_st0_ptr->sign;
+				reg_move(&CONST_Z, FPU_st0_ptr);
+				FPU_st0_ptr->sign = sign;
+			}
+			return;
+		}
+		control_word &= ~CW_RC;
+		control_word |= RC_CHOP;
+		reg_move(st1_ptr, &tmp);
+		round_to_int(&tmp);	/* This can never overflow here */
+		control_word = old_cw;
+		scale = st1_ptr->sign ? -tmp.sigl : tmp.sigl;
+		scale += FPU_st0_ptr->exp;
+		FPU_st0_ptr->exp = scale;
+
+		/* Use round_reg() to properly detect under/overflow etc */
+		round_reg(FPU_st0_ptr, 0, control_word);
+
+		return;
+	} else
+		if (FPU_st0_tag == TW_Valid) {
+			if (st1_tag == TW_Zero) {
+
+#ifdef DENORM_OPERAND
+				if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+					return;
+#endif				/* DENORM_OPERAND */
+
+				return;
+			}
+			if (st1_tag == TW_Infinity) {
+				char    sign = st1_ptr->sign;
+
+#ifdef DENORM_OPERAND
+				if ((FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+					return;
+#endif				/* DENORM_OPERAND */
+
+				if (sign == SIGN_POS) {
+					reg_move(&CONST_INF, FPU_st0_ptr);
+				} else
+					reg_move(&CONST_Z, FPU_st0_ptr);
+				FPU_st0_ptr->sign = sign;
+				return;
+			}
+			if (st1_tag == TW_NaN) {
+				real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
+				return;
+			}
+		} else
+			if (FPU_st0_tag == TW_Zero) {
+				if (st1_tag == TW_Valid) {
+
+#ifdef DENORM_OPERAND
+					if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+						return;
+#endif				/* DENORM_OPERAND */
+
+					return;
+				} else
+					if (st1_tag == TW_Zero) {
+						return;
+					} else
+						if (st1_tag == TW_Infinity) {
+							if (st1_ptr->sign == SIGN_NEG)
+								return;
+							else {
+								arith_invalid(FPU_st0_ptr);	/* Zero scaled by
+												 * +Infinity */
+								return;
+							}
+						} else
+							if (st1_tag == TW_NaN) {
+								real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
+								return;
+							}
+			} else
+				if (FPU_st0_tag == TW_Infinity) {
+					if (st1_tag == TW_Valid) {
+
+#ifdef DENORM_OPERAND
+						if ((st1_ptr->exp <= EXP_UNDER) && (denormal_operand()))
+							return;
+#endif				/* DENORM_OPERAND */
+
+						return;
+					}
+					if (((st1_tag == TW_Infinity) && (st1_ptr->sign == SIGN_POS))
+					    || (st1_tag == TW_Zero))
+						return;
+					else
+						if (st1_tag == TW_Infinity) {
+							arith_invalid(FPU_st0_ptr);	/* Infinity scaled by
+											 * -Infinity */
+							return;
+						} else
+							if (st1_tag == TW_NaN) {
+								real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
+								return;
+							}
+				} else
+					if (FPU_st0_tag == TW_NaN) {
+						if (st1_tag != TW_Empty) {
+							real_2op_NaN(FPU_st0_ptr, st1_ptr, FPU_st0_ptr);
+							return;
+						}
+					}
+#ifdef PARANOID
+	if (!((FPU_st0_tag == TW_Empty) || (st1_tag == TW_Empty))) {
+		EXCEPTION(EX_INTERNAL | 0x115);
+		return;
+	}
+#endif
+
+	/* At least one of st(0), st(1) must be empty */
+	stack_underflow();
+
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+static FUNC trig_table_a[] = {
+	f2xm1, fyl2x, fptan, fpatan, fxtract, fprem1, fdecstp, fincstp
+};
+
+void
+trig_a(void)
+{
+	(trig_table_a[FPU_rm]) ();
+}
+
+
+static FUNC trig_table_b[] =
+{
+	fprem, fyl2xp1, fsqrt_, fsincos, frndint_, emu_fscale, fsin, fcos
+};
+
+void
+trig_b(void)
+{
+	(trig_table_b[FPU_rm]) ();
+}