fpu_trig.c

LINUX1.0内核源代码,学习LINUX编程的一定要看。

	      return;
	    }

	  set_precision_flag_up();  /* Must be up. */
	}
    }
  else if ( FPU_st0_tag == TW_Zero )
    {
      setcc(0);
      return;
    }
  else if ( FPU_st0_tag == TW_Infinity )
    {
      /* The 80486 treats infinity as an invalid operand */
      arith_invalid(FPU_st0_ptr);
      return;
    }
  else
    single_arg_error();
}


static int f_cos(FPU_REG *arg)
{
  char arg_sign = arg->sign;

  if ( arg->tag == TW_Valid )
    {
      FPU_REG rv;
      int q;

      if ( arg->exp > EXP_BIAS - 40 )
	{
	  arg->sign = SIGN_POS;
	  if ( (q = trig_arg(arg, FCOS)) != -1 )
	    {
	      reg_div(arg, &CONST_PI2, arg, FULL_PRECISION);
	      
	      poly_sine(arg, &rv);

	      if ((q+1) & 2)
		rv.sign ^= SIGN_POS ^ SIGN_NEG;
	      reg_move(&rv, arg);

	      /* We do not really know if up or down */
	      set_precision_flag_down();
	  
	      return 0;
	    }
	  else
	    {
	      /* Operand is out of range */
	      arg->sign = arg_sign;         /* restore st(0) */
	      return 1;
	    }
	}
      else
	{
#ifdef DENORM_OPERAND
	  if ( (arg->exp <= EXP_UNDER) && (denormal_operand()) )
	    return 1;
#endif DENORM_OPERAND

	  setcc(0);
	  reg_move(&CONST_1, arg);
#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 ( arg->tag == TW_Zero )
    {
      reg_move(&CONST_1, arg);
      setcc(0);
      return 0;
    }
  else if ( FPU_st0_tag == TW_Infinity )
    {
      /* The 80486 treats infinity as an invalid operand */
      arith_invalid(FPU_st0_ptr);
      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;

  /* 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; }

  if ( FPU_st0_tag == TW_NaN )
    {
      single_arg_2_error();
      return;
    }
  else if ( FPU_st0_tag == TW_Infinity )
    {
      /* The 80486 treats infinity as an invalid operand */
      if ( !arith_invalid(FPU_st0_ptr) )
	{
	  /* unmasked response */
	  push();
	  arith_invalid(FPU_st0_ptr);
	}
      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) */

/* 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)
{
  unsigned long long x;

  x = st0 << n;

  /* Do the required multiplication and subtraction in the one operation */
  asm volatile ("movl %2,%%eax; mull %4; subl %%eax,%0; sbbl %%edx,%1;
                 movl %3,%%eax; mull %4; subl %%eax,%1;
                 movl %2,%%eax; mull %5; subl %%eax,%1;"
		:"=m" (x), "=m" (((unsigned *)&x)[1])
		:"m" (st1),"m" (((unsigned *)&st1)[1]),
		 "m" (q),"m" (((unsigned *)&q)[1])
		:"%ax","%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(int round)
{
  FPU_REG *st1_ptr = &st(1);
  char st1_tag = st1_ptr->tag;
  char sign = FPU_st0_ptr->sign;

  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;
      long long q;
      unsigned short saved_status;
      int cc = 0;

#ifdef DENORM_OPERAND
      if ( ((FPU_st0_ptr->exp <= EXP_UNDER) ||
	    (st1_ptr->exp <= EXP_UNDER)) && (denormal_operand()) )
	return;
#endif DENORM_OPERAND
      
      /* 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 )
	    {
	      reg_div(FPU_st0_ptr, st1_ptr, &tmp, PR_64_BITS | RC_CHOP | 0x3f);

	      if ( tmp.exp >= EXP_BIAS )
		{
		  round_to_int(&tmp);  /* Fortunately, this can't overflow
					  to 2^64 */
		  q = significand(&tmp);

		  rem_kernel(significand(FPU_st0_ptr),
			     &significand(&tmp),
			     significand(st1_ptr),
			     q, expdif);

		  tmp.exp = st1_ptr->exp;
		}
	      else
		{
		  reg_move(FPU_st0_ptr, &tmp);
		  q = 0;
		}
	      tmp.sign = sign;

	      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 = st1_ptr->exp - tmp.exp;
		  if ( expdif <= 1 )
		    {
		      if ( expdif == 0 )
			x = significand(st1_ptr) - significand(&tmp);
		      else /* expdif is 1 */
			x = (significand(st1_ptr) << 1) - significand(&tmp);
		      if ( (x < significand(&tmp)) ||
			  /* or equi-distant (from 0 & st(1)) and q is odd */
			  ((x == significand(&tmp)) && (q & 1) ) )
			{
			  tmp.sign ^= (SIGN_POS^SIGN_NEG);
			  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;

	  /* prevent overflow here */
	  /* N is 'a number between 32 and 63' (p26-113) */
	  reg_move(FPU_st0_ptr, &tmp);
	  tmp.exp = EXP_BIAS + 56;
	  exp_1 = st1_ptr->exp;      st1_ptr->exp = EXP_BIAS;
	  expdif -= 56;

	  reg_div(&tmp, st1_ptr, &tmp, PR_64_BITS | RC_CHOP | 0x3f);
	  st1_ptr->exp = exp_1;

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

	  rem_kernel(significand(FPU_st0_ptr),
		     &significand(&tmp),
		     significand(st1_ptr),
		     significand(&tmp),
		     tmp.exp - EXP_BIAS
		     ); 
	  tmp.exp = exp_1 + expdif;
	  tmp.sign = sign;

	  /* 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;
	      reg_move(&CONST_Z, FPU_st0_ptr);
	      FPU_st0_ptr->sign = 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;
      normalize_nuo(&tmp);
      reg_move(&tmp, FPU_st0_ptr);
      setcc(cc);

      /* The only condition to be looked for is underflow,
	 and it can occur here only if underflow is unmasked. */
      if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (FPU_st0_ptr->tag != TW_Zero)
	  && !(control_word & CW_Underflow) )
	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(st1_ptr, FPU_st0_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;

  clear_C1();
  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 = partial_status;
	  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;
	  partial_status = 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 */
	  if ( !arith_invalid(st1_ptr) )
	    pop();
	  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) )
    {
      if ( !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 )
	{
	  if ( st1_tag == TW_Zero )
	    {
	      /* Both args zero is invalid */
	      if ( !arith_invalid(st1_ptr) )
		pop();
	    }
#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, st1_ptr);
	      pop();
	    }
#endif PECULIAR_486
	  else
	    {
	      if ( !divide_by_zero(st1_ptr->sign^SIGN_NEG^SIGN_POS, st1_ptr) )
		pop();
	    }
	  return;
	}
      else
	{
	  /* st(1) contains zero, st(0) valid <> 0 */
	  /* Zero is the valid answer */
	  char sign = st1_ptr->sign;

	  if ( FPU_st0_ptr->sign == SIGN_NEG )
	    {
	      /* log(negative) */
	      if ( !arith_invalid(st1_ptr) )
		pop();
	      return;
	    }

#ifdef DENORM_OPERAND
	  if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
	    return;
#endif DENORM_OPERAND

	  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) )
	{
	  /* log(-infinity) or 0*log(infinity) */
	  if ( !arith_invalid(st1_ptr) )
	    pop();
	  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) holds 1.0 */
	      /* infinity*log(1) */
	      if ( !arith_invalid(st1_ptr) )
		pop();
	      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 )
	{
	  /* This should be invalid, but a real 80486 is happy with it. */
#ifndef PECULIAR_486
	  if ( !divide_by_zero(st1_ptr->sign, st1_ptr) )
#endif PECULIAR_486
	    {
	      st1_ptr->sign ^= SIGN_NEG^SIGN_POS;
	      pop();
	    }
	}
      else
	{
	  /* log(negative) */
	  if ( !arith_invalid(st1_ptr) )