Home | History | Annotate | Line # | Download | only in spmath
      1  1.4    andvar /*	$NetBSD: dbl_float.h,v 1.4 2022/05/28 22:08:47 andvar Exp $	*/
      2  1.1  fredette 
      3  1.3     skrll /*	$OpenBSD: dbl_float.h,v 1.10 2004/01/02 14:39:01 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.4    andvar /* 32-bit word grabbing 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.3     skrll     Dallp1(dbl_valueA) = ((DBL_EMAX+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.3     skrll     Dallp1(dbl_valueA) = ((DBL_EMAX+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.3     skrll 	(((DBL_EMAX+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.3     skrll 	((DBL_EMAX+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.3     skrll     Dallp1(desta) = ((DBL_EMAX+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.3     skrll     Dallp1(desta) = ((DBL_EMAX+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.2     skrll int dbl_fadd(dbl_floating_point *, dbl_floating_point*, dbl_floating_point*, unsigned int *);
    558  1.2     skrll int dbl_fcmp(dbl_floating_point *, dbl_floating_point*, unsigned int, unsigned int *);
    559  1.2     skrll int dbl_fdiv(dbl_floating_point *, dbl_floating_point *, dbl_floating_point *, unsigned int *);
    560  1.2     skrll int dbl_fmpy(dbl_floating_point *, dbl_floating_point *, dbl_floating_point*, unsigned int *);
    561  1.2     skrll int dbl_frem(dbl_floating_point *, dbl_floating_point *, dbl_floating_point*, unsigned int *);
    562  1.2     skrll int dbl_fsqrt(dbl_floating_point *, dbl_floating_point *, unsigned int *);
    563  1.2     skrll int dbl_fsub(dbl_floating_point *, dbl_floating_point *, dbl_floating_point*, unsigned int *);
    564  1.1  fredette 
    565  1.2     skrll dbl_floating_point dbl_setoverflow(unsigned int);
    566  1.1  fredette 
    567  1.2     skrll int sgl_to_dbl_fcnvff(sgl_floating_point *, dbl_floating_point *, unsigned int *);
    568  1.2     skrll int dbl_to_sgl_fcnvff(dbl_floating_point *, sgl_floating_point *, unsigned int *);
    569  1.1  fredette 
    570  1.2     skrll int dbl_frnd(dbl_floating_point *, dbl_floating_point *, unsigned int *);
    571