code:
/*	$OpenBSD: cnv_float.h,v 1.9 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)cnv_float.h: Revision: 2.9.88.1 Date: 93/12/07 15:05:29 */

/*
 * Some more constants
 */
#define SGL_FX_MAX_EXP 30
#define DBL_FX_MAX_EXP 62
#define QUAD_FX_MAX_EXP 126

#define Dintp1(object) (object)
#define Dintp2(object) (object)

#define Qintp0(object) (object)
#define Qintp1(object) (object)
#define Qintp2(object) (object)
#define Qintp3(object) (object)

/*
 * These macros will be used specifically by the convert instructions.
 *
 *
 * Single format macros
 */

#define Sgl_to_dbl_exponent(src_exponent,dest)			\
    Deposit_dexponent(dest,src_exponent+(DBL_BIAS-SGL_BIAS))

#define Sgl_to_dbl_mantissa(src_mantissa,destA,destB)	\
    Deposit_dmantissap1(destA,src_mantissa>>3);		\
    Dmantissap2(destB) = src_mantissa << 29

#define Sgl_isinexact_to_fix(sgl_value,exponent)	\
    ((exponent < (SGL_P - 1)) ?				\
     (Sall(sgl_value) << (SGL_EXP_LENGTH + 1 + exponent)) : FALSE)

#define Int_isinexact_to_sgl(int_value)	(int_value << (33 - SGL_EXP_LENGTH))

#define Sgl_roundnearest_from_int(int_value,sgl_value)			\
    if (int_value & 1<<(SGL_EXP_LENGTH - 2))   /* round bit */		\
	if ((int_value << (34 - SGL_EXP_LENGTH)) || Slow(sgl_value))	\
		Sall(sgl_value)++

#define Dint_isinexact_to_sgl(dint_valueA,dint_valueB)		\
    ((Dintp1(dint_valueA) << (33 - SGL_EXP_LENGTH)) || Dintp2(dint_valueB))

#define Sgl_roundnearest_from_dint(dint_valueA,dint_valueB,sgl_value)	\
    if (Dintp1(dint_valueA) & 1<<(SGL_EXP_LENGTH - 2))			\
	if ((Dintp1(dint_valueA) << (34 - SGL_EXP_LENGTH)) ||		\
	Dintp2(dint_valueB) || Slow(sgl_value)) Sall(sgl_value)++

#define Dint_isinexact_to_dbl(dint_value)	\
    (Dintp2(dint_value) << (33 - DBL_EXP_LENGTH))

#define Dbl_roundnearest_from_dint(dint_opndB,dbl_opndA,dbl_opndB)	\
    if (Dintp2(dint_opndB) & 1<<(DBL_EXP_LENGTH - 2))			\
       if ((Dintp2(dint_opndB) << (34 - DBL_EXP_LENGTH)) || Dlowp2(dbl_opndB))  \
          if ((++Dallp2(dbl_opndB))==0) Dallp1(dbl_opndA)++

#define Sgl_isone_roundbit(sgl_value,exponent)			\
    ((Sall(sgl_value) << (SGL_EXP_LENGTH + 1 + exponent)) >> 31)

#define Sgl_isone_stickybit(sgl_value,exponent)		\
    (exponent < (SGL_P - 2) ?				\
     Sall(sgl_value) << (SGL_EXP_LENGTH + 2 + exponent) : FALSE)

/*
 * Double format macros
 */

#define Dbl_to_sgl_exponent(src_exponent,dest)			\
    dest = src_exponent + (SGL_BIAS - DBL_BIAS)

#define Dbl_to_sgl_mantissa(srcA,srcB,dest,inexact,guard,sticky,odd)	\
    Shiftdouble(Dmantissap1(srcA),Dmantissap2(srcB),29,dest);	\
    guard = Dbit3p2(srcB);					\
    sticky = Dallp2(srcB)<<4;					\
    inexact = guard | sticky;					\
    odd = Dbit2p2(srcB)

#define Dbl_to_sgl_denormalized(srcA,srcB,exp,dest,inexact,guard,sticky,odd,tiny) \
    Deposit_dexponent(srcA,1);						\
    tiny = TRUE;							\
    if (exp >= -2) {							\
	if (exp == 0) {							\
	    inexact = Dallp2(srcB) << 3;				\
	    guard = inexact >> 31;					\
	    sticky = inexact << 1;					\
	    Shiftdouble(Dmantissap1(srcA),Dmantissap2(srcB),29,dest);	\
	    odd = dest << 31;						\
	    if (inexact) {						\
		switch(Rounding_mode()) {				\
		    case ROUNDPLUS:					\
			if (Dbl_iszero_sign(srcA)) {			\
			    dest++;					\
			    if (Sgl_isone_hidden(dest))	\
				tiny = FALSE;				\
			    dest--;					\
			}						\
			break;						\
		    case ROUNDMINUS:					\
			if (Dbl_isone_sign(srcA)) {			\
			    dest++;					\
			    if (Sgl_isone_hidden(dest))	\
				tiny = FALSE;				\
			    dest--;					\
			}						\
			break;						\
		    case ROUNDNEAREST:					\
			if (guard && (sticky || odd)) {			\
			    dest++;					\
			    if (Sgl_isone_hidden(dest))	\
				tiny = FALSE;				\
			    dest--;					\
			}						\
			break;						\
		}							\
	    }								\
		/* shift right by one to get correct result */		\
		guard = odd;						\
		sticky = inexact;					\
		inexact |= guard;					\
		dest >>= 1;						\
		Deposit_dsign(srcA,0);					\
	        Shiftdouble(Dallp1(srcA),Dallp2(srcB),30,dest);		\
	        odd = dest << 31;					\
	}								\
	else {								\
	    inexact = Dallp2(srcB) << (2 + exp);			\
	    guard = inexact >> 31;					\
	    sticky = inexact << 1;					\
	    Deposit_dsign(srcA,0);					\
	    if (exp == -2) dest = Dallp1(srcA);				\
	    else Variable_shift_double(Dallp1(srcA),Dallp2(srcB),30-exp,dest); \
	    odd = dest << 31;						\
	}								\
    }									\
    else {								\
	Deposit_dsign(srcA,0);						\
	if (exp > (1 - SGL_P)) {					\
	    dest = Dallp1(srcA) >> (- 2 - exp);				\
	    inexact = Dallp1(srcA) << (34 + exp);			\
	    guard = inexact >> 31;					\
	    sticky = (inexact << 1) | Dallp2(srcB);			\
	    inexact |= Dallp2(srcB);					\
	    odd = dest << 31;						\
	}								\
	else {								\
	    dest = 0;							\
	    inexact = Dallp1(srcA) | Dallp2(srcB);			\
	    if (exp == (1 - SGL_P)) {					\
		guard = Dhidden(srcA);					\
		sticky = Dmantissap1(srcA) | Dallp2(srcB);		\
	    }								\
	    else {							\
		guard = 0;						\
		sticky = inexact;					\
	    }								\
	    odd = 0;							\
	}								\
    }									\
    exp = 0

#define Dbl_isinexact_to_fix(dbl_valueA,dbl_valueB,exponent)		\
    (exponent < (DBL_P-33) ?						\
     Dallp2(dbl_valueB) || Dallp1(dbl_valueA) << (DBL_EXP_LENGTH+1+exponent) : \
     (exponent < (DBL_P-1) ? Dallp2(dbl_valueB) << (exponent + (33-DBL_P)) :   \
      FALSE))

#define Dbl_isoverflow_to_int(exponent,dbl_valueA,dbl_valueB)		\
    ((exponent > SGL_FX_MAX_EXP + 1) || Dsign(dbl_valueA)==0 ||		\
     Dmantissap1(dbl_valueA)!=0 || (Dallp2(dbl_valueB)>>21)!=0 )

#define Dbl_isone_roundbit(dbl_valueA,dbl_valueB,exponent)              \
    ((exponent < (DBL_P - 33) ?						\
      Dallp1(dbl_valueA) >> ((30 - DBL_EXP_LENGTH) - exponent) :	\
      Dallp2(dbl_valueB) >> ((DBL_P - 2) - exponent)) & 1)

#define Dbl_isone_stickybit(dbl_valueA,dbl_valueB,exponent)		\
    (exponent < (DBL_P-34) ?						\
     (Dallp2(dbl_valueB) || Dallp1(dbl_valueA)<<(DBL_EXP_LENGTH+2+exponent)) : \
     (exponent<(DBL_P-2) ? (Dallp2(dbl_valueB) << (exponent + (34-DBL_P))) : \
      FALSE))

/* Int macros */

#define Int_from_sgl_mantissa(sgl_value,exponent)	\
    Sall(sgl_value) =				\
	(unsigned)(Sall(sgl_value) << SGL_EXP_LENGTH)>>(31 - exponent)

#define Int_from_dbl_mantissa(dbl_valueA,dbl_valueB,exponent)	\
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),22,Dallp1(dbl_valueA)); \
    if (exponent < 31) Dallp1(dbl_valueA) >>= 30 - exponent;	\
    else Dallp1(dbl_valueA) <<= 1

#define Int_negate(int_value) int_value = -int_value

/* Dint macros */

#define Dint_from_sgl_mantissa(sgl_value,exponent,dresultA,dresultB)	\
    {Sall(sgl_value) <<= SGL_EXP_LENGTH;  /*  left-justify  */		\
    if (exponent <= 31) {						\
	Dintp1(dresultA) = 0;						\
	Dintp2(dresultB) = (unsigned)Sall(sgl_value) >> (31 - exponent); \
    }									\
    else {								\
	Dintp1(dresultA) = Sall(sgl_value) >> (63 - exponent);		\
	Dintp2(dresultB) = Sall(sgl_value) << (exponent - 31);		\
    }}

#define Dint_from_dbl_mantissa(dbl_valueA,dbl_valueB,exponent,destA,destB) \
    {if (exponent < 32) {						\
	Dintp1(destA) = 0;						\
	if (exponent <= 20)						\
	    Dintp2(destB) = Dallp1(dbl_valueA) >> (20-exponent);	\
	else Variable_shift_double(Dallp1(dbl_valueA),Dallp2(dbl_valueB), \
	     52-exponent,Dintp2(destB));					\
    }									\
    else {								\
	if (exponent <= 52) {						\
	    Dintp1(destA) = Dallp1(dbl_valueA) >> (52-exponent);	\
	    if (exponent == 52) Dintp2(destB) = Dallp2(dbl_valueB);	\
	    else Variable_shift_double(Dallp1(dbl_valueA),Dallp2(dbl_valueB), \
	    52-exponent,Dintp2(destB));					\
        }								\
	else {								\
	    Variable_shift_double(Dallp1(dbl_valueA),Dallp2(dbl_valueB), \
	    84-exponent,Dintp1(destA));					\
	    Dintp2(destB) = Dallp2(dbl_valueB) << (exponent-52);	\
	}								\
    }}

#define Dint_setzero(dresultA,dresultB)		\
    Dintp1(dresultA) = 0;			\
    Dintp2(dresultB) = 0

#define Dint_setone_sign(dresultA,dresultB)		\
    Dintp1(dresultA) = ~Dintp1(dresultA);		\
    if ((Dintp2(dresultB) = -Dintp2(dresultB)) == 0) Dintp1(dresultA)++

#define Dint_set_minint(dresultA,dresultB)		\
    Dintp1(dresultA) = 1<<31;				\
    Dintp2(dresultB) = 0

#define Dint_isone_lowp2(dresultB)  (Dintp2(dresultB) & 01)

#define Dint_increment(dresultA,dresultB)		\
    if ((++Dintp2(dresultB))==0) Dintp1(dresultA)++

#define Dint_decrement(dresultA,dresultB)		\
    if ((Dintp2(dresultB)--)==0) Dintp1(dresultA)--

#define Dint_negate(dresultA,dresultB)			\
    Dintp1(dresultA) = ~Dintp1(dresultA);		\
    if ((Dintp2(dresultB) = -Dintp2(dresultB))==0) Dintp1(dresultA)++

#define Dint_copyfromptr(src,destA,destB) \
     Dintp1(destA) = src->wd0;		\
     Dintp2(destB) = src->wd1
#define Dint_copytoptr(srcA,srcB,dest)	\
    dest->wd0 = Dintp1(srcA);		\
    dest->wd1 = Dintp2(srcB)

/* other macros  */

#define Find_ms_one_bit(value, position)	\
    {						\
	int var;				\
	for (var = 8; var >= 1; var >>= 1) {	\
	    if (value >> (32 - position))	\
		position -= var;		\
		else position += var;		\
	}					\
	if ((value >> (32 - position)) == 0)	\
	    position--;				\
	else position -= 2;			\
    }

int sgl_to_sgl_fcnvfx(sgl_floating_point *, sgl_floating_point *, sgl_floating_point *, unsigned int *);
int sgl_to_dbl_fcnvfx(sgl_floating_point *, sgl_floating_point *, dbl_integer *, unsigned int *);
int dbl_to_sgl_fcnvfx(dbl_floating_point *, dbl_floating_point *, int *, unsigned int *);
int dbl_to_dbl_fcnvfx(dbl_floating_point *, dbl_floating_point *, dbl_integer *, unsigned int *);

int sgl_to_sgl_fcnvfxt(sgl_floating_point *, sgl_floating_point *, int *, unsigned int *);
int sgl_to_dbl_fcnvfxt(sgl_floating_point *, sgl_floating_point *, dbl_integer *, unsigned int *);
int dbl_to_sgl_fcnvfxt(dbl_floating_point *, dbl_floating_point *, int *, unsigned int *);
int dbl_to_dbl_fcnvfxt(dbl_floating_point *, dbl_floating_point *, dbl_integer *, unsigned int *);

int sgl_to_sgl_fcnvxf(int *, int *, sgl_floating_point *, unsigned int *);
int sgl_to_dbl_fcnvxf(int *, int *, dbl_floating_point *, unsigned int *);
int dbl_to_sgl_fcnvxf(dbl_integer *, dbl_integer *, sgl_floating_point *, unsigned int *);
int dbl_to_dbl_fcnvxf(dbl_integer *, dbl_integer *, dbl_floating_point *, unsigned int *);


code:
/*	$OpenBSD: dbl_float.h,v 1.12 2013/11/26 20:33:12 deraadt Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)dbl_float.h: Revision: 2.9.88.1 Date: 93/12/07 15:05:32 */

/**************************************
 * Declare double precision functions *
 **************************************/

/* 32-bit word grabing functions */
#define Dbl_firstword(value) Dallp1(value)
#define Dbl_secondword(value) Dallp2(value)
#define Dbl_thirdword(value) dummy_location
#define Dbl_fourthword(value) dummy_location

#define Dbl_sign(object) Dsign(object)
#define Dbl_exponent(object) Dexponent(object)
#define Dbl_signexponent(object) Dsignexponent(object)
#define Dbl_mantissap1(object) Dmantissap1(object)
#define Dbl_mantissap2(object) Dmantissap2(object)
#define Dbl_exponentmantissap1(object) Dexponentmantissap1(object)
#define Dbl_allp1(object) Dallp1(object)
#define Dbl_allp2(object) Dallp2(object)

/* dbl_and_signs ands the sign bits of each argument and puts the result
 * into the first argument. dbl_or_signs ors those same sign bits */
#define Dbl_and_signs( src1dst, src2)		\
    Dallp1(src1dst) = (Dallp1(src2)|~(1<<31)) & Dallp1(src1dst)
#define Dbl_or_signs( src1dst, src2)		\
    Dallp1(src1dst) = (Dallp1(src2)&(1<<31)) | Dallp1(src1dst)

/* The hidden bit is always the low bit of the exponent */
#define Dbl_clear_exponent_set_hidden(srcdst) Deposit_dexponent(srcdst,1)
#define Dbl_clear_signexponent_set_hidden(srcdst) \
    Deposit_dsignexponent(srcdst,1)
#define Dbl_clear_sign(srcdst) Dallp1(srcdst) &= ~(1<<31)
#define Dbl_clear_signexponent(srcdst) \
    Dallp1(srcdst) &= Dmantissap1((unsigned)-1)

/* Exponent field for doubles has already been cleared and may be
 * included in the shift.  Here we need to generate two double width
 * variable shifts.  The insignificant bits can be ignored.
 *      MTSAR f(varamount)
 *      VSHD	srcdst.high,srcdst.low => srcdst.low
 *	VSHD	0,srcdst.high => srcdst.high
 * This is very difficult to model with C expressions since the shift amount
 * could exceed 32.  */
/* varamount must be less than 64 */
#define Dbl_rightshift(srcdstA, srcdstB, varamount)			\
    {if((varamount) >= 32) {						\
	Dallp2(srcdstB) = Dallp1(srcdstA) >> (varamount-32);		\
	Dallp1(srcdstA)=0;						\
    }									\
    else if(varamount > 0) {						\
	Variable_shift_double(Dallp1(srcdstA), Dallp2(srcdstB),		\
	  (varamount), Dallp2(srcdstB));				\
	Dallp1(srcdstA) >>= varamount;					\
    } }
/* varamount must be less than 64 */
#define Dbl_rightshift_exponentmantissa(srcdstA, srcdstB, varamount)	\
    {if((varamount) >= 32) {						\
	Dallp2(srcdstB) = Dexponentmantissap1(srcdstA) >> ((varamount)-32); \
	Dallp1(srcdstA) &= (1<<31);  /* clear exponentmantissa field */ \
    }									\
    else if(varamount > 0) {						\
	Variable_shift_double(Dexponentmantissap1(srcdstA), Dallp2(srcdstB), \
	(varamount), Dallp2(srcdstB));					\
	Deposit_dexponentmantissap1(srcdstA,				\
	    (Dexponentmantissap1(srcdstA)>>(varamount)));			\
    } }
/* varamount must be less than 64 */
#define Dbl_leftshift(srcdstA, srcdstB, varamount)			\
    {if((varamount) >= 32) {						\
	Dallp1(srcdstA) = Dallp2(srcdstB) << (varamount-32);		\
	Dallp2(srcdstB)=0;						\
    }									\
    else {								\
	if ((varamount) > 0) {						\
	    Dallp1(srcdstA) = (Dallp1(srcdstA) << (varamount)) |	\
		(Dallp2(srcdstB) >> (32-(varamount)));			\
	    Dallp2(srcdstB) <<= varamount;				\
	}								\
    } }
#define Dbl_leftshiftby1_withextent(lefta,leftb,right,resulta,resultb)	\
    Shiftdouble(Dallp1(lefta), Dallp2(leftb), 31, Dallp1(resulta));	\
    Shiftdouble(Dallp2(leftb), Extall(right), 31, Dallp2(resultb))

#define Dbl_rightshiftby1_withextent(leftb,right,dst)		\
    Extall(dst) = (Dallp2(leftb) << 31) | ((unsigned)Extall(right) >> 1) | \
		  Extlow(right)

#define Dbl_arithrightshiftby1(srcdstA,srcdstB)			\
    Shiftdouble(Dallp1(srcdstA),Dallp2(srcdstB),1,Dallp2(srcdstB));\
    Dallp1(srcdstA) = (int)Dallp1(srcdstA) >> 1

/* Sign extend the sign bit with an integer destination */
#define Dbl_signextendedsign(value)  Dsignedsign(value)

#define Dbl_isone_hidden(dbl_value) (Is_dhidden(dbl_value)!=0)
/* Singles and doubles may include the sign and exponent fields.  The
 * hidden bit and the hidden overflow must be included. */
#define Dbl_increment(dbl_valueA,dbl_valueB) \
    if( (Dallp2(dbl_valueB) += 1) == 0 )  Dallp1(dbl_valueA) += 1
#define Dbl_increment_mantissa(dbl_valueA,dbl_valueB) \
    if( (Dmantissap2(dbl_valueB) += 1) == 0 )  \
    Deposit_dmantissap1(dbl_valueA,dbl_valueA+1)
#define Dbl_decrement(dbl_valueA,dbl_valueB) \
    if( Dallp2(dbl_valueB) == 0 )  Dallp1(dbl_valueA) -= 1; \
    Dallp2(dbl_valueB) -= 1

#define Dbl_isone_sign(dbl_value) (Is_dsign(dbl_value)!=0)
#define Dbl_isone_hiddenoverflow(dbl_value) (Is_dhiddenoverflow(dbl_value)!=0)
#define Dbl_isone_lowmantissap1(dbl_valueA) (Is_dlowp1(dbl_valueA)!=0)
#define Dbl_isone_lowmantissap2(dbl_valueB) (Is_dlowp2(dbl_valueB)!=0)
#define Dbl_isone_signaling(dbl_value) (Is_dsignaling(dbl_value)!=0)
#define Dbl_is_signalingnan(dbl_value) (Dsignalingnan(dbl_value)==0xfff)
#define Dbl_isnotzero(dbl_valueA,dbl_valueB) \
    (Dallp1(dbl_valueA) || Dallp2(dbl_valueB))
#define Dbl_isnotzero_hiddenhigh7mantissa(dbl_value) \
    (Dhiddenhigh7mantissa(dbl_value)!=0)
#define Dbl_isnotzero_exponent(dbl_value) (Dexponent(dbl_value)!=0)
#define Dbl_isnotzero_mantissa(dbl_valueA,dbl_valueB) \
    (Dmantissap1(dbl_valueA) || Dmantissap2(dbl_valueB))
#define Dbl_isnotzero_mantissap1(dbl_valueA) (Dmantissap1(dbl_valueA)!=0)
#define Dbl_isnotzero_mantissap2(dbl_valueB) (Dmantissap2(dbl_valueB)!=0)
#define Dbl_isnotzero_exponentmantissa(dbl_valueA,dbl_valueB) \
    (Dexponentmantissap1(dbl_valueA) || Dmantissap2(dbl_valueB))
#define Dbl_isnotzero_low4p2(dbl_value) (Dlow4p2(dbl_value)!=0)
#define Dbl_iszero(dbl_valueA,dbl_valueB) (Dallp1(dbl_valueA)==0 && \
    Dallp2(dbl_valueB)==0)
#define Dbl_iszero_allp1(dbl_value) (Dallp1(dbl_value)==0)
#define Dbl_iszero_allp2(dbl_value) (Dallp2(dbl_value)==0)
#define Dbl_iszero_hidden(dbl_value) (Is_dhidden(dbl_value)==0)
#define Dbl_iszero_hiddenoverflow(dbl_value) (Is_dhiddenoverflow(dbl_value)==0)
#define Dbl_iszero_hiddenhigh3mantissa(dbl_value) \
    (Dhiddenhigh3mantissa(dbl_value)==0)
#define Dbl_iszero_hiddenhigh7mantissa(dbl_value) \
    (Dhiddenhigh7mantissa(dbl_value)==0)
#define Dbl_iszero_sign(dbl_value) (Is_dsign(dbl_value)==0)
#define Dbl_iszero_exponent(dbl_value) (Dexponent(dbl_value)==0)
#define Dbl_iszero_mantissa(dbl_valueA,dbl_valueB) \
    (Dmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0)
#define Dbl_iszero_exponentmantissa(dbl_valueA,dbl_valueB) \
    (Dexponentmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0)
#define Dbl_isinfinity_exponent(dbl_value)		\
    (Dexponent(dbl_value)==DBL_INFINITY_EXPONENT)
#define Dbl_isnotinfinity_exponent(dbl_value)		\
    (Dexponent(dbl_value)!=DBL_INFINITY_EXPONENT)
#define Dbl_isinfinity(dbl_valueA,dbl_valueB)			\
    (Dexponent(dbl_valueA)==DBL_INFINITY_EXPONENT &&	\
    Dmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0)
#define Dbl_isnan(dbl_valueA,dbl_valueB)		\
    (Dexponent(dbl_valueA)==DBL_INFINITY_EXPONENT &&	\
    (Dmantissap1(dbl_valueA)!=0 || Dmantissap2(dbl_valueB)!=0))
#define Dbl_isnotnan(dbl_valueA,dbl_valueB)		\
    (Dexponent(dbl_valueA)!=DBL_INFINITY_EXPONENT ||	\
    (Dmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0))

#define Dbl_islessthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
    (Dallp1(dbl_op1a) < Dallp1(dbl_op2a) ||			\
     (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
      Dallp2(dbl_op1b) < Dallp2(dbl_op2b)))
#define Dbl_isgreaterthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
    (Dallp1(dbl_op1a) > Dallp1(dbl_op2a) ||			\
     (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
      Dallp2(dbl_op1b) > Dallp2(dbl_op2b)))
#define Dbl_isnotlessthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
    (Dallp1(dbl_op1a) > Dallp1(dbl_op2a) ||			\
     (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
      Dallp2(dbl_op1b) >= Dallp2(dbl_op2b)))
#define Dbl_isnotgreaterthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b) \
    (Dallp1(dbl_op1a) < Dallp1(dbl_op2a) ||			\
     (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
      Dallp2(dbl_op1b) <= Dallp2(dbl_op2b)))
#define Dbl_isequal(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
     ((Dallp1(dbl_op1a) == Dallp1(dbl_op2a)) &&			\
      (Dallp2(dbl_op1b) == Dallp2(dbl_op2b)))

#define Dbl_leftshiftby8(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),24,Dallp1(dbl_valueA)); \
    Dallp2(dbl_valueB) <<= 8
#define Dbl_leftshiftby7(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),25,Dallp1(dbl_valueA)); \
    Dallp2(dbl_valueB) <<= 7
#define Dbl_leftshiftby4(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),28,Dallp1(dbl_valueA)); \
    Dallp2(dbl_valueB) <<= 4
#define Dbl_leftshiftby3(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),29,Dallp1(dbl_valueA)); \
    Dallp2(dbl_valueB) <<= 3
#define Dbl_leftshiftby2(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),30,Dallp1(dbl_valueA)); \
    Dallp2(dbl_valueB) <<= 2
#define Dbl_leftshiftby1(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),31,Dallp1(dbl_valueA)); \
    Dallp2(dbl_valueB) <<= 1

#define Dbl_rightshiftby8(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),8,Dallp2(dbl_valueB)); \
    Dallp1(dbl_valueA) >>= 8
#define Dbl_rightshiftby4(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),4,Dallp2(dbl_valueB)); \
    Dallp1(dbl_valueA) >>= 4
#define Dbl_rightshiftby2(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),2,Dallp2(dbl_valueB)); \
    Dallp1(dbl_valueA) >>= 2
#define Dbl_rightshiftby1(dbl_valueA,dbl_valueB) \
    Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),1,Dallp2(dbl_valueB)); \
    Dallp1(dbl_valueA) >>= 1

/* This magnitude comparison uses the signless first words and
 * the regular part2 words.  The comparison is graphically:
 *
 *       1st greater?  -------------
 *				   |
 *       1st less?-----------------+---------
 *				   |	    |
 *       2nd greater or equal----->|	    |
 *				 False     True
 */
#define Dbl_ismagnitudeless(leftB,rightB,signlessleft,signlessright)	\
      ((signlessleft <= signlessright) &&				\
       ( (signlessleft < signlessright) || (Dallp2(leftB)<Dallp2(rightB)) ))

#define Dbl_copytoint_exponentmantissap1(src,dest) \
    dest = Dexponentmantissap1(src)

/* A quiet NaN has the high mantissa bit clear and at least on other (in this
 * case the adjacent bit) bit set. */
#define Dbl_set_quiet(dbl_value) Deposit_dhigh2mantissa(dbl_value,1)
#define Dbl_set_exponent(dbl_value, exp) Deposit_dexponent(dbl_value,exp)

#define Dbl_set_mantissa(desta,destb,valuea,valueb)	\
    Deposit_dmantissap1(desta,valuea);			\
    Dmantissap2(destb) = Dmantissap2(valueb)
#define Dbl_set_mantissap1(desta,valuea)		\
    Deposit_dmantissap1(desta,valuea)
#define Dbl_set_mantissap2(destb,valueb)		\
    Dmantissap2(destb) = Dmantissap2(valueb)

#define Dbl_set_exponentmantissa(desta,destb,valuea,valueb)	\
    Deposit_dexponentmantissap1(desta,valuea);			\
    Dmantissap2(destb) = Dmantissap2(valueb)
#define Dbl_set_exponentmantissap1(dest,value)			\
    Deposit_dexponentmantissap1(dest,value)

#define Dbl_copyfromptr(src,desta,destb) \
    Dallp1(desta) = src->wd0;		\
    Dallp2(destb) = src->wd1
#define Dbl_copytoptr(srca,srcb,dest)	\
    dest->wd0 = Dallp1(srca);		\
    dest->wd1 = Dallp2(srcb)

/*  An infinity is represented with the max exponent and a zero mantissa */
#define Dbl_setinfinity_exponent(dbl_value) \
    Deposit_dexponent(dbl_value,DBL_INFINITY_EXPONENT)
#define Dbl_setinfinity_exponentmantissa(dbl_valueA,dbl_valueB)	\
    Deposit_dexponentmantissap1(dbl_valueA,			\
    (DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH))));	\
    Dmantissap2(dbl_valueB) = 0
#define Dbl_setinfinitypositive(dbl_valueA,dbl_valueB)		\
    Dallp1(dbl_valueA)						\
	= (DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH)));	\
    Dmantissap2(dbl_valueB) = 0
#define Dbl_setinfinitynegative(dbl_valueA,dbl_valueB)		\
    Dallp1(dbl_valueA) = (1<<31) |				\
	(DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH)));	\
    Dmantissap2(dbl_valueB) = 0
#define Dbl_setinfinity(dbl_valueA,dbl_valueB,sign)		\
    Dallp1(dbl_valueA) = (sign << 31) |				\
	(DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH)));	\
    Dmantissap2(dbl_valueB) = 0

#define Dbl_sethigh4bits(dbl_value, extsign) Deposit_dhigh4p1(dbl_value,extsign)
#define Dbl_set_sign(dbl_value,sign) Deposit_dsign(dbl_value,sign)
#define Dbl_invert_sign(dbl_value) Deposit_dsign(dbl_value,~Dsign(dbl_value))
#define Dbl_setone_sign(dbl_value) Deposit_dsign(dbl_value,1)
#define Dbl_setone_lowmantissap2(dbl_value) Deposit_dlowp2(dbl_value,1)
#define Dbl_setzero_sign(dbl_value) Dallp1(dbl_value) &= 0x7fffffff
#define Dbl_setzero_exponent(dbl_value)			\
    Dallp1(dbl_value) &= 0x800fffff
#define Dbl_setzero_mantissa(dbl_valueA,dbl_valueB)	\
    Dallp1(dbl_valueA) &= 0xfff00000;			\
    Dallp2(dbl_valueB) = 0
#define Dbl_setzero_mantissap1(dbl_value) Dallp1(dbl_value) &= 0xfff00000
#define Dbl_setzero_mantissap2(dbl_value) Dallp2(dbl_value) = 0
#define Dbl_setzero_exponentmantissa(dbl_valueA,dbl_valueB)	\
    Dallp1(dbl_valueA) &= 0x80000000;		\
    Dallp2(dbl_valueB) = 0
#define Dbl_setzero_exponentmantissap1(dbl_valueA)	\
    Dallp1(dbl_valueA) &= 0x80000000
#define Dbl_setzero(dbl_valueA,dbl_valueB) \
    Dallp1(dbl_valueA) = 0; Dallp2(dbl_valueB) = 0
#define Dbl_setzerop1(dbl_value) Dallp1(dbl_value) = 0
#define Dbl_setzerop2(dbl_value) Dallp2(dbl_value) = 0
#define Dbl_setnegativezero(dbl_value) \
    Dallp1(dbl_value) = 1U << 31; Dallp2(dbl_value) = 0
#define Dbl_setnegativezerop1(dbl_value) Dallp1(dbl_value) = 1U << 31

/* Use the following macro for both overflow & underflow conditions */
#define ovfl -
#define unfl +
#define Dbl_setwrapped_exponent(dbl_value,exponent,op) \
    Deposit_dexponent(dbl_value,(exponent op DBL_WRAP))

#define Dbl_setlargestpositive(dbl_valueA,dbl_valueB)			\
    Dallp1(dbl_valueA) = ((DBL_EMAX+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH))) \
			| ((1<<(32-(1+DBL_EXP_LENGTH))) - 1 );		\
    Dallp2(dbl_valueB) = 0xFFFFFFFF
#define Dbl_setlargestnegative(dbl_valueA,dbl_valueB)			\
    Dallp1(dbl_valueA) = ((DBL_EMAX+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH))) \
			| ((1<<(32-(1+DBL_EXP_LENGTH))) - 1 ) | (1<<31); \
    Dallp2(dbl_valueB) = 0xFFFFFFFF
#define Dbl_setlargest_exponentmantissa(dbl_valueA,dbl_valueB)		\
    Deposit_dexponentmantissap1(dbl_valueA,				\
	(((DBL_EMAX+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH)))		\
			| ((1<<(32-(1+DBL_EXP_LENGTH))) - 1 )));	\
    Dallp2(dbl_valueB) = 0xFFFFFFFF

#define Dbl_setnegativeinfinity(dbl_valueA,dbl_valueB)			\
    Dallp1(dbl_valueA) = ((1<<DBL_EXP_LENGTH) | DBL_INFINITY_EXPONENT)	\
			 << (32-(1+DBL_EXP_LENGTH)) ;			\
    Dallp2(dbl_valueB) = 0
#define Dbl_setlargest(dbl_valueA,dbl_valueB,sign)			\
    Dallp1(dbl_valueA) = (sign << 31) |					\
	((DBL_EMAX+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH))) |		\
	 ((1 << (32-(1+DBL_EXP_LENGTH))) - 1 );				\
    Dallp2(dbl_valueB) = 0xFFFFFFFF

/* The high bit is always zero so arithmetic or logical shifts will work. */
#define Dbl_right_align(srcdstA,srcdstB,shift,extent)			\
    if( shift >= 32 )							\
	{								\
	/* Big shift requires examining the portion shift off		\
	the end to properly set inexact.  */				\
	if(shift < 64)							\
	    {								\
	    if(shift > 32)						\
		{							\
		Variable_shift_double(Dallp1(srcdstA),Dallp2(srcdstB),	\
		 shift-32, Extall(extent));				\
		if(Dallp2(srcdstB) << (64 - (shift))) Ext_setone_low(extent); \
		}							\
	    else Extall(extent) = Dallp2(srcdstB);			\
	    Dallp2(srcdstB) = Dallp1(srcdstA) >> (shift - 32);		\
	    }								\
	else								\
	    {								\
	    Extall(extent) = Dallp1(srcdstA);				\
	    if(Dallp2(srcdstB)) Ext_setone_low(extent);			\
	    Dallp2(srcdstB) = 0;					\
	    }								\
	Dallp1(srcdstA) = 0;						\
	}								\
    else								\
	{								\
	/* Small alignment is simpler.  Extension is easily set. */	\
	if (shift > 0)							\
	    {								\
	    Extall(extent) = Dallp2(srcdstB) << (32 - (shift));		\
	    Variable_shift_double(Dallp1(srcdstA),Dallp2(srcdstB),shift, \
	     Dallp2(srcdstB));						\
	    Dallp1(srcdstA) >>= shift;					\
	    }								\
	else Extall(extent) = 0;					\
	}

/*
 * Here we need to shift the result right to correct for an overshift
 * (due to the exponent becoming negative) during normalization.
 */
#define Dbl_fix_overshift(srcdstA,srcdstB,shift,extent)			\
	    Extall(extent) = Dallp2(srcdstB) << (32 - (shift));		\
	    Dallp2(srcdstB) = (Dallp1(srcdstA) << (32 - (shift))) |	\
		(Dallp2(srcdstB) >> (shift));				\
	    Dallp1(srcdstA) = Dallp1(srcdstA) >> shift

#define Dbl_hiddenhigh3mantissa(dbl_value) Dhiddenhigh3mantissa(dbl_value)
#define Dbl_hidden(dbl_value) Dhidden(dbl_value)
#define Dbl_lowmantissap2(dbl_value) Dlowp2(dbl_value)

/* The left argument is never smaller than the right argument */
#define Dbl_subtract(lefta,leftb,righta,rightb,resulta,resultb)			\
    if( Dallp2(rightb) > Dallp2(leftb) ) Dallp1(lefta)--;	\
    Dallp2(resultb) = Dallp2(leftb) - Dallp2(rightb);		\
    Dallp1(resulta) = Dallp1(lefta) - Dallp1(righta)

/* Subtract right augmented with extension from left augmented with zeros and
 * store into result and extension. */
#define Dbl_subtract_withextension(lefta,leftb,righta,rightb,extent,resulta,resultb)	\
    Dbl_subtract(lefta,leftb,righta,rightb,resulta,resultb);		\
    if( (Extall(extent) = 0-Extall(extent)) )				\
	{								\
	if((Dallp2(resultb)--) == 0) Dallp1(resulta)--;			\
	}

#define Dbl_addition(lefta,leftb,righta,rightb,resulta,resultb)		\
    /* If the sum of the low words is less than either source, then	\
     * an overflow into the next word occurred. */			\
    Dallp1(resulta) = Dallp1(lefta) + Dallp1(righta);			\
    if((Dallp2(resultb) = Dallp2(leftb) + Dallp2(rightb)) < Dallp2(rightb)) \
	Dallp1(resulta)++

#define Dbl_xortointp1(left,right,result)			\
    result = Dallp1(left) XOR Dallp1(right)

#define Dbl_xorfromintp1(left,right,result)			\
    Dallp1(result) = left XOR Dallp1(right)

#define Dbl_swap_lower(left,right)				\
    Dallp2(left)  = Dallp2(left) XOR Dallp2(right);		\
    Dallp2(right) = Dallp2(left) XOR Dallp2(right);		\
    Dallp2(left)  = Dallp2(left) XOR Dallp2(right)

/* Need to Initialize */
#define Dbl_makequietnan(desta,destb)					\
    Dallp1(desta) = ((DBL_EMAX+DBL_BIAS)+1)<< (32-(1+DBL_EXP_LENGTH))	\
		| (1<<(32-(1+DBL_EXP_LENGTH+2)));			\
    Dallp2(destb) = 0
#define Dbl_makesignalingnan(desta,destb)				\
    Dallp1(desta) = ((DBL_EMAX+DBL_BIAS)+1)<< (32-(1+DBL_EXP_LENGTH))	\
		| (1<<(32-(1+DBL_EXP_LENGTH+1)));			\
    Dallp2(destb) = 0

#define Dbl_normalize(dbl_opndA,dbl_opndB,exponent)			\
	while(Dbl_iszero_hiddenhigh7mantissa(dbl_opndA)) {		\
		Dbl_leftshiftby8(dbl_opndA,dbl_opndB);			\
		exponent -= 8;						\
	}								\
	if(Dbl_iszero_hiddenhigh3mantissa(dbl_opndA)) {			\
		Dbl_leftshiftby4(dbl_opndA,dbl_opndB);			\
		exponent -= 4;						\
	}								\
	while(Dbl_iszero_hidden(dbl_opndA)) {				\
		Dbl_leftshiftby1(dbl_opndA,dbl_opndB);			\
		exponent -= 1;						\
	}

#define Twoword_add(src1dstA,src1dstB,src2A,src2B)		\
	/*							\
	 * want this macro to generate:				\
	 *	ADD	src1dstB,src2B,src1dstB;		\
	 *	ADDC	src1dstA,src2A,src1dstA;		\
	 */							\
	if ((src1dstB) + (src2B) < (src1dstB)) Dallp1(src1dstA)++; \
	Dallp1(src1dstA) += (src2A);				\
	Dallp2(src1dstB) += (src2B)

#define Twoword_subtract(src1dstA,src1dstB,src2A,src2B)		\
	/*							\
	 * want this macro to generate:				\
	 *	SUB	src1dstB,src2B,src1dstB;		\
	 *	SUBB	src1dstA,src2A,src1dstA;		\
	 */							\
	if ((src1dstB) < (src2B)) Dallp1(src1dstA)--;		\
	Dallp1(src1dstA) -= (src2A);				\
	Dallp2(src1dstB) -= (src2B)

#define Dbl_setoverflow(resultA,resultB)				\
	/* set result to infinity or largest number */			\
	switch (Rounding_mode()) {					\
		case ROUNDPLUS:						\
			if (Dbl_isone_sign(resultA)) {			\
				Dbl_setlargestnegative(resultA,resultB); \
			}						\
			else {						\
				Dbl_setinfinitypositive(resultA,resultB); \
			}						\
			break;						\
		case ROUNDMINUS:					\
			if (Dbl_iszero_sign(resultA)) {			\
				Dbl_setlargestpositive(resultA,resultB); \
			}						\
			else {						\
				Dbl_setinfinitynegative(resultA,resultB); \
			}						\
			break;						\
		case ROUNDNEAREST:					\
			Dbl_setinfinity_exponentmantissa(resultA,resultB); \
			break;						\
		case ROUNDZERO:						\
			Dbl_setlargest_exponentmantissa(resultA,resultB); \
	}

#define Dbl_denormalize(opndp1,opndp2,exponent,guard,sticky,inexact)	\
    Dbl_clear_signexponent_set_hidden(opndp1);				\
    if (exponent >= (1-DBL_P)) {					\
	if (exponent >= -31) {						\
	    guard = (Dallp2(opndp2) >> (-(exponent))) & 1;		\
	    if (exponent < 0) sticky |= Dallp2(opndp2) << (32+exponent); \
	    if (exponent > -31) {					\
		Variable_shift_double(opndp1,opndp2,1-exponent,opndp2);	\
		Dallp1(opndp1) >>= 1-exponent;				\
	    }								\
	    else {							\
		Dallp2(opndp2) = Dallp1(opndp1);			\
		Dbl_setzerop1(opndp1);					\
	    }								\
	}								\
	else {								\
	    guard = (Dallp1(opndp1) >> (-32-(exponent))) & 1;		\
	    if (exponent == -32) sticky |= Dallp2(opndp2);		\
	    else sticky |= (Dallp2(opndp2) | Dallp1(opndp1) << (64+(exponent))); \
	    Dallp2(opndp2) = Dallp1(opndp1) >> (-31-(exponent));	\
	    Dbl_setzerop1(opndp1);					\
	}								\
	inexact = guard | sticky;					\
    }									\
    else {								\
	guard = 0;							\
	sticky |= (Dallp1(opndp1) | Dallp2(opndp2));			\
	Dbl_setzero(opndp1,opndp2);					\
	inexact = sticky;						\
    }

int dbl_fadd(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);
int dbl_fcmp(dbl_floating_point *, dbl_floating_point *, unsigned int, unsigned int *);
int dbl_fdiv(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);
int dbl_fmpy(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);
int dbl_frem(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);
int dbl_fsqrt(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);
int dbl_fsub(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);

int sgl_to_dbl_fcnvff(sgl_floating_point *, sgl_floating_point *, dbl_floating_point *, unsigned int *);
int dbl_to_sgl_fcnvff(dbl_floating_point *, dbl_floating_point *, sgl_floating_point *, unsigned int *);

int dbl_frnd(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);

license: unknown


code:
/*	$OpenBSD: dfadd.c,v 1.7 2002/11/29 09:27:34 deraadt Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)dfadd.c: Revision: 2.7.88.1 Date: 93/12/07 15:05:35 */

#include "float.h"
#include "dbl_float.h"

/*
 * Double_add: add two double precision values.
 */
int
dbl_fadd(leftptr, rightptr, dstptr, status)
    dbl_floating_point *leftptr, *rightptr, *dstptr;
    unsigned int *status;
{
    register unsigned int signless_upper_left, signless_upper_right, save;
    register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
    register unsigned int resultp1 = 0, resultp2 = 0;

    register int result_exponent, right_exponent, diff_exponent;
    register int sign_save, jumpsize;
    register int inexact = FALSE;
    register int underflowtrap;

    /* Create local copies of the numbers */
    Dbl_copyfromptr(leftptr,leftp1,leftp2);
    Dbl_copyfromptr(rightptr,rightp1,rightp2);

    /* A zero "save" helps discover equal operands (for later),	*
     * and is used in swapping operands (if needed).		*/
    Dbl_xortointp1(leftp1,rightp1,/*to*/save);

    /*
     * check first operand for NaN's or infinity
     */
    if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT)
	{
	if (Dbl_iszero_mantissa(leftp1,leftp2))
	    {
	    if (Dbl_isnotnan(rightp1,rightp2))
		{
		if (Dbl_isinfinity(rightp1,rightp2) && save!=0)
		    {
		    /*
		     * invalid since operands are opposite signed infinity's
		     */
		    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		    Set_invalidflag();
		    Dbl_makequietnan(resultp1,resultp2);
		    Dbl_copytoptr(resultp1,resultp2,dstptr);
		    return(NOEXCEPTION);
		    }
		/*
		 * return infinity
		 */
		Dbl_copytoptr(leftp1,leftp2,dstptr);
		return(NOEXCEPTION);
		}
	    }
	else
	    {
	    /*
	     * is NaN; signaling or quiet?
	     */
	    if (Dbl_isone_signaling(leftp1))
		{
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Dbl_set_quiet(leftp1);
	    }
	    /*
	     * is second operand a signaling NaN?
	     */
	    else if (Dbl_is_signalingnan(rightp1))
		{
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Dbl_set_quiet(rightp1);
		Dbl_copytoptr(rightp1,rightp2,dstptr);
		return(NOEXCEPTION);
		}
	    /*
	     * return quiet NaN
	     */
	    Dbl_copytoptr(leftp1,leftp2,dstptr);
	    return(NOEXCEPTION);
	    }
	} /* End left NaN or Infinity processing */
    /*
     * check second operand for NaN's or infinity
     */
    if (Dbl_isinfinity_exponent(rightp1))
	{
	if (Dbl_iszero_mantissa(rightp1,rightp2))
	    {
	    /* return infinity */
	    Dbl_copytoptr(rightp1,rightp2,dstptr);
	    return(NOEXCEPTION);
	    }
	/*
	 * is NaN; signaling or quiet?
	 */
	if (Dbl_isone_signaling(rightp1))
	    {
	    /* trap if INVALIDTRAP enabled */
	    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
	    /* make NaN quiet */
	    Set_invalidflag();
	    Dbl_set_quiet(rightp1);
	    }
	/*
	 * return quiet NaN
	 */
	Dbl_copytoptr(rightp1,rightp2,dstptr);
	return(NOEXCEPTION);
	} /* End right NaN or Infinity processing */

    /* Invariant: Must be dealing with finite numbers */

    /* Compare operands by removing the sign */
    Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
    Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);

    /* sign difference selects add or sub operation. */
    if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right))
	{
	/* Set the left operand to the larger one by XOR swap	*
	 *  First finish the first word using "save"		*/
	Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
	Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
	Dbl_swap_lower(leftp2,rightp2);
	result_exponent = Dbl_exponent(leftp1);
	}
    /* Invariant:  left is not smaller than right. */

    if((right_exponent = Dbl_exponent(rightp1)) == 0)
	{
	/* Denormalized operands.  First look for zeroes */
	if(Dbl_iszero_mantissa(rightp1,rightp2))
	    {
	    /* right is zero */
	    if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
		{
		/* Both operands are zeros */
		if(Is_rounding_mode(ROUNDMINUS))
		    {
		    Dbl_or_signs(leftp1,/*with*/rightp1);
		    }
		else
		    {
		    Dbl_and_signs(leftp1,/*with*/rightp1);
		    }
		}
	    else
		{
		/* Left is not a zero and must be the result.  Trapped
		 * underflows are signaled if left is denormalized.  Result
		 * is always exact. */
		if( (result_exponent == 0) && Is_underflowtrap_enabled() )
		    {
		    /* need to normalize results mantissa */
		    sign_save = Dbl_signextendedsign(leftp1);
		    Dbl_leftshiftby1(leftp1,leftp2);
		    Dbl_normalize(leftp1,leftp2,result_exponent);
		    Dbl_set_sign(leftp1,/*using*/sign_save);
		    Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
		    Dbl_copytoptr(leftp1,leftp2,dstptr);
		    /* inexact = FALSE */
		    return(UNDERFLOWEXCEPTION);
		    }
		}
	    Dbl_copytoptr(leftp1,leftp2,dstptr);
	    return(NOEXCEPTION);
	    }

	/* Neither are zeroes */
	Dbl_clear_sign(rightp1);	/* Exponent is already cleared */
	if(result_exponent == 0 )
	    {
	    /* Both operands are denormalized.  The result must be exact
	     * and is simply calculated.  A sum could become normalized and a
	     * difference could cancel to a true zero. */
	    if( (/*signed*/int) save < 0 )
		{
		Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
		/*into*/resultp1,resultp2);
		if(Dbl_iszero_mantissa(resultp1,resultp2))
		    {
		    if(Is_rounding_mode(ROUNDMINUS))
			{
			Dbl_setone_sign(resultp1);
			}
		    else
			{
			Dbl_setzero_sign(resultp1);
			}
		    Dbl_copytoptr(resultp1,resultp2,dstptr);
		    return(NOEXCEPTION);
		    }
		}
	    else
		{
		Dbl_addition(leftp1,leftp2,rightp1,rightp2,
		/*into*/resultp1,resultp2);
		if(Dbl_isone_hidden(resultp1))
		    {
		    Dbl_copytoptr(resultp1,resultp2,dstptr);
		    return(NOEXCEPTION);
		    }
		}
	    if(Is_underflowtrap_enabled())
		{
		/* need to normalize result */
		sign_save = Dbl_signextendedsign(resultp1);
		Dbl_leftshiftby1(resultp1,resultp2);
		Dbl_normalize(resultp1,resultp2,result_exponent);
		Dbl_set_sign(resultp1,/*using*/sign_save);
		Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		/* inexact = FALSE */
		return(UNDERFLOWEXCEPTION);
		}
	    Dbl_copytoptr(resultp1,resultp2,dstptr);
	    return(NOEXCEPTION);
	    }
	right_exponent = 1;	/* Set exponent to reflect different bias
				 * with denomalized numbers. */
	}
    else
	{
	Dbl_clear_signexponent_set_hidden(rightp1);
	}
    Dbl_clear_exponent_set_hidden(leftp1);
    diff_exponent = result_exponent - right_exponent;

    /*
     * Special case alignment of operands that would force alignment
     * beyond the extent of the extension.  A further optimization
     * could special case this but only reduces the path length for this
     * infrequent case.
     */
    if(diff_exponent > DBL_THRESHOLD)
	{
	diff_exponent = DBL_THRESHOLD;
	}

    /* Align right operand by shifting to right */
    Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent,
    /*and lower to*/extent);

    /* Treat sum and difference of the operands separately. */
    if( (/*signed*/int) save < 0 )
	{
	/*
	 * Difference of the two operands.  Their can be no overflow.  A
	 * borrow can occur out of the hidden bit and force a post
	 * normalization phase.
	 */
	Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2,
	/*with*/extent,/*into*/resultp1,resultp2);
	if(Dbl_iszero_hidden(resultp1))
	    {
	    /* Handle normalization */
	    /* A straight forward algorithm would now shift the result
	     * and extension left until the hidden bit becomes one.  Not
	     * all of the extension bits need participate in the shift.
	     * Only the two most significant bits (round and guard) are
	     * needed.  If only a single shift is needed then the guard
	     * bit becomes a significant low order bit and the extension
	     * must participate in the rounding.  If more than a single
	     * shift is needed, then all bits to the right of the guard
	     * bit are zeros, and the guard bit may or may not be zero. */
	    sign_save = Dbl_signextendedsign(resultp1);
	    Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2);

	    /* Need to check for a zero result.  The sign and exponent
	     * fields have already been zeroed.  The more efficient test
	     * of the full object can be used.
	     */
	    if(Dbl_iszero(resultp1,resultp2))
		/* Must have been "x-x" or "x+(-x)". */
		{
		if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		return(NOEXCEPTION);
		}
	    result_exponent--;
	    /* Look to see if normalization is finished. */
	    if(Dbl_isone_hidden(resultp1))
		{
		if(result_exponent==0)
		    {
		    /* Denormalized, exponent should be zero.  Left operand *
		     * was normalized, so extent (guard, round) was zero    */
		    goto underflow;
		    }
		else
		    {
		    /* No further normalization is needed. */
		    Dbl_set_sign(resultp1,/*using*/sign_save);
		    Ext_leftshiftby1(extent);
		    goto round;
		    }
		}

	    /* Check for denormalized, exponent should be zero.  Left    *
	     * operand was normalized, so extent (guard, round) was zero */
	    if(!(underflowtrap = Is_underflowtrap_enabled()) &&
	       result_exponent==0) goto underflow;

	    /* Shift extension to complete one bit of normalization and
	     * update exponent. */
	    Ext_leftshiftby1(extent);

	    /* Discover first one bit to determine shift amount.  Use a
	     * modified binary search.  We have already shifted the result
	     * one position right and still not found a one so the remainder
	     * of the extension must be zero and simplifies rounding. */
	    /* Scan bytes */
	    while(Dbl_iszero_hiddenhigh7mantissa(resultp1))
		{
		Dbl_leftshiftby8(resultp1,resultp2);
		if((result_exponent -= 8) <= 0  && !underflowtrap)
		    goto underflow;
		}
	    /* Now narrow it down to the nibble */
	    if(Dbl_iszero_hiddenhigh3mantissa(resultp1))
		{
		/* The lower nibble contains the normalizing one */
		Dbl_leftshiftby4(resultp1,resultp2);
		if((result_exponent -= 4) <= 0 && !underflowtrap)
		    goto underflow;
		}
	    /* Select case were first bit is set (already normalized)
	     * otherwise select the proper shift. */
	    if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7)
		{
		/* Already normalized */
		if(result_exponent <= 0) goto underflow;
		Dbl_set_sign(resultp1,/*using*/sign_save);
		Dbl_set_exponent(resultp1,/*using*/result_exponent);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		return(NOEXCEPTION);
		}
	    Dbl_sethigh4bits(resultp1,/*using*/sign_save);
	    switch(jumpsize)
		{
		case 1:
		    {
		    Dbl_leftshiftby3(resultp1,resultp2);
		    result_exponent -= 3;
		    break;
		    }
		case 2:
		case 3:
		    {
		    Dbl_leftshiftby2(resultp1,resultp2);
		    result_exponent -= 2;
		    break;
		    }
		case 4:
		case 5:
		case 6:
		case 7:
		    {
		    Dbl_leftshiftby1(resultp1,resultp2);
		    result_exponent -= 1;
		    break;
		    }
		}
	    if(result_exponent > 0)
		{
		Dbl_set_exponent(resultp1,/*using*/result_exponent);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		return(NOEXCEPTION);	/* Sign bit is already set */
		}
	    /* Fixup potential underflows */
	  underflow:
	    if(Is_underflowtrap_enabled())
		{
		Dbl_set_sign(resultp1,sign_save);
		Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		/* inexact = FALSE */
		return(UNDERFLOWEXCEPTION);
		}
	    /*
	     * Since we cannot get an inexact denormalized result,
	     * we can now return.
	     */
	    Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent);
	    Dbl_clear_signexponent(resultp1);
	    Dbl_set_sign(resultp1,sign_save);
	    Dbl_copytoptr(resultp1,resultp2,dstptr);
	    return(NOEXCEPTION);
	    } /* end if(hidden...)... */
	/* Fall through and round */
	} /* end if(save < 0)... */
    else
	{
	/* Add magnitudes */
	Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2);
	if(Dbl_isone_hiddenoverflow(resultp1))
	    {
	    /* Prenormalization required. */
	    Dbl_rightshiftby1_withextent(resultp2,extent,extent);
	    Dbl_arithrightshiftby1(resultp1,resultp2);
	    result_exponent++;
	    } /* end if hiddenoverflow... */
	} /* end else ...add magnitudes... */

    /* Round the result.  If the extension is all zeros,then the result is
     * exact.  Otherwise round in the correct direction.  No underflow is
     * possible. If a postnormalization is necessary, then the mantissa is
     * all zeros so no shift is needed. */
  round:
    if(Ext_isnotzero(extent))
	{
	inexact = TRUE;
	switch(Rounding_mode())
	    {
	    case ROUNDNEAREST: /* The default. */
	    if(Ext_isone_sign(extent))
		{
		/* at least 1/2 ulp */
		if(Ext_isnotzero_lower(extent)  ||
		  Dbl_isone_lowmantissap2(resultp2))
		    {
		    /* either exactly half way and odd or more than 1/2ulp */
		    Dbl_increment(resultp1,resultp2);
		    }
		}
	    break;

	    case ROUNDPLUS:
	    if(Dbl_iszero_sign(resultp1))
		{
		/* Round up positive results */
		Dbl_increment(resultp1,resultp2);
		}
	    break;

	    case ROUNDMINUS:
	    if(Dbl_isone_sign(resultp1))
		{
		/* Round down negative results */
		Dbl_increment(resultp1,resultp2);
		}

	    case ROUNDZERO:;
	    /* truncate is simple */
	    } /* end switch... */
	if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
	}
    if(result_exponent == DBL_INFINITY_EXPONENT)
	{
	/* Overflow */
	if(Is_overflowtrap_enabled())
	    {
	    Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
	    Dbl_copytoptr(resultp1,resultp2,dstptr);
	    if (inexact) {
		if (Is_inexacttrap_enabled())
			return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
		else
			Set_inexactflag();
	    }
	    return(OVERFLOWEXCEPTION);
	    }
	else
	    {
	    inexact = TRUE;
	    Set_overflowflag();
	    Dbl_setoverflow(resultp1,resultp2);
	    }
	}
    else Dbl_set_exponent(resultp1,result_exponent);
    Dbl_copytoptr(resultp1,resultp2,dstptr);
    if(inexact) {
	if(Is_inexacttrap_enabled())
	    return(INEXACTEXCEPTION);
	else
	    Set_inexactflag();
    }
    return(NOEXCEPTION);
}


code:
/*	$OpenBSD: dfcmp.c,v 1.6 2002/09/20 19:26:59 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)dfcmp.c: Revision: 1.7.88.2 Date: 93/12/08 13:26:38 */

#include "float.h"
#include "dbl_float.h"

/*
 * dbl_cmp: compare two values
 */
int
dbl_fcmp(leftptr, rightptr, cond, status)
    dbl_floating_point *leftptr, *rightptr;
    unsigned int cond; /* The predicate to be tested */
    unsigned int *status;
{
    register unsigned int leftp1, leftp2, rightp1, rightp2;
    register int xorresult;

    /* Create local copies of the numbers */
    Dbl_copyfromptr(leftptr,leftp1,leftp2);
    Dbl_copyfromptr(rightptr,rightp1,rightp2);
    /*
     * Test for NaN
     */
    if(    (Dbl_exponent(leftp1) == DBL_INFINITY_EXPONENT)
	|| (Dbl_exponent(rightp1) == DBL_INFINITY_EXPONENT) )
	{
	/* Check if a NaN is involved.  Signal an invalid exception when
	 * comparing a signaling NaN or when comparing quiet NaNs and the
	 * low bit of the condition is set */
	if( ((Dbl_exponent(leftp1) == DBL_INFINITY_EXPONENT)
	    && Dbl_isnotzero_mantissa(leftp1,leftp2)
	    && (Exception(cond) || Dbl_isone_signaling(leftp1)))
	   ||
	    ((Dbl_exponent(rightp1) == DBL_INFINITY_EXPONENT)
	    && Dbl_isnotzero_mantissa(rightp1,rightp2)
	    && (Exception(cond) || Dbl_isone_signaling(rightp1))) )
	    {
	    if( Is_invalidtrap_enabled() ) {
		Set_status_cbit(Unordered(cond));
		return(INVALIDEXCEPTION);
	    }
	    else Set_invalidflag();
	    Set_status_cbit(Unordered(cond));
	    return(NOEXCEPTION);
	    }
	/* All the exceptional conditions are handled, now special case
	   NaN compares */
	else if( ((Dbl_exponent(leftp1) == DBL_INFINITY_EXPONENT)
	    && Dbl_isnotzero_mantissa(leftp1,leftp2))
	   ||
	    ((Dbl_exponent(rightp1) == DBL_INFINITY_EXPONENT)
	    && Dbl_isnotzero_mantissa(rightp1,rightp2)) )
	    {
	    /* NaNs always compare unordered. */
	    Set_status_cbit(Unordered(cond));
	    return(NOEXCEPTION);
	    }
	/* infinities will drop down to the normal compare mechanisms */
	}
    /* First compare for unequal signs => less or greater or
     * special equal case */
    Dbl_xortointp1(leftp1,rightp1,xorresult);
    if( xorresult < 0 )
	{
	/* left negative => less, left positive => greater.
	 * equal is possible if both operands are zeros. */
	if( Dbl_iszero_exponentmantissa(leftp1,leftp2)
	  && Dbl_iszero_exponentmantissa(rightp1,rightp2) )
	    {
	    Set_status_cbit(Equal(cond));
	    }
	else if( Dbl_isone_sign(leftp1) )
	    {
	    Set_status_cbit(Lessthan(cond));
	    }
	else
	    {
	    Set_status_cbit(Greaterthan(cond));
	    }
	}
    /* Signs are the same.  Treat negative numbers separately
     * from the positives because of the reversed sense.  */
    else if(Dbl_isequal(leftp1,leftp2,rightp1,rightp2))
	{
	Set_status_cbit(Equal(cond));
	}
    else if( Dbl_iszero_sign(leftp1) )
	{
	/* Positive compare */
	if( Dbl_allp1(leftp1) < Dbl_allp1(rightp1) )
	    {
	    Set_status_cbit(Lessthan(cond));
	    }
	else if( Dbl_allp1(leftp1) > Dbl_allp1(rightp1) )
	    {
	    Set_status_cbit(Greaterthan(cond));
	    }
	else
	    {
	    /* Equal first parts.  Now we must use unsigned compares to
	     * resolve the two possibilities. */
	    if( Dbl_allp2(leftp2) < Dbl_allp2(rightp2) )
		{
		Set_status_cbit(Lessthan(cond));
		}
	    else
		{
		Set_status_cbit(Greaterthan(cond));
		}
	    }
	}
    else
	{
	/* Negative compare.  Signed or unsigned compares
	 * both work the same.  That distinction is only
	 * important when the sign bits differ. */
	if( Dbl_allp1(leftp1) > Dbl_allp1(rightp1) )
	    {
	    Set_status_cbit(Lessthan(cond));
	    }
	else if( Dbl_allp1(leftp1) < Dbl_allp1(rightp1) )
	    {
	    Set_status_cbit(Greaterthan(cond));
	    }
	else
	    {
	    /* Equal first parts.  Now we must use unsigned compares to
	     * resolve the two possibilities. */
	    if( Dbl_allp2(leftp2) > Dbl_allp2(rightp2) )
		{
		Set_status_cbit(Lessthan(cond));
		}
	    else
		{
		Set_status_cbit(Greaterthan(cond));
		}
	    }
	}
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: dfdiv.c,v 1.5 2002/05/07 22:19:30 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)dfdiv.c: Revision: 2.11.88.1 Date: 93/12/07 15:05:39 */

#include "float.h"
#include "dbl_float.h"

/*
 *  Double Precision Floating-point Divide
 */

int
dbl_fdiv(srcptr1,srcptr2,dstptr,status)
	dbl_floating_point *srcptr1, *srcptr2, *dstptr;
	unsigned int *status;
{
	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
	register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
	register int dest_exponent, count;
	register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
	int is_tiny;

	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
	/*
	 * set sign bit of result
	 */
	if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
		Dbl_setnegativezerop1(resultp1);
	else Dbl_setzerop1(resultp1);
	/*
	 * check first operand for NaN's or infinity
	 */
	if (Dbl_isinfinity_exponent(opnd1p1)) {
		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
				if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
					/*
					 * invalid since both operands
					 * are infinity
					 */
					if (Is_invalidtrap_enabled())
						return(INVALIDEXCEPTION);
					Set_invalidflag();
					Dbl_makequietnan(resultp1,resultp2);
					Dbl_copytoptr(resultp1,resultp2,dstptr);
					return(NOEXCEPTION);
				}
				/*
				 * return infinity
				 */
				Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
				Dbl_copytoptr(resultp1,resultp2,dstptr);
				return(NOEXCEPTION);
			}
		}
		else {
			/*
			 * is NaN; signaling or quiet?
			 */
			if (Dbl_isone_signaling(opnd1p1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Dbl_set_quiet(opnd1p1);
			}
			/*
			 * is second operand a signaling NaN?
			 */
			else if (Dbl_is_signalingnan(opnd2p1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Dbl_set_quiet(opnd2p1);
				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
				return(NOEXCEPTION);
			}
			/*
			 * return quiet NaN
			 */
			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
			return(NOEXCEPTION);
		}
	}
	/*
	 * check second operand for NaN's or infinity
	 */
	if (Dbl_isinfinity_exponent(opnd2p1)) {
		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
			/*
			 * return zero
			 */
			Dbl_setzero_exponentmantissa(resultp1,resultp2);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(NOEXCEPTION);
		}
		/*
		 * is NaN; signaling or quiet?
		 */
		if (Dbl_isone_signaling(opnd2p1)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			/* make NaN quiet */
			Set_invalidflag();
			Dbl_set_quiet(opnd2p1);
		}
		/*
		 * return quiet NaN
		 */
		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
		return(NOEXCEPTION);
	}
	/*
	 * check for division by zero
	 */
	if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
		if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
			/* invalid since both operands are zero */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			Set_invalidflag();
			Dbl_makequietnan(resultp1,resultp2);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(NOEXCEPTION);
		}
		if (Is_divisionbyzerotrap_enabled())
			return(DIVISIONBYZEROEXCEPTION);
		Set_divisionbyzeroflag();
		Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		return(NOEXCEPTION);
	}
	/*
	 * Generate exponent
	 */
	dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;

	/*
	 * Generate mantissa
	 */
	if (Dbl_isnotzero_exponent(opnd1p1)) {
		/* set hidden bit */
		Dbl_clear_signexponent_set_hidden(opnd1p1);
	}
	else {
		/* check for zero */
		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
			Dbl_setzero_exponentmantissa(resultp1,resultp2);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(NOEXCEPTION);
		}
		/* is denormalized, want to normalize */
		Dbl_clear_signexponent(opnd1p1);
		Dbl_leftshiftby1(opnd1p1,opnd1p2);
		Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
	}
	/* opnd2 needs to have hidden bit set with msb in hidden bit */
	if (Dbl_isnotzero_exponent(opnd2p1)) {
		Dbl_clear_signexponent_set_hidden(opnd2p1);
	}
	else {
		/* is denormalized; want to normalize */
		Dbl_clear_signexponent(opnd2p1);
		Dbl_leftshiftby1(opnd2p1,opnd2p2);
		while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
			dest_exponent+=8;
			Dbl_leftshiftby8(opnd2p1,opnd2p2);
		}
		if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
			dest_exponent+=4;
			Dbl_leftshiftby4(opnd2p1,opnd2p2);
		}
		while (Dbl_iszero_hidden(opnd2p1)) {
			dest_exponent++;
			Dbl_leftshiftby1(opnd2p1,opnd2p2);
		}
	}

	/* Divide the source mantissas */

	/*
	 * A non-restoring divide algorithm is used.
	 */
	Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
	Dbl_setzero(opnd3p1,opnd3p2);
	for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
		Dbl_leftshiftby1(opnd1p1,opnd1p2);
		Dbl_leftshiftby1(opnd3p1,opnd3p2);
		if (Dbl_iszero_sign(opnd1p1)) {
			Dbl_setone_lowmantissap2(opnd3p2);
			Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
		}
		else {
			Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
		}
	}
	if (count <= DBL_P) {
		Dbl_leftshiftby1(opnd3p1,opnd3p2);
		Dbl_setone_lowmantissap2(opnd3p2);
		Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
		if (Dbl_iszero_hidden(opnd3p1)) {
			Dbl_leftshiftby1(opnd3p1,opnd3p2);
			dest_exponent--;
		}
	}
	else {
		if (Dbl_iszero_hidden(opnd3p1)) {
			/* need to get one more bit of result */
			Dbl_leftshiftby1(opnd1p1,opnd1p2);
			Dbl_leftshiftby1(opnd3p1,opnd3p2);
			if (Dbl_iszero_sign(opnd1p1)) {
				Dbl_setone_lowmantissap2(opnd3p2);
				Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
			}
			else {
				Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
			}
			dest_exponent--;
		}
		if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE;
		stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
	}
	inexact = guardbit | stickybit;

	/*
	 * round result
	 */
	if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
		Dbl_clear_signexponent(opnd3p1);
		switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Dbl_iszero_sign(resultp1))
					Dbl_increment(opnd3p1,opnd3p2);
				break;
			case ROUNDMINUS:
				if (Dbl_isone_sign(resultp1))
					Dbl_increment(opnd3p1,opnd3p2);
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Dbl_isone_lowmantissap2(opnd3p2))) {
					Dbl_increment(opnd3p1,opnd3p2);
				}
		}
		if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
	}
	Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);

	/*
	 * Test for overflow
	 */
	if (dest_exponent >= DBL_INFINITY_EXPONENT) {
		/* trap if OVERFLOWTRAP enabled */
		if (Is_overflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			if (inexact) {
				if (Is_inexacttrap_enabled())
					return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
				else
					Set_inexactflag();
			}
			return(OVERFLOWEXCEPTION);
		}
		Set_overflowflag();
		/* set result to infinity or largest number */
		Dbl_setoverflow(resultp1,resultp2);
		inexact = TRUE;
	}
	/*
	 * Test for underflow
	 */
	else if (dest_exponent <= 0) {
		/* trap if UNDERFLOWTRAP enabled */
		if (Is_underflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			if (inexact) {
				if (Is_inexacttrap_enabled())
					return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
				else
					Set_inexactflag();
			}
			return(UNDERFLOWEXCEPTION);
		}

		/* Determine if should set underflow flag */
		is_tiny = TRUE;
		if (dest_exponent == 0 && inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Dbl_iszero_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
					if (Dbl_isone_hiddenoverflow(opnd3p1))
						is_tiny = FALSE;
					Dbl_decrement(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDMINUS:
				if (Dbl_isone_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
					if (Dbl_isone_hiddenoverflow(opnd3p1))
						is_tiny = FALSE;
					Dbl_decrement(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Dbl_isone_lowmantissap2(opnd3p2))) {
					Dbl_increment(opnd3p1,opnd3p2);
					if (Dbl_isone_hiddenoverflow(opnd3p1))
						is_tiny = FALSE;
					Dbl_decrement(opnd3p1,opnd3p2);
				}
				break;
			}
		}

		/*
		 * denormalize result or set to signed zero
		 */
		stickybit = inexact;
		Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
		 stickybit,inexact);

		/* return rounded number */
		if (inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Dbl_iszero_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDMINUS:
				if (Dbl_isone_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Dbl_isone_lowmantissap2(opnd3p2))) {
					Dbl_increment(opnd3p1,opnd3p2);
				}
				break;
			}
			if (is_tiny)
				Set_underflowflag();
		}
		Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
	}
	else Dbl_set_exponent(resultp1,dest_exponent);
	Dbl_copytoptr(resultp1,resultp2,dstptr);

	/* check for inexact */
	if (inexact) {
		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
		else Set_inexactflag();
	}
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: dfmpy.c,v 1.6 2015/05/07 01:55:43 jsg Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)dfmpy.c: Revision: 2.11.88.1 Date: 93/12/07 15:05:41 */

#include "float.h"
#include "dbl_float.h"

/*
 *  Double Precision Floating-point Multiply
 */

int
dbl_fmpy(srcptr1,srcptr2,dstptr,status)
	dbl_floating_point *srcptr1, *srcptr2, *dstptr;
	unsigned int *status;
{
	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
	register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
	register int dest_exponent, count;
	register int inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
	int is_tiny;

	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);

	/*
	 * set sign bit of result
	 */
	if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
		Dbl_setnegativezerop1(resultp1);
	else Dbl_setzerop1(resultp1);
	/*
	 * check first operand for NaN's or infinity
	 */
	if (Dbl_isinfinity_exponent(opnd1p1)) {
		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
				if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
					/*
					 * invalid since operands are infinity
					 * and zero
					 */
					if (Is_invalidtrap_enabled())
						return(INVALIDEXCEPTION);
					Set_invalidflag();
					Dbl_makequietnan(resultp1,resultp2);
					Dbl_copytoptr(resultp1,resultp2,dstptr);
					return(NOEXCEPTION);
				}
				/*
				 * return infinity
				 */
				Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
				Dbl_copytoptr(resultp1,resultp2,dstptr);
				return(NOEXCEPTION);
			}
		}
		else {
			/*
			 * is NaN; signaling or quiet?
			 */
			if (Dbl_isone_signaling(opnd1p1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Dbl_set_quiet(opnd1p1);
			}
			/*
			 * is second operand a signaling NaN?
			 */
			else if (Dbl_is_signalingnan(opnd2p1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Dbl_set_quiet(opnd2p1);
				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
				return(NOEXCEPTION);
			}
			/*
			 * return quiet NaN
			 */
			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
			return(NOEXCEPTION);
		}
	}
	/*
	 * check second operand for NaN's or infinity
	 */
	if (Dbl_isinfinity_exponent(opnd2p1)) {
		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
			if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
				/* invalid since operands are zero & infinity */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				Set_invalidflag();
				Dbl_makequietnan(opnd2p1,opnd2p2);
				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
				return(NOEXCEPTION);
			}
			/*
			 * return infinity
			 */
			Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(NOEXCEPTION);
		}
		/*
		 * is NaN; signaling or quiet?
		 */
		if (Dbl_isone_signaling(opnd2p1)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			/* make NaN quiet */
			Set_invalidflag();
			Dbl_set_quiet(opnd2p1);
		}
		/*
		 * return quiet NaN
		 */
		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
		return(NOEXCEPTION);
	}
	/*
	 * Generate exponent
	 */
	dest_exponent = Dbl_exponent(opnd1p1) + Dbl_exponent(opnd2p1) -DBL_BIAS;

	/*
	 * Generate mantissa
	 */
	if (Dbl_isnotzero_exponent(opnd1p1)) {
		/* set hidden bit */
		Dbl_clear_signexponent_set_hidden(opnd1p1);
	}
	else {
		/* check for zero */
		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
			Dbl_setzero_exponentmantissa(resultp1,resultp2);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(NOEXCEPTION);
		}
		/* is denormalized, adjust exponent */
		Dbl_clear_signexponent(opnd1p1);
		Dbl_leftshiftby1(opnd1p1,opnd1p2);
		Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
	}
	/* opnd2 needs to have hidden bit set with msb in hidden bit */
	if (Dbl_isnotzero_exponent(opnd2p1)) {
		Dbl_clear_signexponent_set_hidden(opnd2p1);
	}
	else {
		/* check for zero */
		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
			Dbl_setzero_exponentmantissa(resultp1,resultp2);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(NOEXCEPTION);
		}
		/* is denormalized; want to normalize */
		Dbl_clear_signexponent(opnd2p1);
		Dbl_leftshiftby1(opnd2p1,opnd2p2);
		Dbl_normalize(opnd2p1,opnd2p2,dest_exponent);
	}

	/* Multiply two source mantissas together */

	/* make room for guard bits */
	Dbl_leftshiftby7(opnd2p1,opnd2p2);
	Dbl_setzero(opnd3p1,opnd3p2);
	/*
	 * Four bits at a time are inspected in each loop, and a
	 * simple shift and add multiply algorithm is used.
	 */
	for (count=1;count<=DBL_P;count+=4) {
		stickybit |= Dlow4p2(opnd3p2);
		Dbl_rightshiftby4(opnd3p1,opnd3p2);
		if (Dbit28p2(opnd1p2)) {
			/* Twoword_add should be an ADDC followed by an ADD. */
			Twoword_add(opnd3p1, opnd3p2, opnd2p1<<3 | opnd2p2>>29,
				    opnd2p2<<3);
		}
		if (Dbit29p2(opnd1p2)) {
			Twoword_add(opnd3p1, opnd3p2, opnd2p1<<2 | opnd2p2>>30,
				    opnd2p2<<2);
		}
		if (Dbit30p2(opnd1p2)) {
			Twoword_add(opnd3p1, opnd3p2, opnd2p1<<1 | opnd2p2>>31,
				    opnd2p2<<1);
		}
		if (Dbit31p2(opnd1p2)) {
			Twoword_add(opnd3p1, opnd3p2, opnd2p1, opnd2p2);
		}
		Dbl_rightshiftby4(opnd1p1,opnd1p2);
	}
	if (Dbit3p1(opnd3p1)==0) {
		Dbl_leftshiftby1(opnd3p1,opnd3p2);
	}
	else {
		/* result mantissa >= 2. */
		dest_exponent++;
	}
	/* check for denormalized result */
	while (Dbit3p1(opnd3p1)==0) {
		Dbl_leftshiftby1(opnd3p1,opnd3p2);
		dest_exponent--;
	}
	/*
	 * check for guard, sticky and inexact bits
	 */
	stickybit |= Dallp2(opnd3p2) << 25;
	guardbit = (Dallp2(opnd3p2) << 24) >> 31;
	inexact = guardbit | stickybit;

	/* align result mantissa */
	Dbl_rightshiftby8(opnd3p1,opnd3p2);

	/*
	 * round result
	 */
	if (inexact && (dest_exponent>0 || Is_underflowtrap_enabled())) {
		Dbl_clear_signexponent(opnd3p1);
		switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Dbl_iszero_sign(resultp1))
					Dbl_increment(opnd3p1,opnd3p2);
				break;
			case ROUNDMINUS:
				if (Dbl_isone_sign(resultp1))
					Dbl_increment(opnd3p1,opnd3p2);
				break;
			case ROUNDNEAREST:
				if (guardbit &&
				    (stickybit || Dbl_isone_lowmantissap2(opnd3p2)))
					Dbl_increment(opnd3p1,opnd3p2);
				break;
		}
		if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
	}
	Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);

	/*
	 * Test for overflow
	 */
	if (dest_exponent >= DBL_INFINITY_EXPONENT) {
		/* trap if OVERFLOWTRAP enabled */
		if (Is_overflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			if (inexact) {
			    if (Is_inexacttrap_enabled())
				return (OVERFLOWEXCEPTION | INEXACTEXCEPTION);
			    else
				Set_inexactflag();
			}
			return (OVERFLOWEXCEPTION);
		}
		inexact = TRUE;
		Set_overflowflag();
		/* set result to infinity or largest number */
		Dbl_setoverflow(resultp1,resultp2);
	}
	/*
	 * Test for underflow
	 */
	else if (dest_exponent <= 0) {
		/* trap if UNDERFLOWTRAP enabled */
		if (Is_underflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			if (inexact) {
			    if (Is_inexacttrap_enabled())
				return (UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
			    else
				Set_inexactflag();
			}
			return (UNDERFLOWEXCEPTION);
		}

		/* Determine if should set underflow flag */
		is_tiny = TRUE;
		if (dest_exponent == 0 && inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Dbl_iszero_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
					if (Dbl_isone_hiddenoverflow(opnd3p1))
						is_tiny = FALSE;
					Dbl_decrement(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDMINUS:
				if (Dbl_isone_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
					if (Dbl_isone_hiddenoverflow(opnd3p1))
						is_tiny = FALSE;
					Dbl_decrement(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Dbl_isone_lowmantissap2(opnd3p2))) {
					Dbl_increment(opnd3p1,opnd3p2);
					if (Dbl_isone_hiddenoverflow(opnd3p1))
						is_tiny = FALSE;
					Dbl_decrement(opnd3p1,opnd3p2);
				}
				break;
			}
		}

		/*
		 * denormalize result or set to signed zero
		 */
		stickybit = inexact;
		Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
		 stickybit,inexact);

		/* return zero or smallest number */
		if (inexact) {
			switch (Rounding_mode()) {
			case ROUNDPLUS:
				if (Dbl_iszero_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDMINUS:
				if (Dbl_isone_sign(resultp1)) {
					Dbl_increment(opnd3p1,opnd3p2);
				}
				break;
			case ROUNDNEAREST:
				if (guardbit && (stickybit ||
				    Dbl_isone_lowmantissap2(opnd3p2))) {
					Dbl_increment(opnd3p1,opnd3p2);
				}
				break;
			}
			if (is_tiny) Set_underflowflag();
		}
		Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
	}
	else Dbl_set_exponent(resultp1,dest_exponent);
	/* check for inexact */
	Dbl_copytoptr(resultp1,resultp2,dstptr);
	if (inexact) {
		if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
		else Set_inexactflag();
	}
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: dfrem.c,v 1.5 2002/05/07 22:19:30 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)dfrem.c: Revision: 1.7.88.1 Date: 93/12/07 15:05:43 */

#include "float.h"
#include "dbl_float.h"

/*
 *  Double Precision Floating-point Remainder
 */
int
dbl_frem(srcptr1,srcptr2,dstptr,status)
	dbl_floating_point *srcptr1, *srcptr2, *dstptr;
	unsigned int *status;
{
	register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
	register unsigned int resultp1, resultp2;
	register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
	register int roundup = FALSE;

	Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
	Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
	/*
	 * check first operand for NaN's or infinity
	 */
	if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
			if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
				/* invalid since first operand is infinity */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				Set_invalidflag();
				Dbl_makequietnan(resultp1,resultp2);
				Dbl_copytoptr(resultp1,resultp2,dstptr);
				return(NOEXCEPTION);
			}
		}
		else {
			/*
			 * is NaN; signaling or quiet?
			 */
			if (Dbl_isone_signaling(opnd1p1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Dbl_set_quiet(opnd1p1);
			}
			/*
			 * is second operand a signaling NaN?
			 */
			else if (Dbl_is_signalingnan(opnd2p1)) {
				/* trap if INVALIDTRAP enabled */
				if (Is_invalidtrap_enabled())
					return(INVALIDEXCEPTION);
				/* make NaN quiet */
				Set_invalidflag();
				Dbl_set_quiet(opnd2p1);
				Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
				return(NOEXCEPTION);
			}
			/*
			 * return quiet NaN
			 */
			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
			return(NOEXCEPTION);
		}
	}
	/*
	 * check second operand for NaN's or infinity
	 */
	if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
		if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
			/*
			 * return first operand
			 */
			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
			return(NOEXCEPTION);
		}
		/*
		 * is NaN; signaling or quiet?
		 */
		if (Dbl_isone_signaling(opnd2p1)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			/* make NaN quiet */
			Set_invalidflag();
			Dbl_set_quiet(opnd2p1);
		}
		/*
		 * return quiet NaN
		 */
		Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
		return(NOEXCEPTION);
	}
	/*
	 * check second operand for zero
	 */
	if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
		/* invalid since second operand is zero */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		Set_invalidflag();
		Dbl_makequietnan(resultp1,resultp2);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		return(NOEXCEPTION);
	}

	/*
	 * get sign of result
	 */
	resultp1 = opnd1p1;

	/*
	 * check for denormalized operands
	 */
	if (opnd1_exponent == 0) {
		/* check for zero */
		if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
			Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
			return(NOEXCEPTION);
		}
		/* normalize, then continue */
		opnd1_exponent = 1;
		Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
	}
	else {
		Dbl_clear_signexponent_set_hidden(opnd1p1);
	}
	if (opnd2_exponent == 0) {
		/* normalize, then continue */
		opnd2_exponent = 1;
		Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
	}
	else {
		Dbl_clear_signexponent_set_hidden(opnd2p1);
	}

	/* find result exponent and divide step loop count */
	dest_exponent = opnd2_exponent - 1;
	stepcount = opnd1_exponent - opnd2_exponent;

	/*
	 * check for opnd1/opnd2 < 1
	 */
	if (stepcount < 0) {
		/*
		 * check for opnd1/opnd2 > 1/2
		 *
		 * In this case n will round to 1, so
		 *    r = opnd1 - opnd2
		 */
		if (stepcount == -1 &&
		    Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
			/* set sign */
			Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
			/* align opnd2 with opnd1 */
			Dbl_leftshiftby1(opnd2p1,opnd2p2);
			Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
			 opnd2p1,opnd2p2);
			/* now normalize */
			while (Dbl_iszero_hidden(opnd2p1)) {
				Dbl_leftshiftby1(opnd2p1,opnd2p2);
				dest_exponent--;
			}
			Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
			goto testforunderflow;
		}
		/*
		 * opnd1/opnd2 <= 1/2
		 *
		 * In this case n will round to zero, so
		 *    r = opnd1
		 */
		Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
		dest_exponent = opnd1_exponent;
		goto testforunderflow;
	}

	/*
	 * Generate result
	 *
	 * Do iterative subtract until remainder is less than operand 2.
	 */
	while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
		if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
			Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
		}
		Dbl_leftshiftby1(opnd1p1,opnd1p2);
	}
	/*
	 * Do last subtract, then determine which way to round if remainder
	 * is exactly 1/2 of opnd2
	 */
	if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
		Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
		roundup = TRUE;
	}
	if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
		/* division is exact, remainder is zero */
		Dbl_setzero_exponentmantissa(resultp1,resultp2);
		Dbl_copytoptr(resultp1,resultp2,dstptr);
		return(NOEXCEPTION);
	}

	/*
	 * Check for cases where opnd1/opnd2 < n
	 *
	 * In this case the result's sign will be opposite that of
	 * opnd1.  The mantissa also needs some correction.
	 */
	Dbl_leftshiftby1(opnd1p1,opnd1p2);
	if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
		Dbl_invert_sign(resultp1);
		Dbl_leftshiftby1(opnd2p1,opnd2p2);
		Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
	}
	/* check for remainder being exactly 1/2 of opnd2 */
	else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
		Dbl_invert_sign(resultp1);
	}

	/* normalize result's mantissa */
	while (Dbl_iszero_hidden(opnd1p1)) {
		dest_exponent--;
		Dbl_leftshiftby1(opnd1p1,opnd1p2);
	}
	Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);

	/*
	 * Test for underflow
	 */
    testforunderflow:
	if (dest_exponent <= 0) {
		/* trap if UNDERFLOWTRAP enabled */
		if (Is_underflowtrap_enabled()) {
			/*
			 * Adjust bias of result
			 */
			Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
			/* frem is always exact */
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(UNDERFLOWEXCEPTION);
		}
		/*
		 * denormalize result or set to signed zero
		 */
		if (dest_exponent >= (1 - DBL_P)) {
			Dbl_rightshift_exponentmantissa(resultp1,resultp2,
			 1-dest_exponent);
		}
		else {
			Dbl_setzero_exponentmantissa(resultp1,resultp2);
		}
	}
	else Dbl_set_exponent(resultp1,dest_exponent);
	Dbl_copytoptr(resultp1,resultp2,dstptr);
	return(NOEXCEPTION);
}


code:
/*	$OpenBSD: dfsqrt.c,v 1.7 2003/04/10 17:27:58 mickey Exp $	*/
/*
  (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  To anyone who acknowledges that this file is provided "AS IS"
  without any express or implied warranty:
      permission to use, copy, modify, and distribute this file
  for any purpose is hereby granted without fee, provided that
  the above copyright notice and this notice appears in all
  copies, and that the name of Hewlett-Packard Company not be
  used in advertising or publicity pertaining to distribution
  of the software without specific, written prior permission.
  Hewlett-Packard Company makes no representations about the
  suitability of this software for any purpose.
*/
/* @(#)dfsqrt.c: Revision: 1.9.88.1 Date: 93/12/07 15:05:46 */

#include "float.h"
#include "dbl_float.h"

/*
 *  Double Floating-point Square Root
 */

/*ARGSUSED*/
int
dbl_fsqrt(srcptr, null, dstptr, status)
	dbl_floating_point *srcptr, *null, *dstptr;
	unsigned int *status;
{
	register unsigned int srcp1, srcp2, resultp1, resultp2;
	register unsigned int newbitp1, newbitp2, sump1, sump2;
	register int src_exponent;
	register int guardbit = FALSE, even_exponent;

	Dbl_copyfromptr(srcptr,srcp1,srcp2);
	/*
	 * check source operand for NaN or infinity
	 */
	if ((src_exponent = Dbl_exponent(srcp1)) == DBL_INFINITY_EXPONENT) {
		/*
		 * is signaling NaN?
		 */
		if (Dbl_isone_signaling(srcp1)) {
			/* trap if INVALIDTRAP enabled */
			if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
			/* make NaN quiet */
			Set_invalidflag();
			Dbl_set_quiet(srcp1);
		}
		/*
		 * Return quiet NaN or positive infinity.
		 *  Fall thru to negative test if negative infinity.
		 */
		if (Dbl_iszero_sign(srcp1) ||
		    Dbl_isnotzero_mantissa(srcp1,srcp2)) {
			Dbl_copytoptr(srcp1,srcp2,dstptr);
			return(NOEXCEPTION);
		}
	}

	/*
	 * check for zero source operand
	 */
	if (Dbl_iszero_exponentmantissa(srcp1,srcp2)) {
		Dbl_copytoptr(srcp1,srcp2,dstptr);
		return(NOEXCEPTION);
	}

	/*
	 * check for negative source operand
	 */
	if (Dbl_isone_sign(srcp1)) {
		/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Dbl_makequietnan(srcp1,srcp2);
		Dbl_copytoptr(srcp1,srcp2,dstptr);
		return(NOEXCEPTION);
	}

	/*
	 * Generate result
	 */
	if (src_exponent > 0) {
		even_exponent = Dbl_hidden(srcp1);
		Dbl_clear_signexponent_set_hidden(srcp1);
	}
	else {
		/* normalize operand */
		Dbl_clear_signexponent(srcp1);
		src_exponent++;
		Dbl_normalize(srcp1,srcp2,src_exponent);
		even_exponent = src_exponent & 1;
	}
	if (even_exponent) {
		/* exponent is even */
		/* Add comment here.  Explain why odd exponent needs correction */
		Dbl_leftshiftby1(srcp1,srcp2);
	}
	/*
	 * Add comment here.  Explain following algorithm.
	 *
	 * Trust me, it works.
	 *
	 */
	Dbl_setzero(resultp1,resultp2);
	Dbl_allp1(newbitp1) = 1 << (DBL_P - 32);
	Dbl_setzero_mantissap2(newbitp2);
	while (Dbl_isnotzero(newbitp1,newbitp2) && Dbl_isnotzero(srcp1,srcp2)) {
		Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,sump1,sump2);
		if(Dbl_isnotgreaterthan(sump1,sump2,srcp1,srcp2)) {
			Dbl_leftshiftby1(newbitp1,newbitp2);
			/* update result */
			Dbl_addition(resultp1,resultp2,newbitp1,newbitp2,
			 resultp1,resultp2);
			Dbl_subtract(srcp1,srcp2,sump1,sump2,srcp1,srcp2);
			Dbl_rightshiftby2(newbitp1,newbitp2);
		}
		else {
			Dbl_rightshiftby1(newbitp1,newbitp2);
		}
		Dbl_leftshiftby1(srcp1,srcp2);
	}
	/* correct exponent for pre-shift */
	if (even_exponent) {
		Dbl_rightshiftby1(resultp1,resultp2);
	}

	/* check for inexact */
	if (Dbl_isnotzero(srcp1,srcp2)) {
		if (!even_exponent & Dbl_islessthan(resultp1,resultp2,srcp1,srcp2)) {
			Dbl_increment(resultp1,resultp2);
		}
		guardbit = Dbl_lowmantissap2(resultp2);
		Dbl_rightshiftby1(resultp1,resultp2);

		/*  now round result  */
		switch (Rounding_mode()) {
		case ROUNDPLUS:
		     Dbl_increment(resultp1,resultp2);
		     break;
		case ROUNDNEAREST:
		     /* stickybit is always true, so guardbit
		      * is enough to determine rounding */
		     if (guardbit) {
			    Dbl_increment(resultp1,resultp2);
		     }
		     break;
		}
		/* increment result exponent by 1 if mantissa overflowed */
		if (Dbl_isone_hiddenoverflow(resultp1)) src_exponent+=2;

		if (Is_inexacttrap_enabled()) {
			Dbl_set_exponent(resultp1,
			 ((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
			Dbl_copytoptr(resultp1,resultp2,dstptr);
			return(INEXACTEXCEPTION);
		}
		else Set_inexactflag();
	}
	else {
		Dbl_rightshiftby1(resultp1,resultp2);
	}
	Dbl_set_exponent(resultp1,((src_exponent-DBL_BIAS)>>1)+DBL_BIAS);
	Dbl_copytoptr(resultp1,resultp2,dstptr);
	return(NOEXCEPTION);
}