fpu_trig.c

linux 内核源代码

/*---------------------------------------------------------------------------+
 |  fpu_trig.c                                                               |
 |                                                                           |
 | Implementation of the FPU "transcendental" functions.                     |
 |                                                                           |
 | Copyright (C) 1992,1993,1994,1997,1999                                    |
 |                       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

/* 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 *st0_ptr, int even)
{
  FPU_REG tmp;
  u_char tmptag;
  unsigned long long q;
  int old_cw = control_word, saved_status = partial_status;
  int tag, st0_tag = TAG_Valid;

  if ( exponent(st0_ptr) >= 63 )
    {
      partial_status |= SW_C2;     /* Reduction incomplete. */
      return -1;
    }

  control_word &= ~CW_RC;
  control_word |= RC_CHOP;

  setpositive(st0_ptr);
  tag = FPU_u_div(st0_ptr, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f,
		  SIGN_POS);

  FPU_round_to_int(&tmp, tag);  /* Fortunately, this can't overflow
				   to 2^64 */
  q = significand(&tmp);
  if ( q )
    {
      rem_kernel(significand(st0_ptr),
		 &significand(&tmp),
		 significand(&CONST_PI2),
		 q, exponent(st0_ptr) - exponent(&CONST_PI2));
      setexponent16(&tmp, exponent(&CONST_PI2));
      st0_tag = FPU_normalize(&tmp);
      FPU_copy_to_reg0(&tmp, st0_tag);
    }

  if ( (even && !(q & 1)) || (!even && (q & 1)) )
    {
      st0_tag = FPU_sub(REV|LOADED|TAG_Valid, (int)&CONST_PI2, 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 ( (exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64) || (q > 1) )
	{
	  /* This code gives the effect of having pi/2 to better than
	     128 bits precision. */

	  significand(&tmp) = q + 1;
	  setexponent16(&tmp, 63);
	  FPU_normalize(&tmp);
	  tmptag =
	    FPU_u_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION, SIGN_POS,
		      exponent(&CONST_PI2extra) + exponent(&tmp));
	  setsign(&tmp, getsign(&CONST_PI2extra));
	  st0_tag = FPU_add(&tmp, tmptag, 0, FULL_PRECISION);
	  if ( signnegative(st0_ptr) )
	    {
	      /* 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. */
	      setpositive(st0_ptr);
	      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) && (exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64))
	   || (q > 1) )
	{
	  /* This code gives the effect of having p/2 to better than
	     128 bits precision. */

	  significand(&tmp) = q;
	  setexponent16(&tmp, 63);
	  FPU_normalize(&tmp);         /* This must return TAG_Valid */
	  tmptag = FPU_u_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION,
			     SIGN_POS,
			     exponent(&CONST_PI2extra) + exponent(&tmp));
	  setsign(&tmp, getsign(&CONST_PI2extra));
	  st0_tag = FPU_sub(LOADED|(tmptag & 0x0f), (int)&tmp,
			    FULL_PRECISION);
	  if ( (exponent(st0_ptr) == exponent(&CONST_PI2)) &&
	      ((st0_ptr->sigh > CONST_PI2.sigh)
	       || ((st0_ptr->sigh == CONST_PI2.sigh)
		   && (st0_ptr->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. */
	      st0_tag = FPU_sub(REV|LOADED|TAG_Valid, (int)&CONST_PI2,
				FULL_PRECISION);
	      q++;
	    }
	}
    }
#endif /* BETTER_THAN_486 */

  FPU_settag0(st0_tag);
  control_word = old_cw;
  partial_status = saved_status & ~SW_C2;     /* Reduction complete. */

  return (q & 3) | even;
}


/* Convert a long to register */
static void convert_l2reg(long const *arg, int deststnr)
{
  int tag;
  long num = *arg;
  u_char sign;
  FPU_REG *dest = &st(deststnr);

  if (num == 0)
    {
      FPU_copy_to_regi(&CONST_Z, TAG_Zero, deststnr);
      return;
    }

  if (num > 0)
    { sign = SIGN_POS; }
  else
    { num = -num; sign = SIGN_NEG; }

  dest->sigh = num;
  dest->sigl = 0;
  setexponent16(dest, 31);
  tag = FPU_normalize(dest);
  FPU_settagi(deststnr, tag);
  setsign(dest, sign);
  return;
}


static void single_arg_error(FPU_REG *st0_ptr, u_char st0_tag)
{
  if ( st0_tag == TAG_Empty )
    FPU_stack_underflow();  /* Puts a QNaN in st(0) */
  else if ( st0_tag == TW_NaN )
    real_1op_NaN(st0_ptr);       /* return with a NaN in st(0) */
#ifdef PARANOID
  else
    EXCEPTION(EX_INTERNAL|0x0112);
#endif /* PARANOID */
}


static void single_arg_2_error(FPU_REG *st0_ptr, u_char st0_tag)
{
  int isNaN;

  switch ( st0_tag )
    {
    case TW_NaN:
      isNaN = (exponent(st0_ptr) == EXP_OVER) && (st0_ptr->sigh & 0x80000000);
      if ( isNaN && !(st0_ptr->sigh & 0x40000000) )   /* Signaling ? */
	{
	  EXCEPTION(EX_Invalid);
	  if ( control_word & CW_Invalid )
	    {
	      /* The masked response */
	      /* Convert to a QNaN */
	      st0_ptr->sigh |= 0x40000000;
	      push();
	      FPU_copy_to_reg0(st0_ptr, TAG_Special);
	    }
	}
      else if ( isNaN )
	{
	  /* A QNaN */
	  push();
	  FPU_copy_to_reg0(st0_ptr, TAG_Special);
	}
      else
	{
	  /* pseudoNaN or other unsupported */
	  EXCEPTION(EX_Invalid);
	  if ( control_word & CW_Invalid )
	    {
	      /* The masked response */
	      FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
	      push();
	      FPU_copy_to_reg0(&CONST_QNaN, TAG_Special);
	    }
	}
      break;              /* return with a NaN in st(0) */
#ifdef PARANOID
    default:
      EXCEPTION(EX_INTERNAL|0x0112);
#endif /* PARANOID */
    }
}


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

static void f2xm1(FPU_REG *st0_ptr, u_char tag)
{
  FPU_REG a;

  clear_C1();

  if ( tag == TAG_Valid )
    {
      /* For an 80486 FPU, the result is undefined if the arg is >= 1.0 */
      if ( exponent(st0_ptr) < 0 )
	{
	denormal_arg:

	  FPU_to_exp16(st0_ptr, &a);

	  /* poly_2xm1(x) requires 0 < st(0) < 1. */
	  poly_2xm1(getsign(st0_ptr), &a, st0_ptr);
	}
      set_precision_flag_up();   /* 80486 appears to always do this */
      return;
    }

  if ( tag == TAG_Zero )
    return;

  if ( tag == TAG_Special )
    tag = FPU_Special(st0_ptr);

  switch ( tag )
    {
    case TW_Denormal:
      if ( denormal_operand() < 0 )
	return;
      goto denormal_arg;
    case TW_Infinity:
      if ( signnegative(st0_ptr) )
	{
	  /* -infinity gives -1 (p16-10) */
	  FPU_copy_to_reg0(&CONST_1, TAG_Valid);
	  setnegative(st0_ptr);
	}
      return;
    default:
      single_arg_error(st0_ptr, tag);
    }
}


static void fptan(FPU_REG *st0_ptr, u_char st0_tag)
{
  FPU_REG *st_new_ptr;
  int q;
  u_char arg_sign = getsign(st0_ptr);

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

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

  if ( st0_tag == TAG_Valid )
    {
      if ( exponent(st0_ptr) > -40 )
	{
	  if ( (q = trig_arg(st0_ptr, 0)) == -1 )
	    {
	      /* Operand is out of range */
	      return;
	    }

	  poly_tan(st0_ptr);
	  setsign(st0_ptr, (q & 1) ^ (arg_sign != 0));
	  set_precision_flag_up();  /* We do not really know if up or down */
	}
      else
	{
	  /* For a small arg, the result == the argument */
	  /* Underflow may happen */

	denormal_arg:

	  FPU_to_exp16(st0_ptr, st0_ptr);
      
	  st0_tag = FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);
	  FPU_settag0(st0_tag);
	}
      push();
      FPU_copy_to_reg0(&CONST_1, TAG_Valid);
      return;
    }

  if ( st0_tag == TAG_Zero )
    {
      push();
      FPU_copy_to_reg0(&CONST_1, TAG_Valid);
      setcc(0);
      return;
    }

  if ( st0_tag == TAG_Special )
    st0_tag = FPU_Special(st0_ptr);

  if ( st0_tag == TW_Denormal )
    {
      if ( denormal_operand() < 0 )
	return;

      goto denormal_arg;
    }

  if ( st0_tag == TW_Infinity )
    {
      /* The 80486 treats infinity as an invalid operand */
      if ( arith_invalid(0) >= 0 )
	{
	  st_new_ptr = &st(-1);
	  push();
	  arith_invalid(0);
	}
      return;
    }

  single_arg_2_error(st0_ptr, st0_tag);
}


static void fxtract(FPU_REG *st0_ptr, u_char st0_tag)
{
  FPU_REG *st_new_ptr;
  u_char sign;
  register FPU_REG *st1_ptr = st0_ptr;  /* anticipate */

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

  clear_C1();

  if ( st0_tag == TAG_Valid )
    {
      long e;

      push();
      sign = getsign(st1_ptr);
      reg_copy(st1_ptr, st_new_ptr);
      setexponent16(st_new_ptr, exponent(st_new_ptr));

    denormal_arg:

      e = exponent16(st_new_ptr);
      convert_l2reg(&e, 1);
      setexponentpos(st_new_ptr, 0);
      setsign(st_new_ptr, sign);
      FPU_settag0(TAG_Valid);       /* Needed if arg was a denormal */
      return;
    }
  else if ( st0_tag == TAG_Zero )
    {
      sign = getsign(st0_ptr);

      if ( FPU_divide_by_zero(0, SIGN_NEG) < 0 )
	return;

      push();
      FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
      setsign(st_new_ptr, sign);
      return;
    }

  if ( st0_tag == TAG_Special )
    st0_tag = FPU_Special(st0_ptr);

  if ( st0_tag == TW_Denormal )
    {
      if (denormal_operand() < 0 )
	return;

      push();
      sign = getsign(st1_ptr);
      FPU_to_exp16(st1_ptr, st_new_ptr);
      goto denormal_arg;
    }
  else if ( st0_tag == TW_Infinity )
    {
      sign = getsign(st0_ptr);
      setpositive(st0_ptr);
      push();
      FPU_copy_to_reg0(&CONST_INF, TAG_Special);
      setsign(st_new_ptr, sign);
      return;
    }
  else if ( st0_tag == TW_NaN )
    {
      if ( real_1op_NaN(st0_ptr) < 0 )
	return;

      push();
      FPU_copy_to_reg0(st0_ptr, TAG_Special);
      return;
    }
  else if ( st0_tag == TAG_Empty )
    {
      /* Is this the correct behaviour? */
      if ( control_word & EX_Invalid )
	{
	  FPU_stack_underflow();
	  push();
	  FPU_stack_underflow();
	}
      else
	EXCEPTION(EX_StackUnder);
    }
#ifdef PARANOID
  else
    EXCEPTION(EX_INTERNAL | 0x119);
#endif /* PARANOID */
}


static void fdecstp(void)
{
  clear_C1();
  top--;
}

static void fincstp(void)
{
  clear_C1();
  top++;
}


static void fsqrt_(FPU_REG *st0_ptr, u_char st0_tag)
{
  int expon;

  clear_C1();

  if ( st0_tag == TAG_Valid )
    {
      u_char tag;
      
      if (signnegative(st0_ptr))
	{
	  arith_invalid(0);  /* sqrt(negative) is invalid */
	  return;
	}

      /* make st(0) in  [1.0 .. 4.0) */
      expon = exponent(st0_ptr);

    denormal_arg:

      setexponent16(st0_ptr, (expon & 1));

      /* Do the computation, the sign of the result will be positive. */
      tag = wm_sqrt(st0_ptr, 0, 0, control_word, SIGN_POS);
      addexponent(st0_ptr, expon >> 1);
      FPU_settag0(tag);
      return;
    }

  if ( st0_tag == TAG_Zero )
    return;

  if ( st0_tag == TAG_Special )
    st0_tag = FPU_Special(st0_ptr);

  if ( st0_tag == TW_Infinity )
    {
      if ( signnegative(st0_ptr) )
	arith_invalid(0);  /* sqrt(-Infinity) is invalid */
      return;
    }
  else if ( st0_tag == TW_Denormal )
    {
      if (signnegative(st0_ptr))
	{
	  arith_invalid(0);  /* sqrt(negative) is invalid */
	  return;
	}

      if ( denormal_operand() < 0 )
	return;

      FPU_to_exp16(st0_ptr, st0_ptr);

      expon = exponent16(st0_ptr);

      goto denormal_arg;
    }

  single_arg_error(st0_ptr, st0_tag);

}


static void frndint_(FPU_REG *st0_ptr, u_char st0_tag)
{
  int flags, tag;

  if ( st0_tag == TAG_Valid )
    {
      u_char sign;

    denormal_arg:

      sign = getsign(st0_ptr);

      if (exponent(st0_ptr) > 63)
	return;

      if ( st0_tag == TW_Denormal )
	{
	  if (denormal_operand() < 0 )
	    return;
	}

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

      setexponent16(st0_ptr, 63);
      tag = FPU_normalize(st0_ptr);
      setsign(st0_ptr, sign);
      FPU_settag0(tag);
      return;
    }

  if ( st0_tag == TAG_Zero )
    return;

  if ( st0_tag == TAG_Special )
    st0_tag = FPU_Special(st0_ptr);

  if ( st0_tag == TW_Denormal )
    goto denormal_arg;
  else if ( st0_tag == TW_Infinity )
    return;
  else
    single_arg_error(st0_ptr, st0_tag);
}


static int fsin(FPU_REG *st0_ptr, u_char tag)
{
  u_char arg_sign = getsign(st0_ptr);

  if ( tag == TAG_Valid )
    {
      int q;

      if ( exponent(st0_ptr) > -40 )
	{
	  if ( (q = trig_arg(st0_ptr, 0)) == -1 )
	    {
	      /* Operand is out of range */
	      return 1;
	    }

	  poly_sine(st0_ptr);
	  
	  if (q & 2)
	    changesign(st0_ptr);

	  setsign(st0_ptr, getsign(st0_ptr) ^ arg_sign);