sim-fpu.c

这个是LINUX下的GDB调度工具的源码

  }
}


INLINE_SIM_FPU (int)
sim_fpu_mul (sim_fpu *f,
	     const sim_fpu *l,
	     const sim_fpu *r)
{
  if (sim_fpu_is_snan (l))
    {
      *f = *l;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_snan (r))
    {
      *f = *r;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_qnan (l))
    {
      *f = *l;
      return 0;
    }
  if (sim_fpu_is_qnan (r))
    {
      *f = *r;
      return 0;
    }
  if (sim_fpu_is_infinity (l))
    {
      if (sim_fpu_is_zero (r))
	{
	  *f = sim_fpu_qnan;
	  return sim_fpu_status_invalid_imz;
	}
      *f = *l;
      f->sign = l->sign ^ r->sign;
      return 0;
    }
  if (sim_fpu_is_infinity (r))
    {
      if (sim_fpu_is_zero (l))
	{
	  *f = sim_fpu_qnan;
	  return sim_fpu_status_invalid_imz;
	}
      *f = *r;
      f->sign = l->sign ^ r->sign;
      return 0;
    }
  if (sim_fpu_is_zero (l) || sim_fpu_is_zero (r))
    {
      *f = sim_fpu_zero;
      f->sign = l->sign ^ r->sign;
      return 0;
    }
  /* Calculate the mantissa by multiplying both 64bit numbers to get a
     128 bit number */
  {
    unsigned64 low;
    unsigned64 high;
    unsigned64 nl = l->fraction & 0xffffffff;
    unsigned64 nh = l->fraction >> 32;
    unsigned64 ml = r->fraction & 0xffffffff;
    unsigned64 mh = r->fraction >>32;
    unsigned64 pp_ll = ml * nl;
    unsigned64 pp_hl = mh * nl;
    unsigned64 pp_lh = ml * nh;
    unsigned64 pp_hh = mh * nh;
    unsigned64 res2 = 0;
    unsigned64 res0 = 0;
    unsigned64 ps_hh__ = pp_hl + pp_lh;
    if (ps_hh__ < pp_hl)
      res2 += UNSIGNED64 (0x100000000);
    pp_hl = (ps_hh__ << 32) & UNSIGNED64 (0xffffffff00000000);
    res0 = pp_ll + pp_hl;
    if (res0 < pp_ll)
      res2++;
    res2 += ((ps_hh__ >> 32) & 0xffffffff) + pp_hh;
    high = res2;
    low = res0;
    
    f->normal_exp = l->normal_exp + r->normal_exp;
    f->sign = l->sign ^ r->sign;
    f->class = sim_fpu_class_number;

    /* Input is bounded by [1,2)   ;   [2^60,2^61)
       Output is bounded by [1,4)  ;   [2^120,2^122) */
 
    /* Adjust the exponent according to where the decimal point ended
       up in the high 64 bit word.  In the source the decimal point
       was at NR_FRAC_GUARD. */
    f->normal_exp += NR_FRAC_GUARD + 64 - (NR_FRAC_GUARD * 2);

    /* The high word is bounded according to the above.  Consequently
       it has never overflowed into IMPLICIT_2. */
    ASSERT (high < LSBIT64 (((NR_FRAC_GUARD + 1) * 2) - 64));
    ASSERT (high >= LSBIT64 ((NR_FRAC_GUARD * 2) - 64));
    ASSERT (LSBIT64 (((NR_FRAC_GUARD + 1) * 2) - 64) < IMPLICIT_1);

    /* normalize */
    do
      {
	f->normal_exp--;
	high <<= 1;
	if (low & LSBIT64 (63))
	  high |= 1;
	low <<= 1;
      }
    while (high < IMPLICIT_1);

    ASSERT (high >= IMPLICIT_1 && high < IMPLICIT_2);
    if (low != 0)
      {
	f->fraction = (high | 1); /* sticky */
	return sim_fpu_status_inexact;
      }
    else
      {
	f->fraction = high;
	return 0;
      }
    return 0;
  }
}

INLINE_SIM_FPU (int)
sim_fpu_div (sim_fpu *f,
	     const sim_fpu *l,
	     const sim_fpu *r)
{
  if (sim_fpu_is_snan (l))
    {
      *f = *l;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_snan (r))
    {
      *f = *r;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_qnan (l))
    {
      *f = *l;
      f->class = sim_fpu_class_qnan;
      return 0;
    }
  if (sim_fpu_is_qnan (r))
    {
      *f = *r;
      f->class = sim_fpu_class_qnan;
      return 0;
    }
  if (sim_fpu_is_infinity (l))
    {
      if (sim_fpu_is_infinity (r))
	{
	  *f = sim_fpu_qnan;
	  return sim_fpu_status_invalid_idi;
	}
      else
	{
	  *f = *l;
	  f->sign = l->sign ^ r->sign;
	  return 0;
	}
    }
  if (sim_fpu_is_zero (l))
    {
      if (sim_fpu_is_zero (r))
	{
	  *f = sim_fpu_qnan;
	  return sim_fpu_status_invalid_zdz;
	}
      else
	{
	  *f = *l;
	  f->sign = l->sign ^ r->sign;
	  return 0;
	}
    }
  if (sim_fpu_is_infinity (r))
    {
      *f = sim_fpu_zero;
      f->sign = l->sign ^ r->sign;
      return 0;
    }
  if (sim_fpu_is_zero (r))
    {
      f->class = sim_fpu_class_infinity;
      f->sign = l->sign ^ r->sign;
      return sim_fpu_status_invalid_div0;
    }

  /* Calculate the mantissa by multiplying both 64bit numbers to get a
     128 bit number */
  {
    /* quotient =  ( ( numerator / denominator)
                      x 2^(numerator exponent -  denominator exponent)
     */
    unsigned64 numerator;
    unsigned64 denominator;
    unsigned64 quotient;
    unsigned64 bit;

    f->class = sim_fpu_class_number;
    f->sign = l->sign ^ r->sign;
    f->normal_exp = l->normal_exp - r->normal_exp;

    numerator = l->fraction;
    denominator = r->fraction;

    /* Fraction will be less than 1.0 */
    if (numerator < denominator)
      {
	numerator <<= 1;
	f->normal_exp--;
      }
    ASSERT (numerator >= denominator);
    
    /* Gain extra precision, already used one spare bit */
    numerator <<=    NR_SPARE;
    denominator <<=  NR_SPARE;

    /* Does divide one bit at a time.  Optimize???  */
    quotient = 0;
    bit = (IMPLICIT_1 << NR_SPARE);
    while (bit)
      {
	if (numerator >= denominator)
	  {
	    quotient |= bit;
	    numerator -= denominator;
	  }
	bit >>= 1;
	numerator <<= 1;
      }

    /* discard (but save) the extra bits */
    if ((quotient & LSMASK64 (NR_SPARE -1, 0)))
      quotient = (quotient >> NR_SPARE) | 1;
    else
      quotient = (quotient >> NR_SPARE);

    f->fraction = quotient;
    ASSERT (f->fraction >= IMPLICIT_1 && f->fraction < IMPLICIT_2);
    if (numerator != 0)
      {
	f->fraction |= 1; /* stick remaining bits */
	return sim_fpu_status_inexact;
      }
    else
      return 0;
  }
}


INLINE_SIM_FPU (int)
sim_fpu_max (sim_fpu *f,
	     const sim_fpu *l,
	     const sim_fpu *r)
{
  if (sim_fpu_is_snan (l))
    {
      *f = *l;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_snan (r))
    {
      *f = *r;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_qnan (l))
    {
      *f = *l;
      return 0;
    }
  if (sim_fpu_is_qnan (r))
    {
      *f = *r;
      return 0;
    }
  if (sim_fpu_is_infinity (l))
    {
      if (sim_fpu_is_infinity (r)
	  && l->sign == r->sign)
	{
	  *f = sim_fpu_qnan;
	  return sim_fpu_status_invalid_isi;
	}
      if (l->sign)
	*f = *r; /* -inf < anything */
      else
	*f = *l; /* +inf > anthing */
      return 0;
    }
  if (sim_fpu_is_infinity (r))
    {
      if (r->sign)
	*f = *l; /* anything > -inf */
      else
	*f = *r; /* anthing < +inf */
      return 0;
    }
  if (l->sign > r->sign)
    {
      *f = *r; /* -ve < +ve */
      return 0;
    }
  if (l->sign < r->sign)
    {
      *f = *l; /* +ve > -ve */
      return 0;
    }
  ASSERT (l->sign == r->sign);
  if (l->normal_exp > r->normal_exp
      || (l->normal_exp == r->normal_exp && 
	  l->fraction > r->fraction))
    {
      /* |l| > |r| */
      if (l->sign)
	*f = *r; /* -ve < -ve */
      else
	*f = *l; /* +ve > +ve */
      return 0;
    }
  else
    {
      /* |l| <= |r| */
      if (l->sign)
	*f = *l; /* -ve > -ve */
      else
	*f = *r; /* +ve < +ve */
      return 0;
    }
}


INLINE_SIM_FPU (int)
sim_fpu_min (sim_fpu *f,
	     const sim_fpu *l,
	     const sim_fpu *r)
{
  if (sim_fpu_is_snan (l))
    {
      *f = *l;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_snan (r))
    {
      *f = *r;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_qnan (l))
    {
      *f = *l;
      return 0;
    }
  if (sim_fpu_is_qnan (r))
    {
      *f = *r;
      return 0;
    }
  if (sim_fpu_is_infinity (l))
    {
      if (sim_fpu_is_infinity (r)
	  && l->sign == r->sign)
	{
	  *f = sim_fpu_qnan;
	  return sim_fpu_status_invalid_isi;
	}
      if (l->sign)
	*f = *l; /* -inf < anything */
      else
	*f = *r; /* +inf > anthing */
      return 0;
    }
  if (sim_fpu_is_infinity (r))
    {
      if (r->sign)
	*f = *r; /* anything > -inf */
      else
	*f = *l; /* anything < +inf */
      return 0;
    }
  if (l->sign > r->sign)
    {
      *f = *l; /* -ve < +ve */
      return 0;
    }
  if (l->sign < r->sign)
    {
      *f = *r; /* +ve > -ve */
      return 0;
    }
  ASSERT (l->sign == r->sign);
  if (l->normal_exp > r->normal_exp
      || (l->normal_exp == r->normal_exp && 
	  l->fraction > r->fraction))
    {
      /* |l| > |r| */
      if (l->sign)
	*f = *l; /* -ve < -ve */
      else
	*f = *r; /* +ve > +ve */
      return 0;
    }
  else
    {
      /* |l| <= |r| */
      if (l->sign)
	*f = *r; /* -ve > -ve */
      else
	*f = *l; /* +ve < +ve */
      return 0;
    }
}


INLINE_SIM_FPU (int)
sim_fpu_neg (sim_fpu *f,
	     const sim_fpu *r)
{
  if (sim_fpu_is_snan (r))
    {
      *f = *r;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_qnan (r))
    {
      *f = *r;
      return 0;
    }
  *f = *r;
  f->sign = !r->sign;
  return 0;
}


INLINE_SIM_FPU (int)
sim_fpu_abs (sim_fpu *f,
	     const sim_fpu *r)
{
  if (sim_fpu_is_snan (r))
    {
      *f = *r;
      f->class = sim_fpu_class_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_qnan (r))
    {
      *f = *r;
      return 0;
    }
  *f = *r;
  f->sign = 0;
  return 0;
}


INLINE_SIM_FPU (int)
sim_fpu_inv (sim_fpu *f,
	     const sim_fpu *r)
{
  return sim_fpu_div (f, &sim_fpu_one, r);
}


INLINE_SIM_FPU (int)
sim_fpu_sqrt (sim_fpu *f,
	      const sim_fpu *r)
{
  if (sim_fpu_is_snan (r))
    {
      *f = sim_fpu_qnan;
      return sim_fpu_status_invalid_snan;
    }
  if (sim_fpu_is_qnan (r))
    {
      *f = sim_fpu_qnan;
      return 0;
    }
  if (sim_fpu_is_zero (r))
    {
      f->class = sim_fpu_class_zero;
      f->sign = r->sign;
      f->normal_exp = 0;
      return 0;
    }
  if (sim_fpu_is_infinity (r))
    {
      if (r->sign)
	{
	  *f = sim_fpu_qnan;
	  return sim_fpu_status_invalid_sqrt;
	}
      else
	{
	  f->class = sim_fpu_class_infinity;
	  f->sign = 0;
	  f->sign = 0;
	  return 0;
	}
    }
  if (r->sign)
    {
      *f = sim_fpu_qnan;
      return sim_fpu_status_invalid_sqrt;
    }

  /* @(#)e_sqrt.c 5.1 93/09/24 */
  /*
   * ====================================================
   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   *
   * Developed at SunPro, a Sun Microsystems, Inc. business.
   * Permission to use, copy, modify, and distribute this
   * software is freely granted, provided that this notice 
   * is preserved.
   * ====================================================
   */
  
  /* __ieee754_sqrt(x)
   * Return correctly rounded sqrt.
   *           ------------------------------------------
   *           |  Use the hardware sqrt if you have one |
   *           ------------------------------------------
   * Method: 
   *   Bit by bit method using integer arithmetic. (Slow, but portable) 
   *   1. Normalization
   *	Scale x to y in [1,4) with even powers of 2: 
   *	find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
   *		sqrt(x) = 2^k * sqrt(y)
   -
   - Since:
   -   sqrt ( x*2^(2m) )     = sqrt(x).2^m    ; m even
   -   sqrt ( x*2^(2m + 1) ) = sqrt(2.x).2^m  ; m odd
   - Define:
   -   y = ((m even) ? x : 2.x)
   - Then:
   -   y in [1, 4)                            ; [IMPLICIT_1,IMPLICIT_4)
   - And:
   -   sqrt (y) in [1, 2)                     ; [IMPLICIT_1,IMPLICIT_2)
   -
   *   2. Bit by bit computation
   *	Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
   *	     i							 0
   *                                     i+1         2
   *	    s  = 2*q , and	y  =  2   * ( y - q  ).		(1)
   *	     i      i            i                 i
   *                                                        
   *	To compute q    from q , one checks whether 
   *		    i+1       i                       
   *
   *			      -(i+1) 2
   *			(q + 2      ) <= y.			(2)
   *     			  i
   *							      -(i+1)
   *	If (2) is false, then q   = q ; otherwise q   = q  + 2      .
   *		 	       i+1   i             i+1   i
   *
   *	With some algebric manipulation, it is not difficult to see
   *	that (2) is equivalent to 
   *                             -(i+1)
   *			s  +  2       <= y			(3)
   *			 i                i
   *
   *	The advantage of (3) is that s  and y  can be computed by 
   *				      i      i
   *	the following recurrence formula:
   *	    if (3) is false
   *
   *	    s     =  s  ,	y    = y   ;			(4)
   *	     i+1      i		 i+1    i
   *
   -
   -                      NOTE: y    = 2*y
   -                             i+1      i
   -
   *	    otherwise,
   *                       -i                      -(i+1)
   *	    s	  =  s  + 2  ,  y    = y  -  s  - 2  		(5)
   *         i+1      i          i+1    i     i
   *				
   -
   -                                                   -(i+1)
   -                      NOTE: y    = 2 (y  -  s  -  2      ) 		
   -                             i+1       i     i
   -
   *	One may easily use induction to prove (4) and (5). 
   *	Note. Since the left hand side of (3) contain only i+2 bits,
   *	      it does not necessary to do a full (53-bit) comparison 
   *	      in (3).
   *   3. Final rounding
   *	After generating the 53 bits result, we compute one more bit.
   *	Together with the remainder, we can decide whether the
   *	result is exact, bigger than 1/2ulp, or less than 1/2ulp
   *	(it will never equal to 1/2ulp).
   *	The rounding mode can be detected by checking whether
   *	huge + tiny is equal to huge, and whether huge - tiny is
   *	equal to huge for some floating point number "huge" and "tiny".
   *		
   * Special cases:
   *	sqrt(+-0) = +-0 	... exact
   *	sqrt(inf) = inf
   *	sqrt(-ve) = NaN		... with invalid signal
   *	sqrt(NaN) = NaN		... with invalid signal for signaling NaN
   *
   * Other methods : see the appended file at the end of the program below.
   *---------------
   */

  {
    /* generate sqrt(x) bit by bit */
    unsigned64 y;
    unsigned64 q;
    unsigned64 s;
    unsigned64 b;

    f->class = sim_fpu_class_number;
    f->sign = 0;
    y = r->fraction;
    f->normal_exp = (r->normal_exp >> 1);	/* exp = [exp/2] */

    /* odd exp, double x to make it even */
    ASSERT (y >= IMPLICIT_1 && y < IMPLICIT_4);
    if ((r->normal_exp & 1))
      {
	y += y;
      }
    ASSERT (y >= IMPLICIT_1 && y < (IMPLICIT_2 << 1));