fpu_trig.c

linux 内核源代码

	  /* We do not really know if up or down */
	  set_precision_flag_up();
	  return 0;
	}
      else
	{
	  /* For a small arg, the result == the argument */
	  set_precision_flag_up();  /* Must be up. */
	  return 0;
	}
    }

  if ( tag == TAG_Zero )
    {
      setcc(0);
      return 0;
    }

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

  if ( tag == TW_Denormal )
    {
      if ( denormal_operand() < 0 )
	return 1;

      /* For a small arg, the result == the argument */
      /* Underflow may happen */
      FPU_to_exp16(st0_ptr, st0_ptr);
      
      tag = FPU_round(st0_ptr, 1, 0, FULL_PRECISION, arg_sign);

      FPU_settag0(tag);

      return 0;
    }
  else if ( tag == TW_Infinity )
    {
      /* The 80486 treats infinity as an invalid operand */
      arith_invalid(0);
      return 1;
    }
  else
    {
      single_arg_error(st0_ptr, tag);
      return 1;
    }
}


static int f_cos(FPU_REG *st0_ptr, u_char tag)
{
  u_char st0_sign;

  st0_sign = getsign(st0_ptr);

  if ( tag == TAG_Valid )
    {
      int q;

      if ( exponent(st0_ptr) > -40 )
	{
	  if ( (exponent(st0_ptr) < 0)
	      || ((exponent(st0_ptr) == 0)
		  && (significand(st0_ptr) <= 0xc90fdaa22168c234LL)) )
	    {
	      poly_cos(st0_ptr);

	      /* We do not really know if up or down */
	      set_precision_flag_down();
	  
	      return 0;
	    }
	  else if ( (q = trig_arg(st0_ptr, FCOS)) != -1 )
	    {
	      poly_sine(st0_ptr);

	      if ((q+1) & 2)
		changesign(st0_ptr);

	      /* We do not really know if up or down */
	      set_precision_flag_down();
	  
	      return 0;
	    }
	  else
	    {
	      /* Operand is out of range */
	      return 1;
	    }
	}
      else
	{
	denormal_arg:

	  setcc(0);
	  FPU_copy_to_reg0(&CONST_1, TAG_Valid);
#ifdef PECULIAR_486
	  set_precision_flag_down();  /* 80486 appears to do this. */
#else
	  set_precision_flag_up();  /* Must be up. */
#endif /* PECULIAR_486 */
	  return 0;
	}
    }
  else if ( tag == TAG_Zero )
    {
      FPU_copy_to_reg0(&CONST_1, TAG_Valid);
      setcc(0);
      return 0;
    }

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

  if ( tag == TW_Denormal )
    {
      if ( denormal_operand() < 0 )
	return 1;

      goto denormal_arg;
    }
  else if ( tag == TW_Infinity )
    {
      /* The 80486 treats infinity as an invalid operand */
      arith_invalid(0);
      return 1;
    }
  else
    {
      single_arg_error(st0_ptr, tag);  /* requires st0_ptr == &st(0) */
      return 1;
    }
}


static void fcos(FPU_REG *st0_ptr, u_char st0_tag)
{
  f_cos(st0_ptr, st0_tag);
}


static void fsincos(FPU_REG *st0_ptr, u_char st0_tag)
{
  FPU_REG *st_new_ptr;
  FPU_REG arg;
  u_char tag;

  /* 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_Special )
    tag = FPU_Special(st0_ptr);
  else
    tag = st0_tag;

  if ( tag == TW_NaN )
    {
      single_arg_2_error(st0_ptr, TW_NaN);
      return;
    }
  else if ( tag == TW_Infinity )
    {
      /* The 80486 treats infinity as an invalid operand */
      if ( arith_invalid(0) >= 0 )
	{
	  /* Masked response */
	  push();
	  arith_invalid(0);
	}
      return;
    }

  reg_copy(st0_ptr, &arg);
  if ( !fsin(st0_ptr, st0_tag) )
    {
      push();
      FPU_copy_to_reg0(&arg, st0_tag);
      f_cos(&st(0), st0_tag);
    }
  else
    {
      /* An error, so restore st(0) */
      FPU_copy_to_reg0(&arg, st0_tag);
    }
}


/*---------------------------------------------------------------------------*/
/* The following all require two arguments: st(0) and st(1) */

/* A lean, mean kernel for the fprem instructions. This relies upon
   the division and rounding to an integer in do_fprem giving an
   exact result. Because of this, rem_kernel() needs to deal only with
   the least significant 64 bits, the more significant bits of the
   result must be zero.
 */
static void rem_kernel(unsigned long long st0, unsigned long long *y,
		       unsigned long long st1,
		       unsigned long long q, int n)
{
  int dummy;
  unsigned long long x;

  x = st0 << n;

  /* Do the required multiplication and subtraction in the one operation */

  /* lsw x -= lsw st1 * lsw q */
  asm volatile ("mull %4; subl %%eax,%0; sbbl %%edx,%1"
		:"=m" (((unsigned *)&x)[0]), "=m" (((unsigned *)&x)[1]),
		"=a" (dummy)
		:"2" (((unsigned *)&st1)[0]), "m" (((unsigned *)&q)[0])
		:"%dx");
  /* msw x -= msw st1 * lsw q */
  asm volatile ("mull %3; subl %%eax,%0"
		:"=m" (((unsigned *)&x)[1]), "=a" (dummy)
		:"1" (((unsigned *)&st1)[1]), "m" (((unsigned *)&q)[0])
		:"%dx");
  /* msw x -= lsw st1 * msw q */
  asm volatile ("mull %3; subl %%eax,%0"
		:"=m" (((unsigned *)&x)[1]), "=a" (dummy)
		:"1" (((unsigned *)&st1)[0]), "m" (((unsigned *)&q)[1])
		:"%dx");

  *y = x;
}


/* Remainder of st(0) / st(1) */
/* This routine produces exact results, i.e. there is never any
   rounding or truncation, etc of the result. */
static void do_fprem(FPU_REG *st0_ptr, u_char st0_tag, int round)
{
  FPU_REG *st1_ptr = &st(1);
  u_char st1_tag = FPU_gettagi(1);

  if ( !((st0_tag ^ TAG_Valid) | (st1_tag ^ TAG_Valid)) )
    {
      FPU_REG tmp, st0, st1;
      u_char st0_sign, st1_sign;
      u_char tmptag;
      int tag;
      int old_cw;
      int expdif;
      long long q;
      unsigned short saved_status;
      int cc;

    fprem_valid:
      /* Convert registers for internal use. */
      st0_sign = FPU_to_exp16(st0_ptr, &st0);
      st1_sign = FPU_to_exp16(st1_ptr, &st1);
      expdif = exponent16(&st0) - exponent16(&st1);

      old_cw = control_word;
      cc = 0;

      /* We want the status following the denorm tests, but don't want
	 the status changed by the arithmetic operations. */
      saved_status = partial_status;
      control_word &= ~CW_RC;
      control_word |= RC_CHOP;

      if ( expdif < 64 )
	{
	  /* This should be the most common case */

	  if ( expdif > -2 )
	    {
	      u_char sign = st0_sign ^ st1_sign;
	      tag = FPU_u_div(&st0, &st1, &tmp,
			      PR_64_BITS | RC_CHOP | 0x3f,
			      sign);
	      setsign(&tmp, sign);

	      if ( exponent(&tmp) >= 0 )
		{
		  FPU_round_to_int(&tmp, tag);  /* Fortunately, this can't
						   overflow to 2^64 */
		  q = significand(&tmp);

		  rem_kernel(significand(&st0),
			     &significand(&tmp),
			     significand(&st1),
			     q, expdif);

		  setexponent16(&tmp, exponent16(&st1));
		}
	      else
		{
		  reg_copy(&st0, &tmp);
		  q = 0;
		}

	      if ( (round == RC_RND) && (tmp.sigh & 0xc0000000) )
		{
		  /* We may need to subtract st(1) once more,
		     to get a result <= 1/2 of st(1). */
		  unsigned long long x;
		  expdif = exponent16(&st1) - exponent16(&tmp);
		  if ( expdif <= 1 )
		    {
		      if ( expdif == 0 )
			x = significand(&st1) - significand(&tmp);
		      else /* expdif is 1 */
			x = (significand(&st1) << 1) - significand(&tmp);
		      if ( (x < significand(&tmp)) ||
			  /* or equi-distant (from 0 & st(1)) and q is odd */
			  ((x == significand(&tmp)) && (q & 1) ) )
			{
			  st0_sign = ! st0_sign;
			  significand(&tmp) = x;
			  q++;
			}
		    }
		}

	      if (q & 4) cc |= SW_C0;
	      if (q & 2) cc |= SW_C3;
	      if (q & 1) cc |= SW_C1;
	    }
	  else
	    {
	      control_word = old_cw;
	      setcc(0);
	      return;
	    }
	}
      else
	{
	  /* There is a large exponent difference ( >= 64 ) */
	  /* To make much sense, the code in this section should
	     be done at high precision. */
	  int exp_1, N;
	  u_char sign;

	  /* prevent overflow here */
	  /* N is 'a number between 32 and 63' (p26-113) */
	  reg_copy(&st0, &tmp);
	  tmptag = st0_tag;
	  N = (expdif & 0x0000001f) + 32;  /* This choice gives results
					      identical to an AMD 486 */
	  setexponent16(&tmp, N);
	  exp_1 = exponent16(&st1);
	  setexponent16(&st1, 0);
	  expdif -= N;

	  sign = getsign(&tmp) ^ st1_sign;
	  tag = FPU_u_div(&tmp, &st1, &tmp, PR_64_BITS | RC_CHOP | 0x3f,
			  sign);
	  setsign(&tmp, sign);

	  FPU_round_to_int(&tmp, tag);  /* Fortunately, this can't
					   overflow to 2^64 */

	  rem_kernel(significand(&st0),
		     &significand(&tmp),
		     significand(&st1),
		     significand(&tmp),
		     exponent(&tmp)
		     ); 
	  setexponent16(&tmp, exp_1 + expdif);

	  /* It is possible for the operation to be complete here.
	     What does the IEEE standard say? The Intel 80486 manual
	     implies that the operation will never be completed at this
	     point, and the behaviour of a real 80486 confirms this.
	   */
	  if ( !(tmp.sigh | tmp.sigl) )
	    {
	      /* The result is zero */
	      control_word = old_cw;
	      partial_status = saved_status;
	      FPU_copy_to_reg0(&CONST_Z, TAG_Zero);
	      setsign(&st0, st0_sign);
#ifdef PECULIAR_486
	      setcc(SW_C2);
#else
	      setcc(0);
#endif /* PECULIAR_486 */
	      return;
	    }
	  cc = SW_C2;
	}

      control_word = old_cw;
      partial_status = saved_status;
      tag = FPU_normalize_nuo(&tmp);
      reg_copy(&tmp, st0_ptr);

      /* The only condition to be looked for is underflow,
	 and it can occur here only if underflow is unmasked. */
      if ( (exponent16(&tmp) <= EXP_UNDER) && (tag != TAG_Zero)
	  && !(control_word & CW_Underflow) )
	{
	  setcc(cc);
	  tag = arith_underflow(st0_ptr);
	  setsign(st0_ptr, st0_sign);
	  FPU_settag0(tag);
	  return;
	}
      else if ( (exponent16(&tmp) > EXP_UNDER) || (tag == TAG_Zero) )
	{
	  stdexp(st0_ptr);
	  setsign(st0_ptr, st0_sign);
	}
      else
	{
	  tag = FPU_round(st0_ptr, 0, 0, FULL_PRECISION, st0_sign);
	}
      FPU_settag0(tag);
      setcc(cc);

      return;
    }

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

  if ( ((st0_tag == TAG_Valid) && (st1_tag == TW_Denormal))
	    || ((st0_tag == TW_Denormal) && (st1_tag == TAG_Valid))
	    || ((st0_tag == TW_Denormal) && (st1_tag == TW_Denormal)) )
    {
      if ( denormal_operand() < 0 )
	return;
      goto fprem_valid;
    }
  else if ( (st0_tag == TAG_Empty) || (st1_tag == TAG_Empty) )
    {
      FPU_stack_underflow();
      return;
    }
  else if ( st0_tag == TAG_Zero )
    {
      if ( st1_tag == TAG_Valid )
	{
	  setcc(0); return;
	}
      else if ( st1_tag == TW_Denormal )
	{
	  if ( denormal_operand() < 0 )
	    return;
	  setcc(0); return;
	}
      else if ( st1_tag == TAG_Zero )
	{ arith_invalid(0); return; } /* fprem(?,0) always invalid */
      else if ( st1_tag == TW_Infinity )
	{ setcc(0); return; }
    }
  else if ( (st0_tag == TAG_Valid) || (st0_tag == TW_Denormal) )
    {
      if ( st1_tag == TAG_Zero )
	{
	  arith_invalid(0); /* fprem(Valid,Zero) is invalid */
	  return;
	}
      else if ( st1_tag != TW_NaN )
	{
	  if ( ((st0_tag == TW_Denormal) || (st1_tag == TW_Denormal))
	       && (denormal_operand() < 0) )
	    return;

	  if ( st1_tag == TW_Infinity )
	    {
	      /* fprem(Valid,Infinity) is o.k. */
	      setcc(0); return;
	    }
	}
    }
  else if ( st0_tag == TW_Infinity )
    {
      if ( st1_tag != TW_NaN )
	{
	  arith_invalid(0); /* fprem(Infinity,?) is invalid */
	  return;
	}
    }

  /* One of the registers must contain a NaN if we got here. */

#ifdef PARANOID
  if ( (st0_tag != TW_NaN) && (st1_tag != TW_NaN) )
      EXCEPTION(EX_INTERNAL | 0x118);
#endif /* PARANOID */

  real_2op_NaN(st1_ptr, st1_tag, 0, st1_ptr);

}


/* ST(1) <- ST(1) * log ST;  pop ST */
static void fyl2x(FPU_REG *st0_ptr, u_char st0_tag)
{
  FPU_REG *st1_ptr = &st(1), exponent;
  u_char st1_tag = FPU_gettagi(1);
  u_char sign;
  int e, tag;

  clear_C1();

  if ( (st0_tag == TAG_Valid) && (st1_tag == TAG_Valid) )
    {
    both_valid:
      /* Both regs are Valid or Denormal */
      if ( signpositive(st0_ptr) )
	{
	  if ( st0_tag == TW_Denormal )
	    FPU_to_exp16(st0_ptr, st0_ptr);
	  else
	    /* Convert st(0) for internal use. */
	    setexponent16(st0_ptr, exponent(st0_ptr));

	  if ( (st0_ptr->sigh == 0x80000000) && (st0_ptr->sigl == 0) )
	    {
	      /* Special case. The result can be precise. */
	      u_char esign;
	      e = exponent16(st0_ptr);
	      if ( e >= 0 )
		{
		  exponent.sigh = e;
		  esign = SIGN_POS;
		}
	      else
		{
		  exponent.sigh = -e;
		  esign = SIGN_NEG;
		}
	      exponent.sigl = 0;
	      setexponent16(&exponent, 31);
	      tag = FPU_normalize_nuo(&exponent);
	      stdexp(&exponent);
	      setsign(&exponent, esign);
	      tag = FPU_mul(&exponent, tag, 1, FULL_PRECISION);
	      if ( tag >= 0 )
		FPU_settagi(1, tag);
	    }
	  else
	    {
	      /* The usual case */
	      sign = getsign(st1_ptr);
	      if ( st1_tag == TW_Denormal )
		FPU_to_exp16(st1_ptr, st1_ptr);
	      else
		/* Convert st(1) for internal use. */
		setexponent16(st1_ptr, exponent(st1_ptr));
	      poly_l2(st0_ptr, st1_ptr, sign);
	    }
	}
      else
	{
	  /* negative */
	  if ( arith_invalid(1) < 0 )
	    return;
	}

      FPU_pop();

      return;
    }

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

  if ( (st0_tag == TAG_Empty) || (st1_tag == TAG_Empty) )
    {
      FPU_stack_underflow_pop(1);
      return;
    }
  else if ( (st0_tag <= TW_Denormal) && (st1_tag <= TW_Denormal) )
    {
      if ( st0_tag == TAG_Zero )
	{
	  if ( st1_tag == TAG_Zero )
	    {
	      /* Both args zero is invalid */
	      if ( arith_invalid(1) < 0 )
		return;
	    }
	  else
	    {
	      u_char sign;
	      sign = getsign(st1_ptr)^SIGN_NEG;
	      if ( FPU_divide_by_zero(1, sign) < 0 )
		return;

	      setsign(st1_ptr, sign);
	    }
	}
      else if ( st1_tag == TAG_Zero )
	{
	  /* st(1) contains zero, st(0) valid <> 0 */
	  /* Zero is the valid answer */
	  sign = getsign(st1_ptr);
	  
	  if ( signnegative(st0_ptr) )
	    {
	      /* log(negative) */