real.c

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

  endian (e, l, TFmode);
}

/* Convert R to a double extended precision value.  The output array L
   contains three 32-bit pieces of the result, in the order they would
   appear in memory.  */

void 
etarldouble (r, l)
     REAL_VALUE_TYPE r;
     long l[];
{
  unsigned EMUSHORT e[NE];

  GET_REAL (&r, e);
  etoe64 (e, e);
  endian (e, l, XFmode);
}

/* Convert R to a double precision value.  The output array L contains two
   32-bit pieces of the result, in the order they would appear in memory.  */

void 
etardouble (r, l)
     REAL_VALUE_TYPE r;
     long l[];
{
  unsigned EMUSHORT e[NE];

  GET_REAL (&r, e);
  etoe53 (e, e);
  endian (e, l, DFmode);
}

/* Convert R to a single precision float value stored in the least-significant
   bits of a `long'.  */

long
etarsingle (r)
     REAL_VALUE_TYPE r;
{
  unsigned EMUSHORT e[NE];
  long l;

  GET_REAL (&r, e);
  etoe24 (e, e);
  endian (e, &l, SFmode);
  return ((long) l);
}

/* Convert X to a decimal ASCII string S for output to an assembly
   language file.  Note, there is no standard way to spell infinity or
   a NaN, so these values may require special treatment in the tm.h
   macros.  */

void
ereal_to_decimal (x, s)
     REAL_VALUE_TYPE x;
     char *s;
{
  unsigned EMUSHORT e[NE];

  GET_REAL (&x, e);
  etoasc (e, s, 20);
}

/* Compare X and Y.  Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
   or -2 if either is a NaN.   */

int
ereal_cmp (x, y)
     REAL_VALUE_TYPE x, y;
{
  unsigned EMUSHORT ex[NE], ey[NE];

  GET_REAL (&x, ex);
  GET_REAL (&y, ey);
  return (ecmp (ex, ey));
}

/*  Return 1 if the sign bit of X is set, else return 0.  */

int
ereal_isneg (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT ex[NE];

  GET_REAL (&x, ex);
  return (eisneg (ex));
}

/* End of REAL_ARITHMETIC interface */

/*
  Extended precision IEEE binary floating point arithmetic routines

  Numbers are stored in C language as arrays of 16-bit unsigned
  short integers.  The arguments of the routines are pointers to
  the arrays.

  External e type data structure, similar to Intel 8087 chip
  temporary real format but possibly with a larger significand:

	NE-1 significand words	(least significant word first,
				 most significant bit is normally set)
	exponent		(value = EXONE for 1.0,
				top bit is the sign)


  Internal exploded e-type data structure of a number (a "word" is 16 bits):

  ei[0]	sign word	(0 for positive, 0xffff for negative)
  ei[1]	biased exponent	(value = EXONE for the number 1.0)
  ei[2]	high guard word	(always zero after normalization)
  ei[3]
  to ei[NI-2]	significand	(NI-4 significand words,
 				 most significant word first,
 				 most significant bit is set)
  ei[NI-1]	low guard word	(0x8000 bit is rounding place)
 
 
 
 		Routines for external format e-type numbers
 
 	asctoe (string, e)	ASCII string to extended double e type
 	asctoe64 (string, &d)	ASCII string to long double
 	asctoe53 (string, &d)	ASCII string to double
 	asctoe24 (string, &f)	ASCII string to single
 	asctoeg (string, e, prec) ASCII string to specified precision
 	e24toe (&f, e)		IEEE single precision to e type
 	e53toe (&d, e)		IEEE double precision to e type
 	e64toe (&d, e)		IEEE long double precision to e type
 	e113toe (&d, e)		128-bit long double precision to e type
 	eabs (e)			absolute value
 	eadd (a, b, c)		c = b + a
 	eclear (e)		e = 0
 	ecmp (a, b)		Returns 1 if a > b, 0 if a == b,
 				-1 if a < b, -2 if either a or b is a NaN.
 	ediv (a, b, c)		c = b / a
 	efloor (a, b)		truncate to integer, toward -infinity
 	efrexp (a, exp, s)	extract exponent and significand
 	eifrac (e, &l, frac)    e to HOST_WIDE_INT and e type fraction
 	euifrac (e, &l, frac)   e to unsigned HOST_WIDE_INT and e type fraction
 	einfin (e)		set e to infinity, leaving its sign alone
 	eldexp (a, n, b)	multiply by 2**n
 	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 */