fpu_trig.c

freebsd v4.4内核源码

/*
 *  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   [e-mail]
 * 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: src/sys/gnu/i386/fpemul/fpu_trig.c,v 1.6.2.1 1999/09/05 08:09:42 peter Exp $
 *
 */


#include <sys/param.h>
#include <sys/proc.h>
#include <machine/cpu.h>
#include <machine/pcb.h>

#include <gnu/i386/fpemul/fpu_emu.h>
#include <gnu/i386/fpemul/fpu_system.h>
#include <gnu/i386/fpemul/exception.h>
#include <gnu/i386/fpemul/status_w.h>
#include <gnu/i386/fpemul/reg_constant.h>
#include <gnu/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