real.c

来自「GCC编译器源代码」· C语言 代码 · 共 2,875 行 · 第 1/5 页

      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;
	case 113:
	  rw = 10;
	  rmsk = 0x7fff;
	  rmbit = 0x4000;
	  rebit = 0x8000;
	  re = rw;
	  break;
	case 64:
	  rw = 7;
	  rmsk = 0xffff;
	  rmbit = 0x8000;
	  re = rw - 1;
	  rebit = 1;
	  break;
	  /* For DEC or IBM arithmetic */
	case 56:
	  rw = 6;
	  rmsk = 0xff;
	  rmbit = 0x80;
	  rebit = 0x100;
	  re = rw;
	  break;
	case 53:
	  rw = 6;
	  rmsk = 0x7ff;
	  rmbit = 0x0400;
	  rebit = 0x800;
	  re = rw;
	  break;
	case 24:
	  rw = 4;
	  rmsk = 0xff;
	  rmbit = 0x80;
	  rebit = 0x100;
	  re = rw;
	  break;
	}
      rbit[re] = rebit;
      rlast = rndprc;
    }

  /* Shift down 1 temporarily if the data structure has an implied
     most significant bit and the number is denormal.
     Intel long double denormals also lose one bit of precision.  */
  if ((exp <= 0) && (rndprc != NBITS)
      && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
    {
      lost |= s[NI - 1] & 1;
      eshdn1 (s);
    }
  /* Clear out all bits below the rounding bit,
     remembering in r if any were nonzero.  */
  r = s[rw] & rmsk;
  if (rndprc < NBITS)
    {
      i = rw + 1;
      while (i < NI)
	{
	  if (s[i])
	    r |= 1;
	  s[i] = 0;
	  ++i;
	}
    }
  s[rw] &= ~rmsk;
  if ((r & rmbit) != 0)
    {
      if (r == rmbit)
	{
	  if (lost == 0)
	    {			/* round to even */
	      if ((s[re] & rebit) == 0)
		goto mddone;
	    }
	  else
	    {
	      if (subflg != 0)
		goto mddone;
	    }
	}
      eaddm (rbit, s);
    }
 mddone:
/* Undo the temporary shift for denormal values.  */
  if ((exp <= 0) && (rndprc != NBITS)
      && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
    {
      eshup1 (s);
    }
  if (s[2] != 0)
    {				/* overflow on roundoff */
      eshdn1 (s);
      exp += 1;
    }
 mdfin:
  s[NI - 1] = 0;
  if (exp >= 32767L)
    {
#ifndef INFINITY
    overf:
#endif
#ifdef INFINITY
      s[1] = 32767;
      for (i = 2; i < NI - 1; i++)
	s[i] = 0;
      if (extra_warnings)
	warning ("floating point overflow");
#else
      s[1] = 32766;
      s[2] = 0;
      for (i = M + 1; i < NI - 1; i++)
	s[i] = 0xffff;
      s[NI - 1] = 0;
      if ((rndprc < 64) || (rndprc == 113))
	{
	  s[rw] &= ~rmsk;
	  if (rndprc == 24)
	    {
	      s[5] = 0;
	      s[6] = 0;
	    }
	}
#endif
      return;
    }
  if (exp < 0)
    s[1] = 0;
  else
    s[1] = (unsigned EMUSHORT) exp;
}

/*  Subtract.  C = B - A, all e type numbers.  */

static int subflg = 0;

static void 
esub (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{

#ifdef NANS
  if (eisnan (a))
    {
      emov (a, c);
      return;
    }
  if (eisnan (b))
    {
      emov (b, c);
      return;
    }
/* Infinity minus infinity is a NaN.
   Test for subtracting infinities of the same sign.  */
  if (eisinf (a) && eisinf (b)
      && ((eisneg (a) ^ eisneg (b)) == 0))
    {
      mtherr ("esub", INVALID);
      enan (c, 0);
      return;
    }
#endif
  subflg = 1;
  eadd1 (a, b, c);
}

/* Add.  C = A + B, all e type.  */

static void 
eadd (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{

#ifdef NANS
/* NaN plus anything is a NaN.  */
  if (eisnan (a))
    {
      emov (a, c);
      return;
    }
  if (eisnan (b))
    {
      emov (b, c);
      return;
    }
/* Infinity minus infinity is a NaN.
   Test for adding infinities of opposite signs.  */
  if (eisinf (a) && eisinf (b)
      && ((eisneg (a) ^ eisneg (b)) != 0))
    {
      mtherr ("esub", INVALID);
      enan (c, 0);
      return;
    }
#endif
  subflg = 0;
  eadd1 (a, b, c);
}

/* Arithmetic common to both addition and subtraction.  */

static void 
eadd1 (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
  int i, lost, j, k;
  EMULONG lt, lta, ltb;

#ifdef INFINITY
  if (eisinf (a))
    {
      emov (a, c);
      if (subflg)
	eneg (c);
      return;
    }
  if (eisinf (b))
    {
      emov (b, c);
      return;
    }
#endif
  emovi (a, ai);
  emovi (b, bi);
  if (subflg)
    ai[0] = ~ai[0];

  /* compare exponents */
  lta = ai[E];
  ltb = bi[E];
  lt = lta - ltb;
  if (lt > 0L)
    {				/* put the larger number in bi */
      emovz (bi, ci);
      emovz (ai, bi);
      emovz (ci, ai);
      ltb = bi[E];
      lt = -lt;
    }
  lost = 0;
  if (lt != 0L)
    {
      if (lt < (EMULONG) (-NBITS - 1))
	goto done;		/* answer same as larger addend */
      k = (int) lt;
      lost = eshift (ai, k);	/* shift the smaller number down */
    }
  else
    {
      /* exponents were the same, so must compare significands */
      i = ecmpm (ai, bi);
      if (i == 0)
	{			/* the numbers are identical in magnitude */
	  /* if different signs, result is zero */
	  if (ai[0] != bi[0])
	    {
	      eclear (c);
	      return;
	    }
	  /* if same sign, result is double */
	  /* double denormalized tiny number */
	  if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
	    {
	      eshup1 (bi);
	      goto done;
	    }
	  /* add 1 to exponent unless both are zero! */
	  for (j = 1; j < NI - 1; j++)
	    {
	      if (bi[j] != 0)
		{
		  ltb += 1;
		  if (ltb >= 0x7fff)
		    {
		      eclear (c);
		      if (ai[0] != 0)
			eneg (c);
		      einfin (c);
		      return;
		    }
		  break;
		}
	    }
	  bi[E] = (unsigned EMUSHORT) ltb;
	  goto done;
	}
      if (i > 0)
	{			/* put the larger number in bi */
	  emovz (bi, ci);
	  emovz (ai, bi);
	  emovz (ci, ai);
	}
    }
  if (ai[0] == bi[0])
    {
      eaddm (ai, bi);
      subflg = 0;
    }
  else
    {
      esubm (ai, bi);
      subflg = 1;
    }
  emdnorm (bi, lost, subflg, ltb, 64);

 done:
  emovo (bi, c);
}

/* Divide: C = B/A, all e type.  */

static void 
ediv (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  int i, sign;
  EMULONG lt, lta, ltb;

/* IEEE says if result is not a NaN, the sign is "-" if and only if
   operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
  sign = eisneg(a) ^ eisneg(b);

#ifdef NANS
/* Return any NaN input.  */
  if (eisnan (a))
    {
    emov (a, c);
    return;
    }
  if (eisnan (b))
    {
    emov (b, c);
    return;
    }
/* Zero over zero, or infinity over infinity, is a NaN.  */
  if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
      || (eisinf (a) && eisinf (b)))
    {
    mtherr ("ediv", INVALID);
    enan (c, sign);
    return;
    }
#endif
/* Infinity over anything else is infinity.  */
#ifdef INFINITY
  if (eisinf (b))
    {
      einfin (c);
      goto divsign;
    }
/* Anything else over infinity is zero.  */
  if (eisinf (a))
    {
      eclear (c);
      goto divsign;
    }
#endif
  emovi (a, ai);
  emovi (b, bi);
  lta = ai[E];
  ltb = bi[E];
  if (bi[E] == 0)
    {				/* See if numerator is zero.  */
      for (i = 1; i < NI - 1; i++)
	{
	  if (bi[i] != 0)
	    {
	      ltb -= enormlz (bi);
	      goto dnzro1;
	    }
	}
      eclear (c);
      goto divsign;
    }
 dnzro1:

  if (ai[E] == 0)
    {				/* possible divide by zero */
      for (i = 1; i < NI - 1; i++)
	{
	  if (ai[i] != 0)
	    {
	      lta -= enormlz (ai);
	      goto dnzro2;
	    }
	}
/* Divide by zero is not an invalid operation.
   It is a divide-by-zero operation!   */
      einfin (c);
      mtherr ("ediv", SING);
      goto divsign;
    }
 dnzro2:

  i = edivm (ai, bi);
  /* calculate exponent */
  lt = ltb - lta + EXONE;
  emdnorm (bi, i, 0, lt, 64);
  emovo (bi, c);

 divsign:

  if (sign
#ifndef IEEE
      && (ecmp (c, ezero) != 0)
#endif
      )
     *(c+(NE-1)) |= 0x8000;
  else
     *(c+(NE-1)) &= ~0x8000;
}

/* Multiply