real.c

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

	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)
		{
		  /* This could overflow, but let emovo take care of that. */
		  ltb += 1;
		  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 e-types A and B, return e-type product C.   */

static void 
emul (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  int i, j, 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
/* NaN times anything is the same NaN. */
  if (eisnan (a))
    {
    emov (a, c);
    return;
    }
  if (eisnan (b))
    {
    emov (b, c);
    return;
    }
/* Zero times infinity is a NaN. */
  if ((eisinf (a) && (ecmp (b, ezero) == 0))
      || (eisinf (b) && (ecmp (a, ezero) == 0)))
    {
    mtherr ("emul", INVALID);
    enan (c, sign);
    return;
    }
#endif
/* Infinity times anything else is infinity. */
#ifdef INFINITY
  if (eisinf (a) || eisinf (b))
    {
      einfin (c);
      goto mulsign;
    }
#endif
  emovi (a, ai);
  emovi (b, bi);
  lta = ai[E];
  ltb = bi[E];
  if (ai[E] == 0)
    {
      for (i = 1; i < NI - 1; i++)
	{
	  if (ai[i] != 0)
	    {
	      lta -= enormlz (ai);
	      goto mnzer1;
	    }
	}
      eclear (c);
      goto mulsign;
    }
 mnzer1:

  if (bi[E] == 0)
    {
      for (i = 1; i < NI - 1; i++)
	{
	  if (bi[i] != 0)
	    {
	      ltb -= enormlz (bi);
	      goto mnzer2;
	    }
	}
      eclear (c);
      goto mulsign;
    }
 mnzer2:

  /* Multiply significands */
  j = emulm (ai, bi);
  /* calculate exponent */
  lt = lta + ltb - (EXONE - 1);
  emdnorm (bi, j, 0, lt, 64);
  emovo (bi, c);

 mulsign:

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

/* Convert double precision PE to e-type Y.  */

static void
e53toe (pe, y)
     unsigned EMUSHORT *pe, *y;
{
#ifdef DEC

  dectoe (pe, y);

#else
#ifdef IBM

  ibmtoe (pe, y, DFmode);

#else
  register unsigned EMUSHORT r;
  register unsigned EMUSHORT *e, *p;
  unsigned EMUSHORT yy[NI];
  int denorm, k;

  e = pe;
  denorm = 0;			/* flag if denormalized number */
  ecleaz (yy);
  if (! REAL_WORDS_BIG_ENDIAN)
    e += 3;
  r = *e;
  yy[0] = 0;
  if (r & 0x8000)
    yy[0] = 0xffff;
  yy[M] = (r & 0x0f) | 0x10;
  r &= ~0x800f;			/* strip sign and 4 significand bits */
#ifdef INFINITY
  if (r == 0x7ff0)
    {
#ifdef NANS
      if (! REAL_WORDS_BIG_ENDIAN)
	{
	  if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
	      || (pe[1] != 0) || (pe[0] != 0))
	    {
	      enan (y, yy[0] != 0);
	      return;
	    }
	}
      else
	{
	  if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
	      || (pe[2] != 0) || (pe[3] != 0))
	    {
	      enan (y, yy[0] != 0);
	      return;
	    }
	}
#endif  /* NANS */
      eclear (y);
      einfin (y);
      if (yy[0])
	eneg (y);
      return;
    }
#endif  /* INFINITY */
  r >>= 4;
  /* If zero exponent, then the significand is denormalized.
     So take back the understood high significand bit. */

  if (r == 0)
    {
      denorm = 1;
      yy[M] &= ~0x10;
    }
  r += EXONE - 01777;
  yy[E] = r;
  p = &yy[M + 1];
#ifdef IEEE
  if (! REAL_WORDS_BIG_ENDIAN)
    {
      *p++ = *(--e);
      *p++ = *(