[thirdparty/glibc.git] / soft-fp / extended.h

/* Software floating-point emulation.
   Definitions for IEEE Extended Precision.
   Copyright (C) 1999 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Jakub Jelinek (jj@ultra.linux.cz).

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307 USA.  */

#if _FP_W_TYPE_SIZE < 32
#error "Here's a nickel, kid. Go buy yourself a real computer."
#endif

#if _FP_W_TYPE_SIZE < 64
#define _FP_FRACTBITS_E         (4*_FP_W_TYPE_SIZE)
#else
#define _FP_FRACTBITS_E		(2*_FP_W_TYPE_SIZE)
#endif

#define _FP_FRACBITS_E		64
#define _FP_FRACXBITS_E		(_FP_FRACTBITS_E - _FP_FRACBITS_E)
#define _FP_WFRACBITS_E		(_FP_WORKBITS + _FP_FRACBITS_E)
#define _FP_WFRACXBITS_E	(_FP_FRACTBITS_E - _FP_WFRACBITS_E)
#define _FP_EXPBITS_E		15
#define _FP_EXPBIAS_E		16383
#define _FP_EXPMAX_E		32767

#define _FP_QNANBIT_E		\
	((_FP_W_TYPE)1 << (_FP_FRACBITS_E-2) % _FP_W_TYPE_SIZE)
#define _FP_IMPLBIT_E		\
	((_FP_W_TYPE)1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE)
#define _FP_OVERFLOW_E		\
	((_FP_W_TYPE)1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE))

#if _FP_W_TYPE_SIZE < 64

union _FP_UNION_E
{
   long double flt;
   struct 
   {
#if __BYTE_ORDER == __BIG_ENDIAN
      unsigned long pad1 : _FP_W_TYPE_SIZE;
      unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
      unsigned long sign : 1;
      unsigned long exp : _FP_EXPBITS_E;
      unsigned long frac1 : _FP_W_TYPE_SIZE;
      unsigned long frac0 : _FP_W_TYPE_SIZE;
#else
      unsigned long frac0 : _FP_W_TYPE_SIZE;
      unsigned long frac1 : _FP_W_TYPE_SIZE;
      unsigned exp : _FP_EXPBITS_E;
      unsigned sign : 1;
#endif /* not bigendian */
   } bits __attribute__((packed));
};


#define FP_DECL_E(X)		_FP_DECL(4,X)

#define FP_UNPACK_RAW_E(X, val)				\
  do {							\
    union _FP_UNION_E _flo; _flo.flt = (val);		\
							\
    X##_f[2] = 0; X##_f[3] = 0;				\
    X##_f[0] = _flo.bits.frac0;				\
    X##_f[1] = _flo.bits.frac1;				\
    X##_e  = _flo.bits.exp;				\
    X##_s  = _flo.bits.sign;				\
    if (!X##_e && (X##_f[1] || X##_f[0])		\
        && !(X##_f[1] & _FP_IMPLBIT_E))			\
      {							\
        X##_e++;					\
        FP_SET_EXCEPTION(FP_EX_DENORM);			\
      }							\
  } while (0)

#define FP_UNPACK_RAW_EP(X, val)			\
  do {							\
    union _FP_UNION_E *_flo =				\
    (union _FP_UNION_E *)(val);				\
							\
    X##_f[2] = 0; X##_f[3] = 0;				\
    X##_f[0] = _flo->bits.frac0;			\
    X##_f[1] = _flo->bits.frac1;			\
    X##_e  = _flo->bits.exp;				\
    X##_s  = _flo->bits.sign;				\
    if (!X##_e && (X##_f[1] || X##_f[0])		\
        && !(X##_f[1] & _FP_IMPLBIT_E))			\
      {							\
        X##_e++;					\
        FP_SET_EXCEPTION(FP_EX_DENORM);			\
      }							\
  } while (0)

#define FP_PACK_RAW_E(val, X)				\
  do {							\
    union _FP_UNION_E _flo;				\
							\
    if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;		\
    else X##_f[1] &= ~(_FP_IMPLBIT_E);			\
    _flo.bits.frac0 = X##_f[0];				\
    _flo.bits.frac1 = X##_f[1];				\
    _flo.bits.exp   = X##_e;				\
    _flo.bits.sign  = X##_s;				\
							\
    (val) = _flo.flt;					\
  } while (0)

#define FP_PACK_RAW_EP(val, X)				\
  do {							\
    if (!FP_INHIBIT_RESULTS)				\
      {							\
	union _FP_UNION_E *_flo =			\
	  (union _FP_UNION_E *)(val);			\
							\
	if (X##_e) X##_f[1] |= _FP_IMPLBIT_E;		\
	else X##_f[1] &= ~(_FP_IMPLBIT_E);		\
	_flo->bits.frac0 = X##_f[0];			\
	_flo->bits.frac1 = X##_f[1];			\
	_flo->bits.exp   = X##_e;			\
	_flo->bits.sign  = X##_s;			\
      }							\
  } while (0)

#define FP_UNPACK_E(X,val)		\
  do {					\
    FP_UNPACK_RAW_E(X,val);		\
    _FP_UNPACK_CANONICAL(E,4,X);	\
  } while (0)

#define FP_UNPACK_EP(X,val)		\
  do {					\
    FP_UNPACK_RAW_2_P(X,val);		\
    _FP_UNPACK_CANONICAL(E,4,X);	\
  } while (0)

#define FP_PACK_E(val,X)		\
  do {					\
    _FP_PACK_CANONICAL(E,4,X);		\
    FP_PACK_RAW_E(val,X);		\
  } while (0)

#define FP_PACK_EP(val,X)		\
  do {					\
    _FP_PACK_CANONICAL(E,4,X);		\
    FP_PACK_RAW_EP(val,X);		\
  } while (0)

#define FP_ISSIGNAN_E(X)	_FP_ISSIGNAN(E,4,X)
#define FP_NEG_E(R,X)		_FP_NEG(E,4,R,X)
#define FP_ADD_E(R,X,Y)		_FP_ADD(E,4,R,X,Y)
#define FP_SUB_E(R,X,Y)		_FP_SUB(E,4,R,X,Y)
#define FP_MUL_E(R,X,Y)		_FP_MUL(E,4,R,X,Y)
#define FP_DIV_E(R,X,Y)		_FP_DIV(E,4,R,X,Y)
#define FP_SQRT_E(R,X)		_FP_SQRT(E,4,R,X)

/*
 * Square root algorithms:
 * We have just one right now, maybe Newton approximation
 * should be added for those machines where division is fast.
 * This has special _E version because standard _4 square
 * root would not work (it has to start normally with the
 * second word and not the first), but as we have to do it
 * anyway, we optimize it by doing most of the calculations
 * in two UWtype registers instead of four.
 */
 
#define _FP_SQRT_MEAT_E(R, S, T, X, q)			\
  do {							\
    q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);		\
    _FP_FRAC_SRL_4(X, (_FP_WORKBITS));			\
    while (q)						\
      {							\
	T##_f[1] = S##_f[1] + q;			\
	if (T##_f[1] <= X##_f[1])			\
	  {						\
	    S##_f[1] = T##_f[1] + q;			\
	    X##_f[1] -= T##_f[1];			\
	    R##_f[1] += q;				\
	  }						\
	_FP_FRAC_SLL_2(X, 1);				\
	q >>= 1;					\
      }							\
    q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);		\
    while (q)						\
      {							\
	T##_f[0] = S##_f[0] + q;			\
	T##_f[1] = S##_f[1];				\
	if (T##_f[1] < X##_f[1] || 			\
	    (T##_f[1] == X##_f[1] &&			\
	     T##_f[0] <= X##_f[0]))			\
	  {						\
	    S##_f[0] = T##_f[0] + q;			\
	    S##_f[1] += (T##_f[0] > S##_f[0]);		\
	    _FP_FRAC_DEC_2(X, T);			\
	    R##_f[0] += q;				\
	  }						\
	_FP_FRAC_SLL_2(X, 1);				\
	q >>= 1;					\
      }							\
    _FP_FRAC_SLL_4(R, (_FP_WORKBITS));			\
    if (X##_f[0] | X##_f[1])				\
      {							\
	if (S##_f[1] < X##_f[1] || 			\
	    (S##_f[1] == X##_f[1] &&			\
	     S##_f[0] < X##_f[0]))			\
	  R##_f[0] |= _FP_WORK_ROUND;			\
	R##_f[0] |= _FP_WORK_STICKY;			\
      }							\
  } while (0)

#define FP_CMP_E(r,X,Y,un)	_FP_CMP(E,4,r,X,Y,un)
#define FP_CMP_EQ_E(r,X,Y)	_FP_CMP_EQ(E,4,r,X,Y)
#define FP_CMP_UNORD_E(r,X,Y)	_FP_CMP_UNORD(E,4,r,X,Y)

#define FP_TO_INT_E(r,X,rsz,rsg)	_FP_TO_INT(E,4,r,X,rsz,rsg)
#define FP_FROM_INT_E(X,r,rs,rt)	_FP_FROM_INT(E,4,X,r,rs,rt)

#define _FP_FRAC_HIGH_E(X)	(X##_f[2])
#define _FP_FRAC_HIGH_RAW_E(X)	(X##_f[1])

#else   /* not _FP_W_TYPE_SIZE < 64 */
union _FP_UNION_E
{
  long double flt /* __attribute__((mode(TF))) */ ;
  struct {
#if __BYTE_ORDER == __BIG_ENDIAN
    unsigned long pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
    unsigned sign  : 1;
    unsigned exp   : _FP_EXPBITS_E;
    unsigned long frac : _FP_W_TYPE_SIZE;
#else
    unsigned long frac : _FP_W_TYPE_SIZE;
    unsigned exp   : _FP_EXPBITS_E;
    unsigned sign  : 1;
#endif
  } bits;
};

#define FP_DECL_E(X)		_FP_DECL(2,X)

#define FP_UNPACK_RAW_E(X, val)					\
  do {								\
    union _FP_UNION_E _flo; _flo.flt = (val);			\
								\
    X##_f0 = _flo.bits.frac;					\
    X##_f1 = 0;							\
    X##_e = _flo.bits.exp;					\
    X##_s = _flo.bits.sign;					\
    if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))		\
      {								\
        X##_e++;						\
        FP_SET_EXCEPTION(FP_EX_DENORM);				\
      }								\
  } while (0)

#define FP_UNPACK_RAW_EP(X, val)				\
  do {								\
    union _FP_UNION_E *_flo =					\
      (union _FP_UNION_E *)(val);				\
								\
    X##_f0 = _flo->bits.frac;					\
    X##_f1 = 0;							\
    X##_e = _flo->bits.exp;					\
    X##_s = _flo->bits.sign;					\
    if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E))		\
      {								\
        X##_e++;						\
        FP_SET_EXCEPTION(FP_EX_DENORM);				\
      }								\
  } while (0)

#define FP_PACK_RAW_E(val, X)					\
  do {								\
    union _FP_UNION_E _flo;					\
								\
    if (X##_e) X##_f0 |= _FP_IMPLBIT_E;				\
    else X##_f0 &= ~(_FP_IMPLBIT_E);				\
    _flo.bits.frac = X##_f0;					\
    _flo.bits.exp  = X##_e;					\
    _flo.bits.sign = X##_s;					\
								\
    (val) = _flo.flt;						\
  } while (0)

#define FP_PACK_RAW_EP(fs, val, X)				\
  do {								\
    if (!FP_INHIBIT_RESULTS)					\
      {								\
	union _FP_UNION_E *_flo =				\
	  (union _FP_UNION_E *)(val);				\
								\
	if (X##_e) X##_f0 |= _FP_IMPLBIT_E;			\
	else X##_f0 &= ~(_FP_IMPLBIT_E);			\
	_flo->bits.frac = X##_f0;				\
	_flo->bits.exp  = X##_e;				\
	_flo->bits.sign = X##_s;				\
      }								\
  } while (0)


#define FP_UNPACK_E(X,val)		\
  do {					\
    FP_UNPACK_RAW_E(X,val);		\
    _FP_UNPACK_CANONICAL(E,2,X);	\
  } while (0)

#define FP_UNPACK_EP(X,val)		\
  do {					\
    FP_UNPACK_RAW_EP(X,val);		\
    _FP_UNPACK_CANONICAL(E,2,X);	\
  } while (0)

#define FP_PACK_E(val,X)		\
  do {					\
    _FP_PACK_CANONICAL(E,2,X);		\
    FP_PACK_RAW_E(val,X);		\
  } while (0)

#define FP_PACK_EP(val,X)		\
  do {					\
    _FP_PACK_CANONICAL(E,2,X);		\
    FP_PACK_RAW_EP(val,X);		\
  } while (0)

#define FP_ISSIGNAN_E(X)	_FP_ISSIGNAN(E,2,X)
#define FP_NEG_E(R,X)		_FP_NEG(E,2,R,X)
#define FP_ADD_E(R,X,Y)		_FP_ADD(E,2,R,X,Y)
#define FP_SUB_E(R,X,Y)		_FP_SUB(E,2,R,X,Y)
#define FP_MUL_E(R,X,Y)		_FP_MUL(E,2,R,X,Y)
#define FP_DIV_E(R,X,Y)		_FP_DIV(E,2,R,X,Y)
#define FP_SQRT_E(R,X)		_FP_SQRT(E,2,R,X)

/*
 * Square root algorithms:
 * We have just one right now, maybe Newton approximation
 * should be added for those machines where division is fast.
 * We optimize it by doing most of the calculations
 * in one UWtype registers instead of two, although we don't
 * have to.
 */
#define _FP_SQRT_MEAT_E(R, S, T, X, q)			\
  do {							\
    q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);		\
    _FP_FRAC_SRL_2(X, (_FP_WORKBITS));			\
    while (q)						\
      {							\
        T##_f0 = S##_f0 + q;				\
        if (T##_f0 <= X##_f0)				\
          {						\
            S##_f0 = T##_f0 + q;			\
            X##_f0 -= T##_f0;				\
            R##_f0 += q;				\
          }						\
        _FP_FRAC_SLL_1(X, 1);				\
        q >>= 1;					\
      }							\
    _FP_FRAC_SLL_2(R, (_FP_WORKBITS));			\
    if (X##_f0)						\
      {							\
	if (S##_f0 < X##_f0)				\
	  R##_f0 |= _FP_WORK_ROUND;			\
	R##_f0 |= _FP_WORK_STICKY;			\
      }							\
  } while (0)
 
#define FP_CMP_E(r,X,Y,un)	_FP_CMP(E,2,r,X,Y,un)
#define FP_CMP_EQ_E(r,X,Y)	_FP_CMP_EQ(E,2,r,X,Y)
#define FP_CMP_UNORD_E(r,X,Y)	_FP_CMP_UNORD(E,2,r,X,Y)

#define FP_TO_INT_E(r,X,rsz,rsg)	_FP_TO_INT(E,2,r,X,rsz,rsg)
#define FP_FROM_INT_E(X,r,rs,rt)	_FP_FROM_INT(E,2,X,r,rs,rt)

#define _FP_FRAC_HIGH_E(X)	(X##_f1)
#define _FP_FRAC_HIGH_RAW_E(X)	(X##_f0)

#endif /* not _FP_W_TYPE_SIZE < 64 */
Commit	Line	Data
d876f532 UD	1	/* Software floating-point emulation.
d876f532 UD	2	Definitions for IEEE Extended Precision.
a334319f	3	Copyright (C) 1999 Free Software Foundation, Inc.
d876f532 UD	4	This file is part of the GNU C Library.
	5	Contributed by Jakub Jelinek (jj@ultra.linux.cz).
	6
	7	The GNU C Library is free software; you can redistribute it and/or
41bdb6e2 AJ	8	modify it under the terms of the GNU Lesser General Public
	9	License as published by the Free Software Foundation; either
	10	version 2.1 of the License, or (at your option) any later version.
d876f532 UD	11
	12	The GNU C Library is distributed in the hope that it will be useful,
	13	but WITHOUT ANY WARRANTY; without even the implied warranty of
	14	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
41bdb6e2	15	Lesser General Public License for more details.
d876f532	16
41bdb6e2 AJ	17	You should have received a copy of the GNU Lesser General Public
41bdb6e2 AJ	18	License along with the GNU C Library; if not, write to the Free
a334319f UD	19	Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
a334319f UD	20	02111-1307 USA. */
d876f532 UD	21
	22	#if _FP_W_TYPE_SIZE < 32
	23	#error "Here's a nickel, kid. Go buy yourself a real computer."
	24	#endif
	25
	26	#if _FP_W_TYPE_SIZE < 64
	27	#define _FP_FRACTBITS_E (4*_FP_W_TYPE_SIZE)
	28	#else
	29	#define _FP_FRACTBITS_E (2*_FP_W_TYPE_SIZE)
	30	#endif
	31
	32	#define _FP_FRACBITS_E 64
	33	#define _FP_FRACXBITS_E (_FP_FRACTBITS_E - _FP_FRACBITS_E)
	34	#define _FP_WFRACBITS_E (_FP_WORKBITS + _FP_FRACBITS_E)
	35	#define _FP_WFRACXBITS_E (_FP_FRACTBITS_E - _FP_WFRACBITS_E)
	36	#define _FP_EXPBITS_E 15
	37	#define _FP_EXPBIAS_E 16383
	38	#define _FP_EXPMAX_E 32767
	39
	40	#define _FP_QNANBIT_E \
	41	((_FP_W_TYPE)1 << (_FP_FRACBITS_E-2) % _FP_W_TYPE_SIZE)
	42	#define _FP_IMPLBIT_E \
	43	((_FP_W_TYPE)1 << (_FP_FRACBITS_E-1) % _FP_W_TYPE_SIZE)
	44	#define _FP_OVERFLOW_E \
	45	((_FP_W_TYPE)1 << (_FP_WFRACBITS_E % _FP_W_TYPE_SIZE))
	46
	47	#if _FP_W_TYPE_SIZE < 64
	48
	49	union _FP_UNION_E
	50	{
a334319f	51	long double flt;
d876f532 UD	52	struct
	53	{
	54	#if __BYTE_ORDER == __BIG_ENDIAN
	55	unsigned long pad1 : _FP_W_TYPE_SIZE;
	56	unsigned long pad2 : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
	57	unsigned long sign : 1;
	58	unsigned long exp : _FP_EXPBITS_E;
	59	unsigned long frac1 : _FP_W_TYPE_SIZE;
	60	unsigned long frac0 : _FP_W_TYPE_SIZE;
	61	#else
	62	unsigned long frac0 : _FP_W_TYPE_SIZE;
	63	unsigned long frac1 : _FP_W_TYPE_SIZE;
	64	unsigned exp : _FP_EXPBITS_E;
	65	unsigned sign : 1;
	66	#endif /* not bigendian */
	67	} bits __attribute__((packed));
	68	};
	69
	70
	71	#define FP_DECL_E(X) _FP_DECL(4,X)
	72
	73	#define FP_UNPACK_RAW_E(X, val) \
	74	do { \
	75	union _FP_UNION_E _flo; _flo.flt = (val); \
	76	\
	77	X##_f[2] = 0; X##_f[3] = 0; \
	78	X##_f[0] = _flo.bits.frac0; \
	79	X##_f[1] = _flo.bits.frac1; \
	80	X##_e = _flo.bits.exp; \
	81	X##_s = _flo.bits.sign; \
	82	if (!X##_e && (X##_f[1] \|\| X##_f[0]) \
	83	&& !(X##_f[1] & _FP_IMPLBIT_E)) \
	84	{ \
	85	X##_e++; \
	86	FP_SET_EXCEPTION(FP_EX_DENORM); \
	87	} \
	88	} while (0)
	89
	90	#define FP_UNPACK_RAW_EP(X, val) \
	91	do { \
	92	union _FP_UNION_E *_flo = \
	93	(union _FP_UNION_E *)(val); \
	94	\
	95	X##_f[2] = 0; X##_f[3] = 0; \
	96	X##_f[0] = _flo->bits.frac0; \
	97	X##_f[1] = _flo->bits.frac1; \
	98	X##_e = _flo->bits.exp; \
	99	X##_s = _flo->bits.sign; \
	100	if (!X##_e && (X##_f[1] \|\| X##_f[0]) \
	101	&& !(X##_f[1] & _FP_IMPLBIT_E)) \
	102	{ \
	103	X##_e++; \
	104	FP_SET_EXCEPTION(FP_EX_DENORM); \
	105	} \
	106	} while (0)
	107
	108	#define FP_PACK_RAW_E(val, X) \
	109	do { \
	110	union _FP_UNION_E _flo; \
	111	\
	112	if (X##_e) X##_f[1] \|= _FP_IMPLBIT_E; \
	113	else X##_f[1] &= ~(_FP_IMPLBIT_E); \
	114	_flo.bits.frac0 = X##_f[0]; \
	115	_flo.bits.frac1 = X##_f[1]; \
116	_flo.bits.exp = X##_e; \
117	_flo.bits.sign = X##_s; \
118	\
119	(val) = _flo.flt; \
120	} while (0)
121
122	#define FP_PACK_RAW_EP(val, X) \
123	do { \
124	if (!FP_INHIBIT_RESULTS) \
125	{ \
126	union _FP_UNION_E *_flo = \
127	(union _FP_UNION_E *)(val); \
128	\
129	if (X##_e) X##_f[1] \|= _FP_IMPLBIT_E; \
130	else X##_f[1] &= ~(_FP_IMPLBIT_E); \
131	_flo->bits.frac0 = X##_f[0]; \
132	_flo->bits.frac1 = X##_f[1]; \
133	_flo->bits.exp = X##_e; \
134	_flo->bits.sign = X##_s; \
135	} \
136	} while (0)
137
138	#define FP_UNPACK_E(X,val) \
139	do { \
140	FP_UNPACK_RAW_E(X,val); \
141	_FP_UNPACK_CANONICAL(E,4,X); \
142	} while (0)
143
144	#define FP_UNPACK_EP(X,val) \
145	do { \
a334319f	146	FP_UNPACK_RAW_2_P(X,val); \
d876f532 UD	147	_FP_UNPACK_CANONICAL(E,4,X); \
	148	} while (0)
	149
	150	#define FP_PACK_E(val,X) \
	151	do { \
	152	_FP_PACK_CANONICAL(E,4,X); \
	153	FP_PACK_RAW_E(val,X); \
	154	} while (0)
	155
	156	#define FP_PACK_EP(val,X) \
	157	do { \
	158	_FP_PACK_CANONICAL(E,4,X); \
	159	FP_PACK_RAW_EP(val,X); \
	160	} while (0)
	161
	162	#define FP_ISSIGNAN_E(X) _FP_ISSIGNAN(E,4,X)
	163	#define FP_NEG_E(R,X) _FP_NEG(E,4,R,X)
	164	#define FP_ADD_E(R,X,Y) _FP_ADD(E,4,R,X,Y)
	165	#define FP_SUB_E(R,X,Y) _FP_SUB(E,4,R,X,Y)
	166	#define FP_MUL_E(R,X,Y) _FP_MUL(E,4,R,X,Y)
	167	#define FP_DIV_E(R,X,Y) _FP_DIV(E,4,R,X,Y)
	168	#define FP_SQRT_E(R,X) _FP_SQRT(E,4,R,X)
	169
	170	/*
	171	* Square root algorithms:
	172	* We have just one right now, maybe Newton approximation
	173	* should be added for those machines where division is fast.
	174	* This has special _E version because standard _4 square
	175	* root would not work (it has to start normally with the
	176	* second word and not the first), but as we have to do it
	177	* anyway, we optimize it by doing most of the calculations
	178	* in two UWtype registers instead of four.
	179	*/
	180
	181	#define _FP_SQRT_MEAT_E(R, S, T, X, q) \
	182	do { \
	183	q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1); \
	184	_FP_FRAC_SRL_4(X, (_FP_WORKBITS)); \
	185	while (q) \
	186	{ \
	187	T##_f[1] = S##_f[1] + q; \
	188	if (T##_f[1] <= X##_f[1]) \
	189	{ \
	190	S##_f[1] = T##_f[1] + q; \
	191	X##_f[1] -= T##_f[1]; \
	192	R##_f[1] += q; \
	193	} \
	194	_FP_FRAC_SLL_2(X, 1); \
	195	q >>= 1; \
	196	} \
	197	q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1); \
	198	while (q) \
	199	{ \
	200	T##_f[0] = S##_f[0] + q; \
	201	T##_f[1] = S##_f[1]; \
	202	if (T##_f[1] < X##_f[1] \|\| \
	203	(T##_f[1] == X##_f[1] && \
	204	T##_f[0] <= X##_f[0])) \
	205	{ \
	206	S##_f[0] = T##_f[0] + q; \
	207	S##_f[1] += (T##_f[0] > S##_f[0]); \
	208	_FP_FRAC_DEC_2(X, T); \
	209	R##_f[0] += q; \
	210	} \
211	_FP_FRAC_SLL_2(X, 1); \
212	q >>= 1; \
213	} \
214	_FP_FRAC_SLL_4(R, (_FP_WORKBITS)); \
215	if (X##_f[0] \| X##_f[1]) \
216	{ \
217	if (S##_f[1] < X##_f[1] \|\| \
218	(S##_f[1] == X##_f[1] && \
219	S##_f[0] < X##_f[0])) \
220	R##_f[0] \|= _FP_WORK_ROUND; \
221	R##_f[0] \|= _FP_WORK_STICKY; \
222	} \
223	} while (0)
224
225	#define FP_CMP_E(r,X,Y,un) _FP_CMP(E,4,r,X,Y,un)
226	#define FP_CMP_EQ_E(r,X,Y) _FP_CMP_EQ(E,4,r,X,Y)
e7b8c7bc	227	#define FP_CMP_UNORD_E(r,X,Y) _FP_CMP_UNORD(E,4,r,X,Y)
d876f532 UD	228
	229	#define FP_TO_INT_E(r,X,rsz,rsg) _FP_TO_INT(E,4,r,X,rsz,rsg)
	230	#define FP_FROM_INT_E(X,r,rs,rt) _FP_FROM_INT(E,4,X,r,rs,rt)
	231
	232	#define _FP_FRAC_HIGH_E(X) (X##_f[2])
	233	#define _FP_FRAC_HIGH_RAW_E(X) (X##_f[1])
	234
	235	#else /* not _FP_W_TYPE_SIZE < 64 */
	236	union _FP_UNION_E
	237	{
a334319f	238	long double flt /* __attribute__((mode(TF))) */ ;
d876f532 UD	239	struct {
	240	#if __BYTE_ORDER == __BIG_ENDIAN
	241	unsigned long pad : (_FP_W_TYPE_SIZE - 1 - _FP_EXPBITS_E);
	242	unsigned sign : 1;
	243	unsigned exp : _FP_EXPBITS_E;
	244	unsigned long frac : _FP_W_TYPE_SIZE;
	245	#else
	246	unsigned long frac : _FP_W_TYPE_SIZE;
	247	unsigned exp : _FP_EXPBITS_E;
	248	unsigned sign : 1;
	249	#endif
	250	} bits;
	251	};
	252
	253	#define FP_DECL_E(X) _FP_DECL(2,X)
	254
	255	#define FP_UNPACK_RAW_E(X, val) \
	256	do { \
	257	union _FP_UNION_E _flo; _flo.flt = (val); \
	258	\
	259	X##_f0 = _flo.bits.frac; \
	260	X##_f1 = 0; \
	261	X##_e = _flo.bits.exp; \
	262	X##_s = _flo.bits.sign; \
	263	if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E)) \
	264	{ \
	265	X##_e++; \
	266	FP_SET_EXCEPTION(FP_EX_DENORM); \
	267	} \
	268	} while (0)
	269
	270	#define FP_UNPACK_RAW_EP(X, val) \
	271	do { \
	272	union _FP_UNION_E *_flo = \
	273	(union _FP_UNION_E *)(val); \
	274	\
	275	X##_f0 = _flo->bits.frac; \
	276	X##_f1 = 0; \
	277	X##_e = _flo->bits.exp; \
	278	X##_s = _flo->bits.sign; \
	279	if (!X##_e && X##_f0 && !(X##_f0 & _FP_IMPLBIT_E)) \
	280	{ \
	281	X##_e++; \
	282	FP_SET_EXCEPTION(FP_EX_DENORM); \
	283	} \
	284	} while (0)
	285
	286	#define FP_PACK_RAW_E(val, X) \
	287	do { \
	288	union _FP_UNION_E _flo; \
	289	\
	290	if (X##_e) X##_f0 \|= _FP_IMPLBIT_E; \
	291	else X##_f0 &= ~(_FP_IMPLBIT_E); \
	292	_flo.bits.frac = X##_f0; \
	293	_flo.bits.exp = X##_e; \
	294	_flo.bits.sign = X##_s; \
	295	\
	296	(val) = _flo.flt; \
	297	} while (0)
	298
	299	#define FP_PACK_RAW_EP(fs, val, X) \
	300	do { \
	301	if (!FP_INHIBIT_RESULTS) \
	302	{ \
303	union _FP_UNION_E *_flo = \
304	(union _FP_UNION_E *)(val); \
305	\
306	if (X##_e) X##_f0 \|= _FP_IMPLBIT_E; \
307	else X##_f0 &= ~(_FP_IMPLBIT_E); \
308	_flo->bits.frac = X##_f0; \
309	_flo->bits.exp = X##_e; \
310	_flo->bits.sign = X##_s; \
311	} \
312	} while (0)
313
314
315	#define FP_UNPACK_E(X,val) \
316	do { \
317	FP_UNPACK_RAW_E(X,val); \
318	_FP_UNPACK_CANONICAL(E,2,X); \
319	} while (0)
320
321	#define FP_UNPACK_EP(X,val) \
322	do { \
323	FP_UNPACK_RAW_EP(X,val); \
324	_FP_UNPACK_CANONICAL(E,2,X); \
325	} while (0)
326
327	#define FP_PACK_E(val,X) \
328	do { \
329	_FP_PACK_CANONICAL(E,2,X); \
330	FP_PACK_RAW_E(val,X); \
331	} while (0)
332
333	#define FP_PACK_EP(val,X) \
334	do { \
335	_FP_PACK_CANONICAL(E,2,X); \
336	FP_PACK_RAW_EP(val,X); \
337	} while (0)
338
339	#define FP_ISSIGNAN_E(X) _FP_ISSIGNAN(E,2,X)
340	#define FP_NEG_E(R,X) _FP_NEG(E,2,R,X)
341	#define FP_ADD_E(R,X,Y) _FP_ADD(E,2,R,X,Y)
342	#define FP_SUB_E(R,X,Y) _FP_SUB(E,2,R,X,Y)
343	#define FP_MUL_E(R,X,Y) _FP_MUL(E,2,R,X,Y)
344	#define FP_DIV_E(R,X,Y) _FP_DIV(E,2,R,X,Y)
345	#define FP_SQRT_E(R,X) _FP_SQRT(E,2,R,X)
346
347	/*
348	* Square root algorithms:
349	* We have just one right now, maybe Newton approximation
350	* should be added for those machines where division is fast.
351	* We optimize it by doing most of the calculations
352	* in one UWtype registers instead of two, although we don't
353	* have to.
354	*/
355	#define _FP_SQRT_MEAT_E(R, S, T, X, q) \
356	do { \
357	q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1); \
358	_FP_FRAC_SRL_2(X, (_FP_WORKBITS)); \
359	while (q) \
360	{ \
361	T##_f0 = S##_f0 + q; \
362	if (T##_f0 <= X##_f0) \
363	{ \
364	S##_f0 = T##_f0 + q; \
365	X##_f0 -= T##_f0; \
366	R##_f0 += q; \
367	} \
368	_FP_FRAC_SLL_1(X, 1); \
369	q >>= 1; \
370	} \
371	_FP_FRAC_SLL_2(R, (_FP_WORKBITS)); \
372	if (X##_f0) \
373	{ \
374	if (S##_f0 < X##_f0) \
375	R##_f0 \|= _FP_WORK_ROUND; \
376	R##_f0 \|= _FP_WORK_STICKY; \
377	} \
378	} while (0)
379
380	#define FP_CMP_E(r,X,Y,un) _FP_CMP(E,2,r,X,Y,un)
381	#define FP_CMP_EQ_E(r,X,Y) _FP_CMP_EQ(E,2,r,X,Y)
1e832e37	382	#define FP_CMP_UNORD_E(r,X,Y) _FP_CMP_UNORD(E,2,r,X,Y)
d876f532 UD	383
	384	#define FP_TO_INT_E(r,X,rsz,rsg) _FP_TO_INT(E,2,r,X,rsz,rsg)
	385	#define FP_FROM_INT_E(X,r,rs,rt) _FP_FROM_INT(E,2,X,r,rs,rt)
	386
	387	#define _FP_FRAC_HIGH_E(X) (X##_f1)
	388	#define _FP_FRAC_HIGH_RAW_E(X) (X##_f0)
	389
	390	#endif /* not _FP_W_TYPE_SIZE < 64 */