Home | History | Annotate | Line # | Download | only in spmath
dbl_float.h revision 1.1.112.1
      1        1.1  fredette /*	$NetBSD: dbl_float.h,v 1.1.112.1 2008/06/02 13:22:13 mjf Exp $	*/
      2        1.1  fredette 
      3        1.1  fredette /*	$OpenBSD: dbl_float.h,v 1.5 2001/03/29 03:58:17 mickey Exp $	*/
      4        1.1  fredette 
      5        1.1  fredette /*
      6        1.1  fredette  * Copyright 1996 1995 by Open Software Foundation, Inc.
      7        1.1  fredette  *              All Rights Reserved
      8        1.1  fredette  *
      9        1.1  fredette  * Permission to use, copy, modify, and distribute this software and
     10        1.1  fredette  * its documentation for any purpose and without fee is hereby granted,
     11        1.1  fredette  * provided that the above copyright notice appears in all copies and
     12        1.1  fredette  * that both the copyright notice and this permission notice appear in
     13        1.1  fredette  * supporting documentation.
     14        1.1  fredette  *
     15        1.1  fredette  * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE
     16        1.1  fredette  * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
     17        1.1  fredette  * FOR A PARTICULAR PURPOSE.
     18        1.1  fredette  *
     19        1.1  fredette  * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR
     20        1.1  fredette  * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
     21        1.1  fredette  * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT,
     22        1.1  fredette  * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
     23        1.1  fredette  * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
     24        1.1  fredette  */
     25        1.1  fredette /*
     26        1.1  fredette  * pmk1.1
     27        1.1  fredette  */
     28        1.1  fredette /*
     29        1.1  fredette  * (c) Copyright 1986 HEWLETT-PACKARD COMPANY
     30        1.1  fredette  *
     31        1.1  fredette  * To anyone who acknowledges that this file is provided "AS IS"
     32        1.1  fredette  * without any express or implied warranty:
     33        1.1  fredette  *     permission to use, copy, modify, and distribute this file
     34        1.1  fredette  * for any purpose is hereby granted without fee, provided that
     35        1.1  fredette  * the above copyright notice and this notice appears in all
     36        1.1  fredette  * copies, and that the name of Hewlett-Packard Company not be
     37        1.1  fredette  * used in advertising or publicity pertaining to distribution
     38        1.1  fredette  * of the software without specific, written prior permission.
     39        1.1  fredette  * Hewlett-Packard Company makes no representations about the
     40        1.1  fredette  * suitability of this software for any purpose.
     41        1.1  fredette  */
     42        1.1  fredette 
     43        1.1  fredette #include <sys/cdefs.h>
     44        1.1  fredette 
     45        1.1  fredette /**************************************
     46        1.1  fredette  * Declare double precision functions *
     47        1.1  fredette  **************************************/
     48        1.1  fredette 
     49        1.1  fredette /* 32-bit word grabing functions */
     50        1.1  fredette #define Dbl_firstword(value) Dallp1(value)
     51        1.1  fredette #define Dbl_secondword(value) Dallp2(value)
     52        1.1  fredette #define Dbl_thirdword(value) dummy_location
     53        1.1  fredette #define Dbl_fourthword(value) dummy_location
     54        1.1  fredette 
     55        1.1  fredette #define Dbl_sign(object) Dsign(object)
     56        1.1  fredette #define Dbl_exponent(object) Dexponent(object)
     57        1.1  fredette #define Dbl_signexponent(object) Dsignexponent(object)
     58        1.1  fredette #define Dbl_mantissap1(object) Dmantissap1(object)
     59        1.1  fredette #define Dbl_mantissap2(object) Dmantissap2(object)
     60        1.1  fredette #define Dbl_exponentmantissap1(object) Dexponentmantissap1(object)
     61        1.1  fredette #define Dbl_allp1(object) Dallp1(object)
     62        1.1  fredette #define Dbl_allp2(object) Dallp2(object)
     63        1.1  fredette 
     64        1.1  fredette /* dbl_and_signs ands the sign bits of each argument and puts the result
     65        1.1  fredette  * into the first argument. dbl_or_signs ors those same sign bits */
     66        1.1  fredette #define Dbl_and_signs( src1dst, src2)		\
     67        1.1  fredette     Dallp1(src1dst) = (Dallp1(src2)|~(1<<31)) & Dallp1(src1dst)
     68        1.1  fredette #define Dbl_or_signs( src1dst, src2)		\
     69        1.1  fredette     Dallp1(src1dst) = (Dallp1(src2)&(1<<31)) | Dallp1(src1dst)
     70        1.1  fredette 
     71        1.1  fredette /* The hidden bit is always the low bit of the exponent */
     72        1.1  fredette #define Dbl_clear_exponent_set_hidden(srcdst) Deposit_dexponent(srcdst,1)
     73        1.1  fredette #define Dbl_clear_signexponent_set_hidden(srcdst) \
     74        1.1  fredette     Deposit_dsignexponent(srcdst,1)
     75        1.1  fredette #define Dbl_clear_sign(srcdst) Dallp1(srcdst) &= ~(1<<31)
     76        1.1  fredette #define Dbl_clear_signexponent(srcdst) \
     77        1.1  fredette     Dallp1(srcdst) &= Dmantissap1((unsigned)-1)
     78        1.1  fredette 
     79        1.1  fredette /* Exponent field for doubles has already been cleared and may be
     80        1.1  fredette  * included in the shift.  Here we need to generate two double width
     81        1.1  fredette  * variable shifts.  The insignificant bits can be ignored.
     82        1.1  fredette  *      MTSAR f(varamount)
     83        1.1  fredette  *      VSHD	srcdst.high,srcdst.low => srcdst.low
     84        1.1  fredette  *	VSHD	0,srcdst.high => srcdst.high
     85        1.1  fredette  * This is very difficult to model with C expressions since the shift amount
     86        1.1  fredette  * could exceed 32.  */
     87        1.1  fredette /* varamount must be less than 64 */
     88        1.1  fredette #define Dbl_rightshift(srcdstA, srcdstB, varamount)			\
     89        1.1  fredette     {if((varamount) >= 32) {						\
     90        1.1  fredette 	Dallp2(srcdstB) = Dallp1(srcdstA) >> (varamount-32);		\
     91        1.1  fredette 	Dallp1(srcdstA)=0;						\
     92        1.1  fredette     }									\
     93        1.1  fredette     else if(varamount > 0) {						\
     94        1.1  fredette 	Variable_shift_double(Dallp1(srcdstA), Dallp2(srcdstB),		\
     95        1.1  fredette 	  (varamount), Dallp2(srcdstB));				\
     96        1.1  fredette 	Dallp1(srcdstA) >>= varamount;					\
     97        1.1  fredette     } }
     98        1.1  fredette /* varamount must be less than 64 */
     99        1.1  fredette #define Dbl_rightshift_exponentmantissa(srcdstA, srcdstB, varamount)	\
    100        1.1  fredette     {if((varamount) >= 32) {						\
    101        1.1  fredette 	Dallp2(srcdstB) = Dexponentmantissap1(srcdstA) >> ((varamount)-32); \
    102        1.1  fredette 	Dallp1(srcdstA) &= (1<<31);  /* clear exponentmantissa field */ \
    103        1.1  fredette     }									\
    104        1.1  fredette     else if(varamount > 0) {						\
    105        1.1  fredette 	Variable_shift_double(Dexponentmantissap1(srcdstA), Dallp2(srcdstB), \
    106        1.1  fredette 	(varamount), Dallp2(srcdstB));					\
    107        1.1  fredette 	Deposit_dexponentmantissap1(srcdstA,				\
    108        1.1  fredette 	    (Dexponentmantissap1(srcdstA)>>(varamount)));			\
    109        1.1  fredette     } }
    110        1.1  fredette /* varamount must be less than 64 */
    111        1.1  fredette #define Dbl_leftshift(srcdstA, srcdstB, varamount)			\
    112        1.1  fredette     {if((varamount) >= 32) {						\
    113        1.1  fredette 	Dallp1(srcdstA) = Dallp2(srcdstB) << (varamount-32);		\
    114        1.1  fredette 	Dallp2(srcdstB)=0;						\
    115        1.1  fredette     }									\
    116        1.1  fredette     else {								\
    117        1.1  fredette 	if ((varamount) > 0) {						\
    118        1.1  fredette 	    Dallp1(srcdstA) = (Dallp1(srcdstA) << (varamount)) |	\
    119        1.1  fredette 		(Dallp2(srcdstB) >> (32-(varamount)));			\
    120        1.1  fredette 	    Dallp2(srcdstB) <<= varamount;				\
    121        1.1  fredette 	}								\
    122        1.1  fredette     } }
    123        1.1  fredette #define Dbl_leftshiftby1_withextent(lefta,leftb,right,resulta,resultb)	\
    124        1.1  fredette     Shiftdouble(Dallp1(lefta), Dallp2(leftb), 31, Dallp1(resulta));	\
    125        1.1  fredette     Shiftdouble(Dallp2(leftb), Extall(right), 31, Dallp2(resultb))
    126        1.1  fredette 
    127        1.1  fredette #define Dbl_rightshiftby1_withextent(leftb,right,dst)		\
    128        1.1  fredette     Extall(dst) = (Dallp2(leftb) << 31) | ((unsigned)Extall(right) >> 1) | \
    129        1.1  fredette 		  Extlow(right)
    130        1.1  fredette 
    131        1.1  fredette #define Dbl_arithrightshiftby1(srcdstA,srcdstB)			\
    132        1.1  fredette     Shiftdouble(Dallp1(srcdstA),Dallp2(srcdstB),1,Dallp2(srcdstB));\
    133        1.1  fredette     Dallp1(srcdstA) = (int)Dallp1(srcdstA) >> 1
    134        1.1  fredette 
    135        1.1  fredette /* Sign extend the sign bit with an integer destination */
    136        1.1  fredette #define Dbl_signextendedsign(value)  Dsignedsign(value)
    137        1.1  fredette 
    138        1.1  fredette #define Dbl_isone_hidden(dbl_value) (Is_dhidden(dbl_value)!=0)
    139        1.1  fredette /* Singles and doubles may include the sign and exponent fields.  The
    140        1.1  fredette  * hidden bit and the hidden overflow must be included. */
    141        1.1  fredette #define Dbl_increment(dbl_valueA,dbl_valueB) \
    142        1.1  fredette     if( (Dallp2(dbl_valueB) += 1) == 0 )  Dallp1(dbl_valueA) += 1
    143        1.1  fredette #define Dbl_increment_mantissa(dbl_valueA,dbl_valueB) \
    144        1.1  fredette     if( (Dmantissap2(dbl_valueB) += 1) == 0 )  \
    145        1.1  fredette     Deposit_dmantissap1(dbl_valueA,dbl_valueA+1)
    146        1.1  fredette #define Dbl_decrement(dbl_valueA,dbl_valueB) \
    147        1.1  fredette     if( Dallp2(dbl_valueB) == 0 )  Dallp1(dbl_valueA) -= 1; \
    148        1.1  fredette     Dallp2(dbl_valueB) -= 1
    149        1.1  fredette 
    150        1.1  fredette #define Dbl_isone_sign(dbl_value) (Is_dsign(dbl_value)!=0)
    151        1.1  fredette #define Dbl_isone_hiddenoverflow(dbl_value) (Is_dhiddenoverflow(dbl_value)!=0)
    152        1.1  fredette #define Dbl_isone_lowmantissap1(dbl_valueA) (Is_dlowp1(dbl_valueA)!=0)
    153        1.1  fredette #define Dbl_isone_lowmantissap2(dbl_valueB) (Is_dlowp2(dbl_valueB)!=0)
    154        1.1  fredette #define Dbl_isone_signaling(dbl_value) (Is_dsignaling(dbl_value)!=0)
    155        1.1  fredette #define Dbl_is_signalingnan(dbl_value) (Dsignalingnan(dbl_value)==0xfff)
    156        1.1  fredette #define Dbl_isnotzero(dbl_valueA,dbl_valueB) \
    157        1.1  fredette     (Dallp1(dbl_valueA) || Dallp2(dbl_valueB))
    158        1.1  fredette #define Dbl_isnotzero_hiddenhigh7mantissa(dbl_value) \
    159        1.1  fredette     (Dhiddenhigh7mantissa(dbl_value)!=0)
    160        1.1  fredette #define Dbl_isnotzero_exponent(dbl_value) (Dexponent(dbl_value)!=0)
    161        1.1  fredette #define Dbl_isnotzero_mantissa(dbl_valueA,dbl_valueB) \
    162        1.1  fredette     (Dmantissap1(dbl_valueA) || Dmantissap2(dbl_valueB))
    163        1.1  fredette #define Dbl_isnotzero_mantissap1(dbl_valueA) (Dmantissap1(dbl_valueA)!=0)
    164        1.1  fredette #define Dbl_isnotzero_mantissap2(dbl_valueB) (Dmantissap2(dbl_valueB)!=0)
    165        1.1  fredette #define Dbl_isnotzero_exponentmantissa(dbl_valueA,dbl_valueB) \
    166        1.1  fredette     (Dexponentmantissap1(dbl_valueA) || Dmantissap2(dbl_valueB))
    167        1.1  fredette #define Dbl_isnotzero_low4p2(dbl_value) (Dlow4p2(dbl_value)!=0)
    168        1.1  fredette #define Dbl_iszero(dbl_valueA,dbl_valueB) (Dallp1(dbl_valueA)==0 && \
    169        1.1  fredette     Dallp2(dbl_valueB)==0)
    170        1.1  fredette #define Dbl_iszero_allp1(dbl_value) (Dallp1(dbl_value)==0)
    171        1.1  fredette #define Dbl_iszero_allp2(dbl_value) (Dallp2(dbl_value)==0)
    172        1.1  fredette #define Dbl_iszero_hidden(dbl_value) (Is_dhidden(dbl_value)==0)
    173        1.1  fredette #define Dbl_iszero_hiddenoverflow(dbl_value) (Is_dhiddenoverflow(dbl_value)==0)
    174        1.1  fredette #define Dbl_iszero_hiddenhigh3mantissa(dbl_value) \
    175        1.1  fredette     (Dhiddenhigh3mantissa(dbl_value)==0)
    176        1.1  fredette #define Dbl_iszero_hiddenhigh7mantissa(dbl_value) \
    177        1.1  fredette     (Dhiddenhigh7mantissa(dbl_value)==0)
    178        1.1  fredette #define Dbl_iszero_sign(dbl_value) (Is_dsign(dbl_value)==0)
    179        1.1  fredette #define Dbl_iszero_exponent(dbl_value) (Dexponent(dbl_value)==0)
    180        1.1  fredette #define Dbl_iszero_mantissa(dbl_valueA,dbl_valueB) \
    181        1.1  fredette     (Dmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0)
    182        1.1  fredette #define Dbl_iszero_exponentmantissa(dbl_valueA,dbl_valueB) \
    183        1.1  fredette     (Dexponentmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0)
    184        1.1  fredette #define Dbl_isinfinity_exponent(dbl_value)		\
    185        1.1  fredette     (Dexponent(dbl_value)==DBL_INFINITY_EXPONENT)
    186        1.1  fredette #define Dbl_isnotinfinity_exponent(dbl_value)		\
    187        1.1  fredette     (Dexponent(dbl_value)!=DBL_INFINITY_EXPONENT)
    188        1.1  fredette #define Dbl_isinfinity(dbl_valueA,dbl_valueB)			\
    189        1.1  fredette     (Dexponent(dbl_valueA)==DBL_INFINITY_EXPONENT &&	\
    190        1.1  fredette     Dmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0)
    191        1.1  fredette #define Dbl_isnan(dbl_valueA,dbl_valueB)		\
    192        1.1  fredette     (Dexponent(dbl_valueA)==DBL_INFINITY_EXPONENT &&	\
    193        1.1  fredette     (Dmantissap1(dbl_valueA)!=0 || Dmantissap2(dbl_valueB)!=0))
    194        1.1  fredette #define Dbl_isnotnan(dbl_valueA,dbl_valueB)		\
    195        1.1  fredette     (Dexponent(dbl_valueA)!=DBL_INFINITY_EXPONENT ||	\
    196        1.1  fredette     (Dmantissap1(dbl_valueA)==0 && Dmantissap2(dbl_valueB)==0))
    197        1.1  fredette 
    198        1.1  fredette #define Dbl_islessthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
    199        1.1  fredette     (Dallp1(dbl_op1a) < Dallp1(dbl_op2a) ||			\
    200        1.1  fredette      (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
    201        1.1  fredette       Dallp2(dbl_op1b) < Dallp2(dbl_op2b)))
    202        1.1  fredette #define Dbl_isgreaterthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
    203        1.1  fredette     (Dallp1(dbl_op1a) > Dallp1(dbl_op2a) ||			\
    204        1.1  fredette      (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
    205        1.1  fredette       Dallp2(dbl_op1b) > Dallp2(dbl_op2b)))
    206        1.1  fredette #define Dbl_isnotlessthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
    207        1.1  fredette     (Dallp1(dbl_op1a) > Dallp1(dbl_op2a) ||			\
    208        1.1  fredette      (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
    209        1.1  fredette       Dallp2(dbl_op1b) >= Dallp2(dbl_op2b)))
    210        1.1  fredette #define Dbl_isnotgreaterthan(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b) \
    211        1.1  fredette     (Dallp1(dbl_op1a) < Dallp1(dbl_op2a) ||			\
    212        1.1  fredette      (Dallp1(dbl_op1a) == Dallp1(dbl_op2a) &&			\
    213        1.1  fredette       Dallp2(dbl_op1b) <= Dallp2(dbl_op2b)))
    214        1.1  fredette #define Dbl_isequal(dbl_op1a,dbl_op1b,dbl_op2a,dbl_op2b)	\
    215        1.1  fredette      ((Dallp1(dbl_op1a) == Dallp1(dbl_op2a)) &&			\
    216        1.1  fredette       (Dallp2(dbl_op1b) == Dallp2(dbl_op2b)))
    217        1.1  fredette 
    218        1.1  fredette #define Dbl_leftshiftby8(dbl_valueA,dbl_valueB) \
    219        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),24,Dallp1(dbl_valueA)); \
    220        1.1  fredette     Dallp2(dbl_valueB) <<= 8
    221        1.1  fredette #define Dbl_leftshiftby7(dbl_valueA,dbl_valueB) \
    222        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),25,Dallp1(dbl_valueA)); \
    223        1.1  fredette     Dallp2(dbl_valueB) <<= 7
    224        1.1  fredette #define Dbl_leftshiftby4(dbl_valueA,dbl_valueB) \
    225        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),28,Dallp1(dbl_valueA)); \
    226        1.1  fredette     Dallp2(dbl_valueB) <<= 4
    227        1.1  fredette #define Dbl_leftshiftby3(dbl_valueA,dbl_valueB) \
    228        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),29,Dallp1(dbl_valueA)); \
    229        1.1  fredette     Dallp2(dbl_valueB) <<= 3
    230        1.1  fredette #define Dbl_leftshiftby2(dbl_valueA,dbl_valueB) \
    231        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),30,Dallp1(dbl_valueA)); \
    232        1.1  fredette     Dallp2(dbl_valueB) <<= 2
    233        1.1  fredette #define Dbl_leftshiftby1(dbl_valueA,dbl_valueB) \
    234        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),31,Dallp1(dbl_valueA)); \
    235        1.1  fredette     Dallp2(dbl_valueB) <<= 1
    236        1.1  fredette 
    237        1.1  fredette #define Dbl_rightshiftby8(dbl_valueA,dbl_valueB) \
    238        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),8,Dallp2(dbl_valueB)); \
    239        1.1  fredette     Dallp1(dbl_valueA) >>= 8
    240        1.1  fredette #define Dbl_rightshiftby4(dbl_valueA,dbl_valueB) \
    241        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),4,Dallp2(dbl_valueB)); \
    242        1.1  fredette     Dallp1(dbl_valueA) >>= 4
    243        1.1  fredette #define Dbl_rightshiftby2(dbl_valueA,dbl_valueB) \
    244        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),2,Dallp2(dbl_valueB)); \
    245        1.1  fredette     Dallp1(dbl_valueA) >>= 2
    246        1.1  fredette #define Dbl_rightshiftby1(dbl_valueA,dbl_valueB) \
    247        1.1  fredette     Shiftdouble(Dallp1(dbl_valueA),Dallp2(dbl_valueB),1,Dallp2(dbl_valueB)); \
    248        1.1  fredette     Dallp1(dbl_valueA) >>= 1
    249        1.1  fredette 
    250        1.1  fredette /* This magnitude comparison uses the signless first words and
    251        1.1  fredette  * the regular part2 words.  The comparison is graphically:
    252        1.1  fredette  *
    253        1.1  fredette  *       1st greater?  -------------
    254        1.1  fredette  *				   |
    255        1.1  fredette  *       1st less?-----------------+---------
    256        1.1  fredette  *				   |	    |
    257        1.1  fredette  *       2nd greater or equal----->|	    |
    258        1.1  fredette  *				 False     True
    259        1.1  fredette  */
    260        1.1  fredette #define Dbl_ismagnitudeless(leftB,rightB,signlessleft,signlessright)	\
    261        1.1  fredette       ((signlessleft <= signlessright) &&				\
    262        1.1  fredette        ( (signlessleft < signlessright) || (Dallp2(leftB)<Dallp2(rightB)) ))
    263        1.1  fredette 
    264        1.1  fredette #define Dbl_copytoint_exponentmantissap1(src,dest) \
    265        1.1  fredette     dest = Dexponentmantissap1(src)
    266        1.1  fredette 
    267        1.1  fredette /* A quiet NaN has the high mantissa bit clear and at least on other (in this
    268        1.1  fredette  * case the adjacent bit) bit set. */
    269        1.1  fredette #define Dbl_set_quiet(dbl_value) Deposit_dhigh2mantissa(dbl_value,1)
    270        1.1  fredette #define Dbl_set_exponent(dbl_value, exp) Deposit_dexponent(dbl_value,exp)
    271        1.1  fredette 
    272        1.1  fredette #define Dbl_set_mantissa(desta,destb,valuea,valueb)	\
    273        1.1  fredette     Deposit_dmantissap1(desta,valuea);			\
    274        1.1  fredette     Dmantissap2(destb) = Dmantissap2(valueb)
    275        1.1  fredette #define Dbl_set_mantissap1(desta,valuea)		\
    276        1.1  fredette     Deposit_dmantissap1(desta,valuea)
    277        1.1  fredette #define Dbl_set_mantissap2(destb,valueb)		\
    278        1.1  fredette     Dmantissap2(destb) = Dmantissap2(valueb)
    279        1.1  fredette 
    280        1.1  fredette #define Dbl_set_exponentmantissa(desta,destb,valuea,valueb)	\
    281        1.1  fredette     Deposit_dexponentmantissap1(desta,valuea);			\
    282        1.1  fredette     Dmantissap2(destb) = Dmantissap2(valueb)
    283        1.1  fredette #define Dbl_set_exponentmantissap1(dest,value)			\
    284        1.1  fredette     Deposit_dexponentmantissap1(dest,value)
    285        1.1  fredette 
    286        1.1  fredette #define Dbl_copyfromptr(src,desta,destb) \
    287        1.1  fredette     Dallp1(desta) = src->wd0;		\
    288        1.1  fredette     Dallp2(destb) = src->wd1
    289        1.1  fredette #define Dbl_copytoptr(srca,srcb,dest)	\
    290        1.1  fredette     dest->wd0 = Dallp1(srca);		\
    291        1.1  fredette     dest->wd1 = Dallp2(srcb)
    292        1.1  fredette 
    293        1.1  fredette /*  An infinity is represented with the max exponent and a zero mantissa */
    294        1.1  fredette #define Dbl_setinfinity_exponent(dbl_value) \
    295        1.1  fredette     Deposit_dexponent(dbl_value,DBL_INFINITY_EXPONENT)
    296        1.1  fredette #define Dbl_setinfinity_exponentmantissa(dbl_valueA,dbl_valueB)	\
    297        1.1  fredette     Deposit_dexponentmantissap1(dbl_valueA,			\
    298        1.1  fredette     (DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH))));	\
    299        1.1  fredette     Dmantissap2(dbl_valueB) = 0
    300        1.1  fredette #define Dbl_setinfinitypositive(dbl_valueA,dbl_valueB)		\
    301        1.1  fredette     Dallp1(dbl_valueA)						\
    302        1.1  fredette 	= (DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH)));	\
    303        1.1  fredette     Dmantissap2(dbl_valueB) = 0
    304        1.1  fredette #define Dbl_setinfinitynegative(dbl_valueA,dbl_valueB)		\
    305        1.1  fredette     Dallp1(dbl_valueA) = (1<<31) |				\
    306        1.1  fredette 	(DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH)));	\
    307        1.1  fredette     Dmantissap2(dbl_valueB) = 0
    308        1.1  fredette #define Dbl_setinfinity(dbl_valueA,dbl_valueB,sign)		\
    309        1.1  fredette     Dallp1(dbl_valueA) = (sign << 31) |				\
    310        1.1  fredette 	(DBL_INFINITY_EXPONENT << (32-(1+DBL_EXP_LENGTH)));	\
    311        1.1  fredette     Dmantissap2(dbl_valueB) = 0
    312        1.1  fredette 
    313        1.1  fredette #define Dbl_sethigh4bits(dbl_value, extsign) Deposit_dhigh4p1(dbl_value,extsign)
    314        1.1  fredette #define Dbl_set_sign(dbl_value,sign) Deposit_dsign(dbl_value,sign)
    315        1.1  fredette #define Dbl_invert_sign(dbl_value) Deposit_dsign(dbl_value,~Dsign(dbl_value))
    316        1.1  fredette #define Dbl_setone_sign(dbl_value) Deposit_dsign(dbl_value,1)
    317        1.1  fredette #define Dbl_setone_lowmantissap2(dbl_value) Deposit_dlowp2(dbl_value,1)
    318        1.1  fredette #define Dbl_setzero_sign(dbl_value) Dallp1(dbl_value) &= 0x7fffffff
    319        1.1  fredette #define Dbl_setzero_exponent(dbl_value)			\
    320        1.1  fredette     Dallp1(dbl_value) &= 0x800fffff
    321        1.1  fredette #define Dbl_setzero_mantissa(dbl_valueA,dbl_valueB)	\
    322        1.1  fredette     Dallp1(dbl_valueA) &= 0xfff00000;			\
    323        1.1  fredette     Dallp2(dbl_valueB) = 0
    324        1.1  fredette #define Dbl_setzero_mantissap1(dbl_value) Dallp1(dbl_value) &= 0xfff00000
    325        1.1  fredette #define Dbl_setzero_mantissap2(dbl_value) Dallp2(dbl_value) = 0
    326        1.1  fredette #define Dbl_setzero_exponentmantissa(dbl_valueA,dbl_valueB)	\
    327        1.1  fredette     Dallp1(dbl_valueA) &= 0x80000000;		\
    328        1.1  fredette     Dallp2(dbl_valueB) = 0
    329        1.1  fredette #define Dbl_setzero_exponentmantissap1(dbl_valueA)	\
    330        1.1  fredette     Dallp1(dbl_valueA) &= 0x80000000
    331        1.1  fredette #define Dbl_setzero(dbl_valueA,dbl_valueB) \
    332        1.1  fredette     Dallp1(dbl_valueA) = 0; Dallp2(dbl_valueB) = 0
    333        1.1  fredette #define Dbl_setzerop1(dbl_value) Dallp1(dbl_value) = 0
    334        1.1  fredette #define Dbl_setzerop2(dbl_value) Dallp2(dbl_value) = 0
    335        1.1  fredette #define Dbl_setnegativezero(dbl_value) \
    336        1.1  fredette     Dallp1(dbl_value) = 1 << 31; Dallp2(dbl_value) = 0
    337        1.1  fredette #define Dbl_setnegativezerop1(dbl_value) Dallp1(dbl_value) = 1 << 31
    338        1.1  fredette 
    339        1.1  fredette /* Use the following macro for both overflow & underflow conditions */
    340        1.1  fredette #define ovfl -
    341        1.1  fredette #define unfl +
    342        1.1  fredette #define Dbl_setwrapped_exponent(dbl_value,exponent,op) \
    343        1.1  fredette     Deposit_dexponent(dbl_value,(exponent op DBL_WRAP))
    344        1.1  fredette 
    345        1.1  fredette #define Dbl_setlargestpositive(dbl_valueA,dbl_valueB)			\
    346        1.1  fredette     Dallp1(dbl_valueA) = ((DBL_MAX_EXP+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH))) \
    347        1.1  fredette 			| ((1<<(32-(1+DBL_EXP_LENGTH))) - 1 );		\
    348        1.1  fredette     Dallp2(dbl_valueB) = 0xFFFFFFFF
    349        1.1  fredette #define Dbl_setlargestnegative(dbl_valueA,dbl_valueB)			\
    350        1.1  fredette     Dallp1(dbl_valueA) = ((DBL_MAX_EXP+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH))) \
    351        1.1  fredette 			| ((1<<(32-(1+DBL_EXP_LENGTH))) - 1 ) | (1<<31); \
    352        1.1  fredette     Dallp2(dbl_valueB) = 0xFFFFFFFF
    353        1.1  fredette #define Dbl_setlargest_exponentmantissa(dbl_valueA,dbl_valueB)		\
    354        1.1  fredette     Deposit_dexponentmantissap1(dbl_valueA,				\
    355        1.1  fredette 	(((DBL_MAX_EXP+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH)))		\
    356        1.1  fredette 			| ((1<<(32-(1+DBL_EXP_LENGTH))) - 1 )));	\
    357        1.1  fredette     Dallp2(dbl_valueB) = 0xFFFFFFFF
    358        1.1  fredette 
    359        1.1  fredette #define Dbl_setnegativeinfinity(dbl_valueA,dbl_valueB)			\
    360        1.1  fredette     Dallp1(dbl_valueA) = ((1<<DBL_EXP_LENGTH) | DBL_INFINITY_EXPONENT)	\
    361        1.1  fredette 			 << (32-(1+DBL_EXP_LENGTH)) ;			\
    362        1.1  fredette     Dallp2(dbl_valueB) = 0
    363        1.1  fredette #define Dbl_setlargest(dbl_valueA,dbl_valueB,sign)			\
    364        1.1  fredette     Dallp1(dbl_valueA) = (sign << 31) |					\
    365        1.1  fredette 	((DBL_MAX_EXP+DBL_BIAS) << (32-(1+DBL_EXP_LENGTH))) |		\
    366        1.1  fredette 	 ((1 << (32-(1+DBL_EXP_LENGTH))) - 1 );				\
    367        1.1  fredette     Dallp2(dbl_valueB) = 0xFFFFFFFF
    368        1.1  fredette 
    369        1.1  fredette 
    370        1.1  fredette /* The high bit is always zero so arithmetic or logical shifts will work. */
    371        1.1  fredette #define Dbl_right_align(srcdstA,srcdstB,shift,extent)			\
    372        1.1  fredette     if( shift >= 32 )							\
    373        1.1  fredette 	{								\
    374        1.1  fredette 	/* Big shift requires examining the portion shift off		\
    375        1.1  fredette 	the end to properly set inexact.  */				\
    376        1.1  fredette 	if(shift < 64)							\
    377        1.1  fredette 	    {								\
    378        1.1  fredette 	    if(shift > 32)						\
    379        1.1  fredette 		{							\
    380        1.1  fredette 		Variable_shift_double(Dallp1(srcdstA),Dallp2(srcdstB),	\
    381        1.1  fredette 		 shift-32, Extall(extent));				\
    382        1.1  fredette 		if(Dallp2(srcdstB) << (64 - (shift))) Ext_setone_low(extent); \
    383        1.1  fredette 		}							\
    384        1.1  fredette 	    else Extall(extent) = Dallp2(srcdstB);			\
    385        1.1  fredette 	    Dallp2(srcdstB) = Dallp1(srcdstA) >> (shift - 32);		\
    386        1.1  fredette 	    }								\
    387        1.1  fredette 	else								\
    388        1.1  fredette 	    {								\
    389        1.1  fredette 	    Extall(extent) = Dallp1(srcdstA);				\
    390        1.1  fredette 	    if(Dallp2(srcdstB)) Ext_setone_low(extent);			\
    391        1.1  fredette 	    Dallp2(srcdstB) = 0;					\
    392        1.1  fredette 	    }								\
    393        1.1  fredette 	Dallp1(srcdstA) = 0;						\
    394        1.1  fredette 	}								\
    395        1.1  fredette     else								\
    396        1.1  fredette 	{								\
    397        1.1  fredette 	/* Small alignment is simpler.  Extension is easily set. */	\
    398        1.1  fredette 	if (shift > 0)							\
    399        1.1  fredette 	    {								\
    400        1.1  fredette 	    Extall(extent) = Dallp2(srcdstB) << (32 - (shift));		\
    401        1.1  fredette 	    Variable_shift_double(Dallp1(srcdstA),Dallp2(srcdstB),shift, \
    402        1.1  fredette 	     Dallp2(srcdstB));						\
    403        1.1  fredette 	    Dallp1(srcdstA) >>= shift;					\
    404        1.1  fredette 	    }								\
    405        1.1  fredette 	else Extall(extent) = 0;					\
    406        1.1  fredette 	}
    407        1.1  fredette 
    408        1.1  fredette /*
    409        1.1  fredette  * Here we need to shift the result right to correct for an overshift
    410        1.1  fredette  * (due to the exponent becoming negative) during normalization.
    411        1.1  fredette  */
    412        1.1  fredette #define Dbl_fix_overshift(srcdstA,srcdstB,shift,extent)			\
    413        1.1  fredette 	    Extall(extent) = Dallp2(srcdstB) << (32 - (shift));		\
    414        1.1  fredette 	    Dallp2(srcdstB) = (Dallp1(srcdstA) << (32 - (shift))) |	\
    415        1.1  fredette 		(Dallp2(srcdstB) >> (shift));				\
    416        1.1  fredette 	    Dallp1(srcdstA) = Dallp1(srcdstA) >> shift
    417        1.1  fredette 
    418        1.1  fredette #define Dbl_hiddenhigh3mantissa(dbl_value) Dhiddenhigh3mantissa(dbl_value)
    419        1.1  fredette #define Dbl_hidden(dbl_value) Dhidden(dbl_value)
    420        1.1  fredette #define Dbl_lowmantissap2(dbl_value) Dlowp2(dbl_value)
    421        1.1  fredette 
    422        1.1  fredette /* The left argument is never smaller than the right argument */
    423        1.1  fredette #define Dbl_subtract(lefta,leftb,righta,rightb,resulta,resultb)			\
    424        1.1  fredette     if( Dallp2(rightb) > Dallp2(leftb) ) Dallp1(lefta)--;	\
    425        1.1  fredette     Dallp2(resultb) = Dallp2(leftb) - Dallp2(rightb);		\
    426        1.1  fredette     Dallp1(resulta) = Dallp1(lefta) - Dallp1(righta)
    427        1.1  fredette 
    428        1.1  fredette /* Subtract right augmented with extension from left augmented with zeros and
    429        1.1  fredette  * store into result and extension. */
    430        1.1  fredette #define Dbl_subtract_withextension(lefta,leftb,righta,rightb,extent,resulta,resultb)	\
    431        1.1  fredette     Dbl_subtract(lefta,leftb,righta,rightb,resulta,resultb);		\
    432        1.1  fredette     if( (Extall(extent) = 0-Extall(extent)) )				\
    433        1.1  fredette 	{								\
    434        1.1  fredette 	if((Dallp2(resultb)--) == 0) Dallp1(resulta)--;			\
    435        1.1  fredette 	}
    436        1.1  fredette 
    437        1.1  fredette #define Dbl_addition(lefta,leftb,righta,rightb,resulta,resultb)		\
    438        1.1  fredette     /* If the sum of the low words is less than either source, then	\
    439        1.1  fredette      * an overflow into the next word occurred. */			\
    440        1.1  fredette     Dallp1(resulta) = Dallp1(lefta) + Dallp1(righta);			\
    441        1.1  fredette     if((Dallp2(resultb) = Dallp2(leftb) + Dallp2(rightb)) < Dallp2(rightb)) \
    442        1.1  fredette 	Dallp1(resulta)++
    443        1.1  fredette 
    444        1.1  fredette #define Dbl_xortointp1(left,right,result)			\
    445        1.1  fredette     result = Dallp1(left) XOR Dallp1(right)
    446        1.1  fredette 
    447        1.1  fredette #define Dbl_xorfromintp1(left,right,result)			\
    448        1.1  fredette     Dallp1(result) = left XOR Dallp1(right)
    449        1.1  fredette 
    450        1.1  fredette #define Dbl_swap_lower(left,right)				\
    451        1.1  fredette     Dallp2(left)  = Dallp2(left) XOR Dallp2(right);		\
    452        1.1  fredette     Dallp2(right) = Dallp2(left) XOR Dallp2(right);		\
    453        1.1  fredette     Dallp2(left)  = Dallp2(left) XOR Dallp2(right)
    454        1.1  fredette 
    455        1.1  fredette /* Need to Initialize */
    456        1.1  fredette #define Dbl_makequietnan(desta,destb)					\
    457        1.1  fredette     Dallp1(desta) = ((DBL_MAX_EXP+DBL_BIAS)+1)<< (32-(1+DBL_EXP_LENGTH))	\
    458        1.1  fredette 		| (1<<(32-(1+DBL_EXP_LENGTH+2)));			\
    459        1.1  fredette     Dallp2(destb) = 0
    460        1.1  fredette #define Dbl_makesignalingnan(desta,destb)				\
    461        1.1  fredette     Dallp1(desta) = ((DBL_MAX_EXP+DBL_BIAS)+1)<< (32-(1+DBL_EXP_LENGTH))	\
    462        1.1  fredette 		| (1<<(32-(1+DBL_EXP_LENGTH+1)));			\
    463        1.1  fredette     Dallp2(destb) = 0
    464        1.1  fredette 
    465        1.1  fredette #define Dbl_normalize(dbl_opndA,dbl_opndB,exponent)			\
    466        1.1  fredette 	while(Dbl_iszero_hiddenhigh7mantissa(dbl_opndA)) {		\
    467        1.1  fredette 		Dbl_leftshiftby8(dbl_opndA,dbl_opndB);			\
    468        1.1  fredette 		exponent -= 8;						\
    469        1.1  fredette 	}								\
    470        1.1  fredette 	if(Dbl_iszero_hiddenhigh3mantissa(dbl_opndA)) {			\
    471        1.1  fredette 		Dbl_leftshiftby4(dbl_opndA,dbl_opndB);			\
    472        1.1  fredette 		exponent -= 4;						\
    473        1.1  fredette 	}								\
    474        1.1  fredette 	while(Dbl_iszero_hidden(dbl_opndA)) {				\
    475        1.1  fredette 		Dbl_leftshiftby1(dbl_opndA,dbl_opndB);			\
    476        1.1  fredette 		exponent -= 1;						\
    477        1.1  fredette 	}
    478        1.1  fredette 
    479        1.1  fredette #define Twoword_add(src1dstA,src1dstB,src2A,src2B)		\
    480        1.1  fredette 	/*							\
    481        1.1  fredette 	 * want this macro to generate:				\
    482        1.1  fredette 	 *	ADD	src1dstB,src2B,src1dstB;		\
    483        1.1  fredette 	 *	ADDC	src1dstA,src2A,src1dstA;		\
    484        1.1  fredette 	 */							\
    485        1.1  fredette 	if ((src1dstB) + (src2B) < (src1dstB)) Dallp1(src1dstA)++; \
    486        1.1  fredette 	Dallp1(src1dstA) += (src2A);				\
    487        1.1  fredette 	Dallp2(src1dstB) += (src2B)
    488        1.1  fredette 
    489        1.1  fredette #define Twoword_subtract(src1dstA,src1dstB,src2A,src2B)		\
    490        1.1  fredette 	/*							\
    491        1.1  fredette 	 * want this macro to generate:				\
    492        1.1  fredette 	 *	SUB	src1dstB,src2B,src1dstB;		\
    493        1.1  fredette 	 *	SUBB	src1dstA,src2A,src1dstA;		\
    494        1.1  fredette 	 */							\
    495        1.1  fredette 	if ((src1dstB) < (src2B)) Dallp1(src1dstA)--;		\
    496        1.1  fredette 	Dallp1(src1dstA) -= (src2A);				\
    497        1.1  fredette 	Dallp2(src1dstB) -= (src2B)
    498        1.1  fredette 
    499        1.1  fredette #define Dbl_setoverflow(resultA,resultB)				\
    500        1.1  fredette 	/* set result to infinity or largest number */			\
    501        1.1  fredette 	switch (Rounding_mode()) {					\
    502        1.1  fredette 		case ROUNDPLUS:						\
    503        1.1  fredette 			if (Dbl_isone_sign(resultA)) {			\
    504        1.1  fredette 				Dbl_setlargestnegative(resultA,resultB); \
    505        1.1  fredette 			}						\
    506        1.1  fredette 			else {						\
    507        1.1  fredette 				Dbl_setinfinitypositive(resultA,resultB); \
    508        1.1  fredette 			}						\
    509        1.1  fredette 			break;						\
    510        1.1  fredette 		case ROUNDMINUS:					\
    511        1.1  fredette 			if (Dbl_iszero_sign(resultA)) {			\
    512        1.1  fredette 				Dbl_setlargestpositive(resultA,resultB); \
    513        1.1  fredette 			}						\
    514        1.1  fredette 			else {						\
    515        1.1  fredette 				Dbl_setinfinitynegative(resultA,resultB); \
    516        1.1  fredette 			}						\
    517        1.1  fredette 			break;						\
    518        1.1  fredette 		case ROUNDNEAREST:					\
    519        1.1  fredette 			Dbl_setinfinity_exponentmantissa(resultA,resultB); \
    520        1.1  fredette 			break;						\
    521        1.1  fredette 		case ROUNDZERO:						\
    522        1.1  fredette 			Dbl_setlargest_exponentmantissa(resultA,resultB); \
    523        1.1  fredette 	}
    524        1.1  fredette 
    525        1.1  fredette #define Dbl_denormalize(opndp1,opndp2,exponent,guard,sticky,inexact)	\
    526        1.1  fredette     Dbl_clear_signexponent_set_hidden(opndp1);				\
    527        1.1  fredette     if (exponent >= (1-DBL_P)) {					\
    528        1.1  fredette 	if (exponent >= -31) {						\
    529        1.1  fredette 	    guard = (Dallp2(opndp2) >> (-(exponent))) & 1;		\
    530        1.1  fredette 	    if (exponent < 0) sticky |= Dallp2(opndp2) << (32+exponent); \
    531        1.1  fredette 	    if (exponent > -31) {					\
    532        1.1  fredette 		Variable_shift_double(opndp1,opndp2,1-exponent,opndp2);	\
    533        1.1  fredette 		Dallp1(opndp1) >>= 1-exponent;				\
    534        1.1  fredette 	    }								\
    535        1.1  fredette 	    else {							\
    536        1.1  fredette 		Dallp2(opndp2) = Dallp1(opndp1);			\
    537        1.1  fredette 		Dbl_setzerop1(opndp1);					\
    538        1.1  fredette 	    }								\
    539        1.1  fredette 	}								\
    540        1.1  fredette 	else {								\
    541        1.1  fredette 	    guard = (Dallp1(opndp1) >> (-32-(exponent))) & 1;		\
    542        1.1  fredette 	    if (exponent == -32) sticky |= Dallp2(opndp2);		\
    543        1.1  fredette 	    else sticky |= (Dallp2(opndp2) | Dallp1(opndp1) << (64+(exponent))); \
    544        1.1  fredette 	    Dallp2(opndp2) = Dallp1(opndp1) >> (-31-(exponent));	\
    545        1.1  fredette 	    Dbl_setzerop1(opndp1);					\
    546        1.1  fredette 	}								\
    547        1.1  fredette 	inexact = guard | sticky;					\
    548        1.1  fredette     }									\
    549        1.1  fredette     else {								\
    550        1.1  fredette 	guard = 0;							\
    551        1.1  fredette 	sticky |= (Dallp1(opndp1) | Dallp2(opndp2));			\
    552        1.1  fredette 	Dbl_setzero(opndp1,opndp2);					\
    553        1.1  fredette 	inexact = sticky;						\
    554        1.1  fredette     }
    555        1.1  fredette 
    556        1.1  fredette 
    557  1.1.112.1       mjf int dbl_fadd(dbl_floating_point *, dbl_floating_point*, dbl_floating_point*, unsigned int *);
    558  1.1.112.1       mjf int dbl_fcmp(dbl_floating_point *, dbl_floating_point*, unsigned int, unsigned int *);
    559  1.1.112.1       mjf int dbl_fdiv(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);
    560  1.1.112.1       mjf int dbl_fmpy(dbl_floating_point *, dbl_floating_point *, dbl_floating_point*, unsigned int *);
    561  1.1.112.1       mjf int dbl_frem(dbl_floating_point *, dbl_floating_point *, dbl_floating_point*, unsigned int *);
    562  1.1.112.1       mjf int dbl_fsqrt(dbl_floating_point *, dbl_floating_point *, unsigned int *);
    563  1.1.112.1       mjf int dbl_fsub(dbl_floating_point *, dbl_floating_point *, dbl_floating_point*, unsigned int *);
    564        1.1  fredette 
    565  1.1.112.1       mjf dbl_floating_point dbl_setoverflow(unsigned int);
    566        1.1  fredette 
    567  1.1.112.1       mjf int sgl_to_dbl_fcnvff(sgl_floating_point *, dbl_floating_point *, unsigned int *);
    568  1.1.112.1       mjf int dbl_to_sgl_fcnvff(dbl_floating_point *, sgl_floating_point *, unsigned int *);
    569        1.1  fredette 
    570  1.1.112.1       mjf int dbl_frnd(dbl_floating_point *, dbl_floating_point *, unsigned int *);
    571