real.c

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

 	emov (a, b)		b = a
 	emul (a, b, c)		c = b * a
 	eneg (e)			e = -e
 	eround (a, b)		b = nearest integer value to a
 	esub (a, b, c)		c = b - a
 	e24toasc (&f, str, n)	single to ASCII string, n digits after decimal
 	e53toasc (&d, str, n)	double to ASCII string, n digits after decimal
 	e64toasc (&d, str, n)	80-bit long double to ASCII string
 	e113toasc (&d, str, n)	128-bit long double to ASCII string
 	etoasc (e, str, n)	e to ASCII string, n digits after decimal
 	etoe24 (e, &f)		convert e type to IEEE single precision
 	etoe53 (e, &d)		convert e type to IEEE double precision
 	etoe64 (e, &d)		convert e type to IEEE long double precision
 	ltoe (&l, e)		HOST_WIDE_INT to e type
 	ultoe (&l, e)		unsigned HOST_WIDE_INT to e type
	eisneg (e)              1 if sign bit of e != 0, else 0
	eisinf (e)              1 if e has maximum exponent (non-IEEE)
 				or is infinite (IEEE)
        eisnan (e)              1 if e is a NaN
 

 		Routines for internal format exploded e-type numbers
 
 	eaddm (ai, bi)		add significands, bi = bi + ai
 	ecleaz (ei)		ei = 0
 	ecleazs (ei)		set ei = 0 but leave its sign alone
 	ecmpm (ai, bi)		compare significands, return 1, 0, or -1
 	edivm (ai, bi)		divide  significands, bi = bi / ai
 	emdnorm (ai,l,s,exp)	normalize and round off
 	emovi (a, ai)		convert external a to internal ai
 	emovo (ai, a)		convert internal ai to external a
 	emovz (ai, bi)		bi = ai, low guard word of bi = 0
 	emulm (ai, bi)		multiply significands, bi = bi * ai
 	enormlz (ei)		left-justify the significand
 	eshdn1 (ai)		shift significand and guards down 1 bit
 	eshdn8 (ai)		shift down 8 bits
 	eshdn6 (ai)		shift down 16 bits
 	eshift (ai, n)		shift ai n bits up (or down if n < 0)
 	eshup1 (ai)		shift significand and guards up 1 bit
 	eshup8 (ai)		shift up 8 bits
 	eshup6 (ai)		shift up 16 bits
 	esubm (ai, bi)		subtract significands, bi = bi - ai
        eiisinf (ai)            1 if infinite
        eiisnan (ai)            1 if a NaN
 	eiisneg (ai)		1 if sign bit of ai != 0, else 0
        einan (ai)              set ai = NaN
        eiinfin (ai)            set ai = infinity

  The result is always normalized and rounded to NI-4 word precision
  after each arithmetic operation.

  Exception flags are NOT fully supported.
 
  Signaling NaN's are NOT supported; they are treated the same
  as quiet NaN's.
 
  Define INFINITY for support of infinity; otherwise a
  saturation arithmetic is implemented.
 
  Define NANS for support of Not-a-Number items; otherwise the
  arithmetic will never produce a NaN output, and might be confused
  by a NaN input.
  If NaN's are supported, the output of `ecmp (a,b)' is -2 if
  either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
  may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
  if in doubt.
 
  Denormals are always supported here where appropriate (e.g., not
  for conversion to DEC numbers).  */

/* Definitions for error codes that are passed to the common error handling
   routine mtherr.

   For Digital Equipment PDP-11 and VAX computers, certain
  IBM systems, and others that use numbers with a 56-bit
  significand, the symbol DEC should be defined.  In this
  mode, most floating point constants are given as arrays
  of octal integers to eliminate decimal to binary conversion
  errors that might be introduced by the compiler.
 
  For computers, such as IBM PC, that follow the IEEE
  Standard for Binary Floating Point Arithmetic (ANSI/IEEE
  Std 754-1985), the symbol IEEE should be defined.
  These numbers have 53-bit significands.  In this mode, constants
  are provided as arrays of hexadecimal 16 bit integers.
  The endian-ness of generated values is controlled by
  REAL_WORDS_BIG_ENDIAN.
 
  To accommodate other types of computer arithmetic, all
  constants are also provided in a normal decimal radix
  which one can hope are correctly converted to a suitable
  format by the available C language compiler.  To invoke
  this mode, the symbol UNK is defined.
 
  An important difference among these modes is a predefined
  set of machine arithmetic constants for each.  The numbers
  MACHEP (the machine roundoff error), MAXNUM (largest number
  represented), and several other parameters are preset by
  the configuration symbol.  Check the file const.c to
  ensure that these values are correct for your computer.
 
  For ANSI C compatibility, define ANSIC equal to 1.  Currently
  this affects only the atan2 function and others that use it. */

/* Constant definitions for math error conditions.  */

#define DOMAIN		1	/* argument domain error */
#define SING		2	/* argument singularity */
#define OVERFLOW	3	/* overflow range error */
#define UNDERFLOW	4	/* underflow range error */
#define TLOSS		5	/* total loss of precision */
#define PLOSS		6	/* partial loss of precision */
#define INVALID		7	/* NaN-producing operation */

/*  e type constants used by high precision check routines */

#if LONG_DOUBLE_TYPE_SIZE == 128
/* 0.0 */
unsigned EMUSHORT ezero[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
extern unsigned EMUSHORT ezero[];

/* 5.0E-1 */
unsigned EMUSHORT ehalf[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
extern unsigned EMUSHORT ehalf[];

/* 1.0E0 */
unsigned EMUSHORT eone[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
extern unsigned EMUSHORT eone[];

/* 2.0E0 */
unsigned EMUSHORT etwo[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
extern unsigned EMUSHORT etwo[];

/* 3.2E1 */
unsigned EMUSHORT e32[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
extern unsigned EMUSHORT e32[];

/* 6.93147180559945309417232121458176568075500134360255E-1 */
unsigned EMUSHORT elog2[NE] =
 {0x40f3, 0xf6af, 0x03f2, 0xb398,
  0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
extern unsigned EMUSHORT elog2[];

/* 1.41421356237309504880168872420969807856967187537695E0 */
unsigned EMUSHORT esqrt2[NE] =
 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
  0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
extern unsigned EMUSHORT esqrt2[];

/* 3.14159265358979323846264338327950288419716939937511E0 */
unsigned EMUSHORT epi[NE] =
 {0x2902, 0x1cd1, 0x80dc, 0x628b,
  0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
extern unsigned EMUSHORT epi[];

#else
/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
unsigned EMUSHORT ezero[NE] =
 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
unsigned EMUSHORT ehalf[NE] =
 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
unsigned EMUSHORT eone[NE] =
 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
unsigned EMUSHORT etwo[NE] =
 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
unsigned EMUSHORT e32[NE] =
 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
unsigned EMUSHORT elog2[NE] =
 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
unsigned EMUSHORT esqrt2[NE] =
 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
unsigned EMUSHORT epi[NE] =
 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
#endif

/* Control register for rounding precision.
   This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits.  */

int rndprc = NBITS;
extern int rndprc;

/*  Clear out entire e-type number X.  */

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

  for (i = 0; i < NE; i++)
    *x++ = 0;
}

/* Move e-type number from A to B.  */

static void 
emov (a, b)
     register unsigned EMUSHORT *a, *b;
{
  register int i;

  for (i = 0; i < NE; i++)
    *b++ = *a++;
}


/* Absolute value of e-type X.  */

static void 
eabs (x)
     unsigned EMUSHORT x[];
{
  /* sign is top bit of last word of external format */
  x[NE - 1] &= 0x7fff;		
}

/* Negate the e-type number X.  */

static void 
eneg (x)
     unsigned EMUSHORT x[];
{

  x[NE - 1] ^= 0x8000;		/* Toggle the sign bit */
}

/* Return 1 if sign bit of e-type number X is nonzero, else zero.  */

static int 
eisneg (x)
     unsigned EMUSHORT x[];
{

  if (x[NE - 1] & 0x8000)
    return (1);
  else
    return (0);
}

/* Return 1 if e-type number X is infinity, else return zero.  */

static int 
eisinf (x)
     unsigned EMUSHORT x[];
{

#ifdef NANS
  if (eisnan (x))
    return (0);
#endif
  if ((x[NE - 1] & 0x7fff) == 0x7fff)
    return (1);
  else
    return (0);
}

/* Check if e-type number is not a number.  The bit pattern is one that we
   defined, so we know for sure how to detect it.  */

static int 
eisnan (x)
     unsigned EMUSHORT x[];
{
#ifdef NANS
  int i;

  /* NaN has maximum exponent */
  if ((x[NE - 1] & 0x7fff) != 0x7fff)
    return (0);
  /* ... and non-zero significand field. */
  for (i = 0; i < NE - 1; i++)
    {
      if (*x++ != 0)
        return (1);
    }
#endif

  return (0);
}

/*  Fill e-type number X with infinity pattern (IEEE)
    or largest possible number (non-IEEE). */

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

#ifdef INFINITY
  for (i = 0; i < NE - 1; i++)
    *x++ = 0;
  *x |= 32767;
#else
  for (i = 0; i < NE - 1; i++)
    *x++ = 0xffff;
  *x |= 32766;
  if (rndprc < NBITS)
    {
      if (rndprc == 113)
	{
	  *(x - 9) = 0;
	  *(x - 8) = 0;
	}
      if (rndprc == 64)
	{
	  *(x - 5) = 0;
	}
      if (rndprc == 53)
	{
	  *(x - 4) = 0xf800;
	}
      else
	{
	  *(x - 4) = 0;
	  *(x - 3) = 0;
	  *(x - 2) = 0xff00;
	}
    }
#endif
}

/* Output an e-type NaN.
   This generates Intel's quiet NaN pattern for extended real.
   The exponent is 7fff, the leading mantissa word is c000.  */

static void 
enan (x, sign)
     register unsigned EMUSHORT *x;
     int sign;
{
  register int i;

  for (i = 0; i < NE - 2; i++)
    *x++ = 0;
  *x++ = 0xc000;
  *x = (sign << 15) | 0x7fff;
}

/* Move in an e-type number A, converting it to exploded e-type B.  */

static void 
emovi (a, b)
     unsigned EMUSHORT *a, *b;
{
  register unsigned EMUSHORT *p, *q;
  int i;

  q = b;
  p = a + (NE - 1);		/* point to last word of external number */
  /* get the sign bit */
  if (*p & 0x8000)
    *q++ = 0xffff;
  else
    *q++ = 0;
  /* get the exponent */
  *q = *p--;
  *q++ &= 0x7fff;		/* delete the sign bit */
#ifdef INFINITY
  if ((*(q - 1) & 0x7fff) == 0x7fff)
    {
#ifdef NANS
      if (eisnan (a))
	{
	  *q++ = 0;
	  for (i = 3; i < NI; i++)
	    *q++ = *p--;
	  return;
	}
#endif

      for (i = 2; i < NI; i++)
	*q++ = 0;
      return;
    }
#endif

  /* clear high guard word */
  *q++ = 0;
  /* move in the significand */
  for (i = 0; i < NE - 1; i++)
    *q++ = *p--;
  /* clear low guard word */
  *q = 0;
}

/* Move out exploded e-type number A, converting it to e type B.  */

static void 
emovo (a, b)
     unsigned EMUSHORT *a, *b;
{
  register unsigned EMUSHORT *p, *q;
  unsigned EMUSHORT i;
  int j;

  p = a;
  q = b + (NE - 1);		/* point to output exponent */
  /* combine sign and exponent */
  i = *p++;
  if (i)
    *q-- = *p++ | 0x8000;
  else
    *q-- = *p++;
#ifdef INFINITY
  if (*(p - 1) == 0x7fff)
    {
#ifdef NANS
      if (eiisnan (a))
	{
	  enan (b, eiisneg (a));
	  return;
	}
#endif
      einfin (b);
	return;
    }
#endif
  /* skip over guard word */
  ++p;
  /* move the significand */
  for (j = 0; j < NE - 1; j++)
    *q-- = *p++;
}

/* Clear out exploded e-type number XI.  */

static void 
ecleaz (xi)
     register unsigned EMUSHORT *xi;
{
  register int i;

  for (i = 0; i < NI; i++)
    *xi++ = 0;
}

/* Clear out exploded e-type XI, but don't touch the sign. */

static void 
ecleazs (xi)
     register unsigned EMUSHORT *xi;
{
  register int i;

  ++xi;
  for (i = 0; i < NI - 1; i++)
    *xi++ = 0;
}

/* Move exploded e-type number from A to B.  */

static void 
emovz (a, b)
     register unsigned EMUSHORT *a, *b;
{
  register int i;

  for (i = 0; i < NI - 1; i++)
    *b++ = *a++;
  /* clear low guard word */
  *b = 0;
}

/* Generate exploded e-type NaN.
   The explicit pattern for this is maximum exponent and
   top two significant bits set.  */

static void
einan (x)
     unsigned EMUSHORT x[];
{

  ecleaz (x);
  x[E] = 0x7fff;
  x[M + 1] = 0xc000;
}

/* Return nonzero if exploded e-type X is a NaN. */

static int 
eiisnan (x)
     unsigned EMUSHORT x[];
{
  int i;

  if ((x[E] & 0x7fff) == 0x7fff)
    {
      for (i = M + 1; i < NI; i++)
	{
	  if (x[i] != 0)
	    return (1);
	}
    }
  return (0);
}

/* Return nonzero if sign of exploded e-type X is nonzero.  */

static int 
eiisneg (x)
     unsigned EMUSHORT x[];
{

  return x[0] != 0;
}

/* Fill exploded e-type X with infinity pattern.
   This has maximum exponent and significand all zeros.  */

static void
eiinfin (x)
     unsigned EMUSHORT x[];
{

  ecleaz (x);
  x[E] = 0x7fff;
}

/* Return nonzero if exploded e-type X is infinite. */

static int 
eiisinf (x)
     unsigned EMUSHORT x[];
{

#ifdef NANS
  if (eiisnan (x))
    return (0);
#endif
  if ((x[E] & 0x7fff) == 0x7fff)
    return (1);
  return (0);
}


/* Compare significands of numbers in internal exploded e-type format.
   Guard words are included in the comparison.

   Returns	+1 if a > b
		 0 if a == b
		-1 if a < b   */

static int
ecmpm (a, b)
     register unsigned EMUSHORT *a, *b;
{
  int i;

  a += M;			/* skip up to significand area */
  b += M;
  for (i = M; i < NI; i++)
    {
      if (*a++ != *b++)
	goto difrnt;
    }
  return (0);

 difrnt:
  if (*(--a) > *(--b))
    return (1);
  else
    return (-1);
}

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

static void 
eshdn1 (x)
     register unsigned EMUSHORT *x;
{
  register unsigned EMUSHORT bits;