include/math-emu/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 Library General Public License as
   published by the Free Software Foundation; either version 2 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
   Library General Public License for more details.

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


#ifndef    __MATH_EMU_EXTENDED_H__
#define    __MATH_EMU_EXTENDED_H__

#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_TO_INT_E(r,X,rsz,rsg)	_FP_TO_INT(E,4,r,X,rsz,rsg)
#define FP_TO_INT_ROUND_E(r,X,rsz,rsg)	_FP_TO_INT_ROUND(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_TO_INT_E(r,X,rsz,rsg)	_FP_TO_INT(E,2,r,X,rsz,rsg)
#define FP_TO_INT_ROUND_E(r,X,rsz,rsg)	_FP_TO_INT_ROUND(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 */

#endif /* __MATH_EMU_EXTENDED_H__ */