fpu_trig.c

freebsd v4.4内核源码

			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) {