real.c

gcc库的原代码,对编程有很大帮助.

  int i;

  x += M;			/* point to significand area */

  bits = 0;
  for (i = M; i < NI; i++)
    {
      if (*x & 1)
	bits |= 1;
      *x >>= 1;
      if (bits & 2)
	*x |= 0x8000;
      bits <<= 1;
      ++x;
    }
}

/* Shift significand of exploded e-type X up by 1 bit.  */

static void 
eshup1 (x)
     register unsigned EMUSHORT *x;
{
  register unsigned EMUSHORT bits;
  int i;

  x += NI - 1;
  bits = 0;

  for (i = M; i < NI; i++)
    {
      if (*x & 0x8000)
	bits |= 1;
      *x <<= 1;
      if (bits & 2)
	*x |= 1;
      bits <<= 1;
      --x;
    }
}


/* Shift significand of exploded e-type X down by 8 bits.  */

static void 
eshdn8 (x)
     register unsigned EMUSHORT *x;
{
  register unsigned EMUSHORT newbyt, oldbyt;
  int i;

  x += M;
  oldbyt = 0;
  for (i = M; i < NI; i++)
    {
      newbyt = *x << 8;
      *x >>= 8;
      *x |= oldbyt;
      oldbyt = newbyt;
      ++x;
    }
}

/* Shift significand of exploded e-type X up by 8 bits.  */

static void 
eshup8 (x)
     register unsigned EMUSHORT *x;
{
  int i;
  register unsigned EMUSHORT newbyt, oldbyt;

  x += NI - 1;
  oldbyt = 0;

  for (i = M; i < NI; i++)
    {
      newbyt = *x >> 8;
      *x <<= 8;
      *x |= oldbyt;
      oldbyt = newbyt;
      --x;
    }
}

/* Shift significand of exploded e-type X up by 16 bits.  */

static void 
eshup6 (x)
     register unsigned EMUSHORT *x;
{
  int i;
  register unsigned EMUSHORT *p;

  p = x + M;
  x += M + 1;

  for (i = M; i < NI - 1; i++)
    *p++ = *x++;

  *p = 0;
}

/* Shift significand of exploded e-type X down by 16 bits.  */

static void 
eshdn6 (x)
     register unsigned EMUSHORT *x;
{
  int i;
  register unsigned EMUSHORT *p;

  x += NI - 1;
  p = x + 1;

  for (i = M; i < NI - 1; i++)
    *(--p) = *(--x);

  *(--p) = 0;
}

/* Add significands of exploded e-type X and Y.  X + Y replaces Y.  */

static void 
eaddm (x, y)
     unsigned EMUSHORT *x, *y;
{
  register unsigned EMULONG a;
  int i;
  unsigned int carry;

  x += NI - 1;
  y += NI - 1;
  carry = 0;
  for (i = M; i < NI; i++)
    {
      a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
      if (a & 0x10000)
	carry = 1;
      else
	carry = 0;
      *y = (unsigned EMUSHORT) a;
      --x;
      --y;
    }
}

/* Subtract significands of exploded e-type X and Y.  Y - X replaces Y.  */

static void 
esubm (x, y)
     unsigned EMUSHORT *x, *y;
{
  unsigned EMULONG a;
  int i;
  unsigned int carry;

  x += NI - 1;
  y += NI - 1;
  carry = 0;
  for (i = M; i < NI; i++)
    {
      a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
      if (a & 0x10000)
	carry = 1;
      else
	carry = 0;
      *y = (unsigned EMUSHORT) a;
      --x;
      --y;
    }
}


static unsigned EMUSHORT equot[NI];


#if 0
/* Radix 2 shift-and-add versions of multiply and divide  */


/* Divide significands */

int 
edivm (den, num)
     unsigned EMUSHORT den[], num[];
{
  int i;
  register unsigned EMUSHORT *p, *q;
  unsigned EMUSHORT j;

  p = &equot[0];
  *p++ = num[0];
  *p++ = num[1];

  for (i = M; i < NI; i++)
    {
      *p++ = 0;
    }

  /* Use faster compare and subtraction if denominator has only 15 bits of
     significance.  */

  p = &den[M + 2];
  if (*p++ == 0)
    {
      for (i = M + 3; i < NI; i++)
	{
	  if (*p++ != 0)
	    goto fulldiv;
	}
      if ((den[M + 1] & 1) != 0)
	goto fulldiv;
      eshdn1 (num);
      eshdn1 (den);

      p = &den[M + 1];
      q = &num[M + 1];

      for (i = 0; i < NBITS + 2; i++)
	{
	  if (*p <= *q)
	    {
	      *q -= *p;
	      j = 1;
	    }
	  else
	    {
	      j = 0;
	    }
	  eshup1 (equot);
	  equot[NI - 2] |= j;
	  eshup1 (num);
	}
      goto divdon;
    }

  /* The number of quotient bits to calculate is NBITS + 1 scaling guard
     bit + 1 roundoff bit.  */

 fulldiv:

  p = &equot[NI - 2];
  for (i = 0; i < NBITS + 2; i++)
    {
      if (ecmpm (den, num) <= 0)
	{
	  esubm (den, num);
	  j = 1;		/* quotient bit = 1 */
	}
      else
	j = 0;
      eshup1 (equot);
      *p |= j;
      eshup1 (num);
    }

 divdon:

  eshdn1 (equot);
  eshdn1 (equot);

  /* test for nonzero remainder after roundoff bit */
  p = &num[M];
  j = 0;
  for (i = M; i < NI; i++)
    {
      j |= *p++;
    }
  if (j)
    j = 1;


  for (i = 0; i < NI; i++)
    num[i] = equot[i];
  return ((int) j);
}


/* Multiply significands */
int 
emulm (a, b)
     unsigned EMUSHORT a[], b[];
{
  unsigned EMUSHORT *p, *q;
  int i, j, k;

  equot[0] = b[0];
  equot[1] = b[1];
  for (i = M; i < NI; i++)
    equot[i] = 0;

  p = &a[NI - 2];
  k = NBITS;
  while (*p == 0)		/* significand is not supposed to be zero */
    {
      eshdn6 (a);
      k -= 16;
    }
  if ((*p & 0xff) == 0)
    {
      eshdn8 (a);
      k -= 8;
    }

  q = &equot[NI - 1];
  j = 0;
  for (i = 0; i < k; i++)
    {
      if (*p & 1)
	eaddm (b, equot);
      /* remember if there were any nonzero bits shifted out */
      if (*q & 1)
	j |= 1;
      eshdn1 (a);
      eshdn1 (equot);
    }

  for (i = 0; i < NI; i++)
    b[i] = equot[i];

  /* return flag for lost nonzero bits */
  return (j);
}

#else

/* Radix 65536 versions of multiply and divide.  */

/* Multiply significand of e-type number B
   by 16-bit quantity A, return e-type result to C. */

static void
m16m (a, b, c)
     unsigned int a;
     unsigned EMUSHORT b[], c[];
{
  register unsigned EMUSHORT *pp;
  register unsigned EMULONG carry;
  unsigned EMUSHORT *ps;
  unsigned EMUSHORT p[NI];
  unsigned EMULONG aa, m;
  int i;

  aa = a;
  pp = &p[NI-2];
  *pp++ = 0;
  *pp = 0;
  ps = &b[NI-1];

  for (i=M+1; i<NI; i++)
    {
      if (*ps == 0)
	{
	  --ps;
	  --pp;
	  *(pp-1) = 0;
	}
      else
	{
	  m = (unsigned EMULONG) aa * *ps--;
	  carry = (m & 0xffff) + *pp;
	  *pp-- = (unsigned EMUSHORT)carry;
	  carry = (carry >> 16) + (m >> 16) + *pp;
	  *pp = (unsigned EMUSHORT)carry;
	  *(pp-1) = carry >> 16;
	}
    }
  for (i=M; i<NI; i++)
    c[i] = p[i];
}

/* Divide significands of exploded e-types NUM / DEN.  Neither the
   numerator NUM nor the denominator DEN is permitted to have its high guard
   word nonzero.  */

static int
edivm (den, num)
     unsigned EMUSHORT den[], num[];
{
  int i;
  register unsigned EMUSHORT *p;
  unsigned EMULONG tnum;
  unsigned EMUSHORT j, tdenm, tquot;
  unsigned EMUSHORT tprod[NI+1];

  p = &equot[0];
  *p++ = num[0];
  *p++ = num[1];

  for (i=M; i<NI; i++)
    {
      *p++ = 0;
    }
  eshdn1 (num);
  tdenm = den[M+1];
  for (i=M; i<NI; i++)
    {
      /* Find trial quotient digit (the radix is 65536). */
      tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];

      /* Do not execute the divide instruction if it will overflow. */
      if ((tdenm * 0xffffL) < tnum)
	tquot = 0xffff;
      else
	tquot = tnum / tdenm;
      /* Multiply denominator by trial quotient digit. */
      m16m ((unsigned int)tquot, den, tprod);
      /* The quotient digit may have been overestimated. */
      if (ecmpm (tprod, num) > 0)
	{
	  tquot -= 1;
	  esubm (den, tprod);
	  if (ecmpm (tprod, num) > 0)
	    {
	      tquot -= 1;
	      esubm (den, tprod);
	    }
	}
      esubm (tprod, num);
      equot[i] = tquot;
      eshup6(num);
    }
  /* test for nonzero remainder after roundoff bit */
  p = &num[M];
  j = 0;
  for (i=M; i<NI; i++)
    {
      j |= *p++;
    }
  if (j)
    j = 1;

  for (i=0; i<NI; i++)
    num[i] = equot[i];

  return ((int)j);
}

/* Multiply significands of exploded e-type A and B, result in B.  */

static int
emulm (a, b)
     unsigned EMUSHORT a[], b[];
{
  unsigned EMUSHORT *p, *q;
  unsigned EMUSHORT pprod[NI];
  unsigned EMUSHORT j;
  int i;

  equot[0] = b[0];
  equot[1] = b[1];
  for (i=M; i<NI; i++)
    equot[i] = 0;

  j = 0;
  p = &a[NI-1];
  q = &equot[NI-1];
  for (i=M+1; i<NI; i++)
    {
      if (*p == 0)
	{
	  --p;
	}
      else
	{
	  m16m ((unsigned int) *p--, b, pprod);
	  eaddm(pprod, equot);
	}
      j |= *q;
      eshdn6(equot);
    }

  for (i=0; i<NI; i++)
    b[i] = equot[i];

  /* return flag for lost nonzero bits */
  return ((int)j);
}
#endif


/* Normalize and round off.

  The internal format number to be rounded is S.
  Input LOST is 0 if the value is exact.  This is the so-called sticky bit.
 
  Input SUBFLG indicates whether the number was obtained
  by a subtraction operation.  In that case if LOST is nonzero
  then the number is slightly smaller than indicated.
 
  Input EXP is the biased exponent, which may be negative.
  the exponent field of S is ignored but is replaced by
  EXP as adjusted by normalization and rounding.
 
  Input RCNTRL is the rounding control.  If it is nonzero, the
  returned value will be rounded to RNDPRC bits.

  For future reference:  In order for emdnorm to round off denormal
   significands at the right point, the input exponent must be
   adjusted to be the actual value it would have after conversion to
   the final floating point type.  This adjustment has been
   implemented for all type conversions (etoe53, etc.) and decimal
   conversions, but not for the arithmetic functions (eadd, etc.). 
   Data types having standard 15-bit exponents are not affected by
   this, but SFmode and DFmode are affected. For example, ediv with
   rndprc = 24 will not round correctly to 24-bit precision if the
   result is denormal.   */

static int rlast = -1;
static int rw = 0;
static unsigned EMUSHORT rmsk = 0;
static unsigned EMUSHORT rmbit = 0;
static unsigned EMUSHORT rebit = 0;
static int re = 0;
static unsigned EMUSHORT rbit[NI];

static void 
emdnorm (s, lost, subflg, exp, rcntrl)
     unsigned EMUSHORT s[];
     int lost;
     int subflg;
     EMULONG exp;
     int rcntrl;
{
  int i, j;
  unsigned EMUSHORT r;

  /* Normalize */
  j = enormlz (s);

  /* a blank significand could mean either zero or infinity. */
#ifndef INFINITY
  if (j > NBITS)
    {
      ecleazs (s);
      return;
    }
#endif
  exp -= j;
#ifndef INFINITY
  if (exp >= 32767L)
    goto overf;
#else
  if ((j > NBITS) && (exp < 32767))
    {
      ecleazs (s);
      return;
    }
#endif
  if (exp < 0L)
    {
      if (exp > (EMULONG) (-NBITS - 1))
	{
	  j = (int) exp;
	  i = eshift (s, j);
	  if (i)
	    lost = 1;
	}
      else
	{
	  ecleazs (s);
	  return;
	}
    }
  /* Round off, unless told not to by rcntrl. */
  if (rcntrl == 0)
    goto mdfin;
  /* Set up rounding parameters if the control register changed. */
  if (rndprc != rlast)
    {
      ecleaz (rbit);
      switch (rndprc)
	{
	default:
	case NBITS:
	  rw = NI - 1;		/* low guard word */
	  rmsk = 0xffff;
	  rmbit = 0x8000;
	  re = rw - 1;
	  rebit = 1;
	  break;