fpu_trig.c

汇编源码大全 有各种汇编源码 希望对你有所帮助

/*---------------------------------------------------------------------------+
 |  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]    |
 |                                                                           |
 |                                                                           |
 +---------------------------------------------------------------------------*/

#include "fpu_system.h"
#include "exception.h"
#include "fpu_emu.h"
#include "status_w.h"
#include "control_w.h"
#include "reg_constant.h"	


static void rem_kernel(unsigned long long st0, unsigned long long *y,
		       unsigned long long st1,
		       unsigned long long q, int n);

#define BETTER_THAN_486

#define FCOS  4
#define FPTAN 1


/* Used only by fptan, fsin, fcos, and fsincos. */
/* This routine produces very accurate results, similar to
   using a value of pi with more than 128 bits precision. */
/* Limited measurements show no results worse than 64 bit precision
   except for the results for arguments close to 2^63, where the
   precision of the result sometimes degrades to about 63.9 bits */
static int trig_arg(FPU_REG *X, int even)
{
  FPU_REG tmp;
  unsigned long long q;
  int old_cw = control_word, saved_status = partial_status;

  if ( X->exp >= EXP_BIAS + 63 )
    {
      partial_status |= SW_C2;     /* Reduction incomplete. */
      return -1;
    }

  control_word &= ~CW_RC;
  control_word |= RC_CHOP;

  reg_div(X, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f);
  round_to_int(&tmp);  /* Fortunately, this can't overflow
			  to 2^64 */
  q = significand(&tmp);
  if ( q )
    {
      rem_kernel(significand(X),
		 &significand(&tmp),
		 significand(&CONST_PI2),
		 q, X->exp - CONST_PI2.exp);
      tmp.exp = CONST_PI2.exp;
      normalize(&tmp);
      reg_move(&tmp, X);
    }

  if ( even == FPTAN )
    {
      if ( ((X->exp >= EXP_BIAS) ||
	    ((X->exp == EXP_BIAS-1)
	     && (X->sigh >= 0xc90fdaa2))) ^ (q & 1) )
	even = FCOS;
      else
	even = 0;
    }

  if ( (even && !(q & 1)) || (!even && (q & 1)) )
    {
      reg_sub(&CONST_PI2, X, X, FULL_PRECISION);
#ifdef BETTER_THAN_486
      /* So far, the results are exact but based upon a 64 bit
	 precision approximation to pi/2. The technique used
	 now is equivalent to using an approximation to pi/2 which
	 is accurate to about 128 bits. */
      if ( (X->exp <= CONST_PI2extra.exp + 64) || (q > 1) )
	{
	  /* This code gives the effect of having p/2 to better than
	     128 bits precision. */
	  significand(&tmp) = q + 1;
	  tmp.exp = EXP_BIAS + 63;
	  tmp.tag = TW_Valid;
	  normalize(&tmp);
	  reg_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION);
	  reg_add(X, &tmp,  X, FULL_PRECISION);
	  if ( X->sign == SIGN_NEG )
	    {
	      /* CONST_PI2extra is negative, so the result of the addition
		 can be negative. This means that the argument is actually
		 in a different quadrant. The correction is always < pi/2,
		 so it can't overflow into yet another quadrant. */
	      X->sign = SIGN_POS;
	      q++;
	    }
	}
#endif BETTER_THAN_486
    }
#ifdef BETTER_THAN_486
  else
    {
      /* So far, the results are exact but based upon a 64 bit
	 precision approximation to pi/2. The technique used
	 now is equivalent to using an approximation to pi/2 which
	 is accurate to about 128 bits. */
      if ( ((q > 0) && (X->exp <= CONST_PI2extra.exp + 64)) || (q > 1) )
	{
	  /* This code gives the effect of having p/2 to better than
	     128 bits precision. */
	  significand(&tmp) = q;
	  tmp.exp = EXP_BIAS + 63;
	  tmp.tag = TW_Valid;
	  normalize(&tmp);
	  reg_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION);
	  reg_sub(X, &tmp, X, FULL_PRECISION);
	  if ( (X->exp == CONST_PI2.exp) &&
	      ((X->sigh > CONST_PI2.sigh)
	       || ((X->sigh == CONST_PI2.sigh)
		   && (X->sigl > CONST_PI2.sigl))) )
	    {
	      /* CONST_PI2extra is negative, so the result of the
		 subtraction can be larger than pi/2. This means
		 that the argument is actually in a different quadrant.
		 The correction is always < pi/2, so it can't overflow
		 into yet another quadrant. */
	      reg_sub(&CONST_PI, X, X, FULL_PRECISION);
	      q++;
	    }
	}
    }
#endif BETTER_THAN_486

  control_word = old_cw;
  partial_status = saved_status & ~SW_C2;     /* Reduction complete. */

  return (q & 3) | even;
}


/* Convert a long to register */
void convert_l2reg(long const *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);
	  if ( control_word & CW_Invalid )
	    FPU_st0_ptr->sigh |= 0x40000000;	  /* Convert to a QNaN */
	}
      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 single_arg_2_error(void)
{
  FPU_REG *st_new_ptr;

  switch ( FPU_st0_tag )
    {
    case TW_NaN:
      if ( !(FPU_st0_ptr->sigh & 0x40000000) )   /* Signaling ? */
	{
	  EXCEPTION(EX_Invalid);
	  if ( control_word & CW_Invalid )
	    {
	      /* The masked response */
	      /* Convert to a QNaN */
	      FPU_st0_ptr->sigh |= 0x40000000;
	      st_new_ptr = &st(-1);
	      push();
	      reg_move(&st(1), FPU_st0_ptr);
	    }
	}
      else
	{
	  /* A QNaN */
	  st_new_ptr = &st(-1);
	  push();
	  reg_move(&st(1), FPU_st0_ptr);
	}
      break;              /* return with a NaN in st(0) */
#ifdef PARANOID
    default:
      EXCEPTION(EX_INTERNAL|0x0112);
#endif PARANOID
    }
}


/*---------------------------------------------------------------------------*/

static void f2xm1(void)
{
  clear_C1();
  switch ( FPU_st0_tag )
    {
    case TW_Valid:
      {
	FPU_REG rv, tmp;

	if ( FPU_st0_ptr->exp >= 0 )
	  {
	    /* For an 80486 FPU, the result is undefined. */
	  }
	else if ( FPU_st0_ptr->exp >= -64 )
	  {
	    if ( FPU_st0_ptr->sign == SIGN_POS )
	      {
		/* poly_2xm1(x) requires 0 < x < 1. */
		poly_2xm1(FPU_st0_ptr, &rv);
		reg_mul(&rv, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
	      }
	    else
	      {
		/* poly_2xm1(x) doesn't handle negative numbers yet. */
		/* So we compute z=poly_2xm1(-x), and the answer is
		   then -z/(1+z) */
		FPU_st0_ptr->sign = SIGN_POS;
		poly_2xm1(FPU_st0_ptr, &rv);
		reg_mul(&rv, FPU_st0_ptr, &rv, FULL_PRECISION);
		reg_add(&rv, &CONST_1, &tmp, FULL_PRECISION);
		reg_div(&rv, &tmp, FPU_st0_ptr, FULL_PRECISION);
		FPU_st0_ptr->sign = SIGN_NEG;
	      }
	  }
	else
	  {
#ifdef DENORM_OPERAND
	    if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
	      return;
#endif DENORM_OPERAND
	    /* For very small arguments, this is accurate enough. */
	    reg_mul(&CONST_LN2, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
	  }
	set_precision_flag_up();
	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;

  /* Stack underflow has higher priority */
  if ( FPU_st0_tag == TW_Empty )
    {
      stack_underflow();  /* Puts a QNaN in st(0) */
      if ( control_word & CW_Invalid )
	{
	  st_new_ptr = &st(-1);
	  push();
	  stack_underflow();  /* Puts a QNaN in the new st(0) */
	}
      return;
    }

  if ( STACK_OVERFLOW )
    { stack_overflow(); return; }

  switch ( FPU_st0_tag )
    {
    case TW_Valid:

      if ( FPU_st0_ptr->exp > EXP_BIAS - 40 )
	{
	  FPU_st0_ptr->sign = SIGN_POS;
	  if ( (q = trig_arg(FPU_st0_ptr, FPTAN)) != -1 )
	    {
	      reg_div(FPU_st0_ptr, &CONST_PI2, FPU_st0_ptr,
		      FULL_PRECISION);
	      poly_tan(FPU_st0_ptr, FPU_st0_ptr, q & FCOS);
	      FPU_st0_ptr->sign = (q & 1) ^ arg_sign;
	    }
	  else
	    {
	      /* Operand is out of range */
	      FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
	      return;
	    }
	}
      else
	{
	  /* For a small arg, the result == the argument */
	  /* Underflow may happen */

	  if ( FPU_st0_ptr->exp <= EXP_UNDER )
	    {
#ifdef DENORM_OPERAND
	      if ( denormal_operand() )
		return;
#endif DENORM_OPERAND
	      /* A denormal result has been produced.
		 Precision must have been lost, this is always
		 an underflow. */
	      if ( arith_underflow(FPU_st0_ptr) )
		return;
	    }
	  else
	    set_precision_flag_up();  /* Must be up. */
	}
      push();
      reg_move(&CONST_1, FPU_st0_ptr);
      return;
      break;
    case TW_Infinity:
      /* The 80486 treats infinity as an invalid operand */
      arith_invalid(FPU_st0_ptr);
      if ( control_word & CW_Invalid )
	{
	  st_new_ptr = &st(-1);
	  push();
	  arith_invalid(FPU_st0_ptr);
	}
      return;
    case TW_Zero:
      push();
      reg_move(&CONST_1, FPU_st0_ptr);
      setcc(0);
      break;
    default:
      single_arg_2_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; }
  clear_C1();
  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;
      if ( divide_by_zero(SIGN_NEG, FPU_st0_ptr) )
	return;
      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 ( real_2op_NaN(FPU_st0_ptr, FPU_st0_ptr, FPU_st0_ptr) )
	return;
      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)
{
  clear_C1();
  top--;  /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
}

static void fincstp(void)
{
  clear_C1();
  top++;  /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
}


static void fsqrt_(void)
{
  clear_C1();
  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)
{
  int flags;

  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

      /* Fortunately, this can't overflow to 2^64 */
      if ( (flags = round_to_int(FPU_st0_ptr)) )
	set_precision_flag(flags);

      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 )
    {
      FPU_REG rv;
      int q;

      if ( FPU_st0_ptr->exp > EXP_BIAS - 40 )
	{
	  FPU_st0_ptr->sign = SIGN_POS;
	  if ( (q = trig_arg(FPU_st0_ptr, 0)) != -1 )
	    {
	      reg_div(FPU_st0_ptr, &CONST_PI2, FPU_st0_ptr, FULL_PRECISION);

	      poly_sine(FPU_st0_ptr, &rv);

	      if (q & 2)
		rv.sign ^= SIGN_POS ^ SIGN_NEG;
	      rv.sign ^= arg_sign;
	      reg_move(&rv, FPU_st0_ptr);

	      /* We do not really know if up or down */
	      set_precision_flag_up();
	      return;
	    }
	  else
	    {
	      /* Operand is out of range */
	      FPU_st0_ptr->sign = arg_sign;         /* restore st(0) */
	      return;
	    }
	}
      else
	{
	  /* For a small arg, the result == the argument */
	  /* Underflow may happen */

	  if ( FPU_st0_ptr->exp <= EXP_UNDER )
	    {
#ifdef DENORM_OPERAND
	      if ( denormal_operand() )
		return;
#endif DENORM_OPERAND
	      /* A denormal result has been produced.
		 Precision must have been lost, this is always
		 an underflow. */
	      arith_underflow(FPU_st0_ptr);