]>
git.proxmox.com Git - mirror_edk2.git/blob - StdLib/LibC/Math/e_fmod.c
1 /* @(#)e_fmod.c 5.1 93/09/24 */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
10 * ====================================================
12 #include <LibConfig.h>
13 #include <sys/EfiCdefs.h>
14 #if defined(LIBM_SCCS) && !defined(lint)
15 __RCSID("$NetBSD: e_fmod.c,v 1.11 2002/05/26 22:01:49 wiz Exp $");
20 * Return x mod y in exact arithmetic
21 * Method: shift and subtract
25 #include "math_private.h"
27 #if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */
28 // unary minus operator applied to unsigned type, result still unsigned
29 #pragma warning ( disable : 4146 )
32 static const double one
= 1.0, Zero
[] = {0.0, -0.0,};
35 __ieee754_fmod(double x
, double y
)
37 int32_t n
,hx
,hy
,hz
,ix
,iy
,sx
,i
;
40 EXTRACT_WORDS(hx
,lx
,x
);
41 EXTRACT_WORDS(hy
,ly
,y
);
42 sx
= hx
&0x80000000; /* sign of x */
44 hy
&= 0x7fffffff; /* |y| */
46 /* purge off exception values */
47 if((hy
|ly
)==0||(hx
>=0x7ff00000)|| /* y=0,or x not finite */
48 ((hy
|((ly
|-ly
)>>31))>0x7ff00000)) /* or y is NaN */
51 if((hx
<hy
)||(lx
<ly
)) return x
; /* |x|<|y| return x */
53 return Zero
[(u_int32_t
)sx
>>31]; /* |x|=|y| return x*0*/
56 /* determine ix = ilogb(x) */
57 if(hx
<0x00100000) { /* subnormal x */
59 for (ix
= -1043, i
=lx
; i
>0; i
<<=1) ix
-=1;
61 for (ix
= -1022,i
=(hx
<<11); i
>0; i
<<=1) ix
-=1;
63 } else ix
= (hx
>>20)-1023;
65 /* determine iy = ilogb(y) */
66 if(hy
<0x00100000) { /* subnormal y */
68 for (iy
= -1043, i
=ly
; i
>0; i
<<=1) iy
-=1;
70 for (iy
= -1022,i
=(hy
<<11); i
>0; i
<<=1) iy
-=1;
72 } else iy
= (hy
>>20)-1023;
74 /* set up {hx,lx}, {hy,ly} and align y to x */
76 hx
= 0x00100000|(0x000fffff&hx
);
77 else { /* subnormal x, shift x to normal */
80 hx
= (hx
<<n
)|(lx
>>(32-n
));
88 hy
= 0x00100000|(0x000fffff&hy
);
89 else { /* subnormal y, shift y to normal */
92 hy
= (hy
<<n
)|(ly
>>(32-n
));
103 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
104 if(hz
<0){hx
= hx
+hx
+(lx
>>31); lx
= lx
+lx
;}
106 if((hz
|lz
)==0) /* return sign(x)*0 */
107 return Zero
[(u_int32_t
)sx
>>31];
108 hx
= hz
+hz
+(lz
>>31); lx
= lz
+lz
;
111 hz
=hx
-hy
;lz
=lx
-ly
; if(lx
<ly
) hz
-= 1;
112 if(hz
>=0) {hx
=hz
;lx
=lz
;}
114 /* convert back to floating value and restore the sign */
115 if((hx
|lx
)==0) /* return sign(x)*0 */
116 return Zero
[(u_int32_t
)sx
>>31];
117 while(hx
<0x00100000) { /* normalize x */
118 hx
= hx
+hx
+(lx
>>31); lx
= lx
+lx
;
121 if(iy
>= -1022) { /* normalize output */
122 hx
= ((hx
-0x00100000)|((iy
+1023)<<20));
123 INSERT_WORDS(x
,hx
|sx
,lx
);
124 } else { /* subnormal output */
127 lx
= (lx
>>n
)|((u_int32_t
)hx
<<(32-n
));
130 lx
= (hx
<<(32-n
))|(lx
>>n
); hx
= sx
;
132 lx
= hx
>>(n
-32); hx
= sx
;
134 INSERT_WORDS(x
,hx
|sx
,lx
);
135 x
*= one
; /* create necessary signal */
137 return x
; /* exact output */