]>
Commit | Line | Data |
---|---|---|
2aa62f2b | 1 | /* @(#)e_fmod.c 5.1 93/09/24 */\r |
2 | /*\r | |
3 | * ====================================================\r | |
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\r | |
5 | *\r | |
6 | * Developed at SunPro, a Sun Microsystems, Inc. business.\r | |
7 | * Permission to use, copy, modify, and distribute this\r | |
8 | * software is freely granted, provided that this notice\r | |
9 | * is preserved.\r | |
10 | * ====================================================\r | |
11 | */\r | |
12 | #include <LibConfig.h>\r | |
13 | #include <sys/EfiCdefs.h>\r | |
14 | #if defined(LIBM_SCCS) && !defined(lint)\r | |
15 | __RCSID("$NetBSD: e_fmod.c,v 1.11 2002/05/26 22:01:49 wiz Exp $");\r | |
16 | #endif\r | |
17 | \r | |
18 | /*\r | |
19 | * __ieee754_fmod(x,y)\r | |
20 | * Return x mod y in exact arithmetic\r | |
21 | * Method: shift and subtract\r | |
22 | */\r | |
23 | \r | |
24 | #include "math.h"\r | |
25 | #include "math_private.h"\r | |
26 | \r | |
27 | #if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */\r | |
28 | // unary minus operator applied to unsigned type, result still unsigned\r | |
29 | #pragma warning ( disable : 4146 )\r | |
30 | #endif\r | |
31 | \r | |
32 | static const double one = 1.0, Zero[] = {0.0, -0.0,};\r | |
33 | \r | |
34 | double\r | |
35 | __ieee754_fmod(double x, double y)\r | |
36 | {\r | |
37 | int32_t n,hx,hy,hz,ix,iy,sx,i;\r | |
38 | u_int32_t lx,ly,lz;\r | |
39 | \r | |
40 | EXTRACT_WORDS(hx,lx,x);\r | |
41 | EXTRACT_WORDS(hy,ly,y);\r | |
42 | sx = hx&0x80000000; /* sign of x */\r | |
43 | hx ^=sx; /* |x| */\r | |
44 | hy &= 0x7fffffff; /* |y| */\r | |
45 | \r | |
46 | /* purge off exception values */\r | |
47 | if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */\r | |
48 | ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */\r | |
49 | return (x*y)/(x*y);\r | |
50 | if(hx<=hy) {\r | |
51 | if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */\r | |
52 | if(lx==ly)\r | |
53 | return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/\r | |
54 | }\r | |
55 | \r | |
56 | /* determine ix = ilogb(x) */\r | |
57 | if(hx<0x00100000) { /* subnormal x */\r | |
58 | if(hx==0) {\r | |
59 | for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;\r | |
60 | } else {\r | |
61 | for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;\r | |
62 | }\r | |
63 | } else ix = (hx>>20)-1023;\r | |
64 | \r | |
65 | /* determine iy = ilogb(y) */\r | |
66 | if(hy<0x00100000) { /* subnormal y */\r | |
67 | if(hy==0) {\r | |
68 | for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;\r | |
69 | } else {\r | |
70 | for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;\r | |
71 | }\r | |
72 | } else iy = (hy>>20)-1023;\r | |
73 | \r | |
74 | /* set up {hx,lx}, {hy,ly} and align y to x */\r | |
75 | if(ix >= -1022)\r | |
76 | hx = 0x00100000|(0x000fffff&hx);\r | |
77 | else { /* subnormal x, shift x to normal */\r | |
78 | n = -1022-ix;\r | |
79 | if(n<=31) {\r | |
80 | hx = (hx<<n)|(lx>>(32-n));\r | |
81 | lx <<= n;\r | |
82 | } else {\r | |
83 | hx = lx<<(n-32);\r | |
84 | lx = 0;\r | |
85 | }\r | |
86 | }\r | |
87 | if(iy >= -1022)\r | |
88 | hy = 0x00100000|(0x000fffff&hy);\r | |
89 | else { /* subnormal y, shift y to normal */\r | |
90 | n = -1022-iy;\r | |
91 | if(n<=31) {\r | |
92 | hy = (hy<<n)|(ly>>(32-n));\r | |
93 | ly <<= n;\r | |
94 | } else {\r | |
95 | hy = ly<<(n-32);\r | |
96 | ly = 0;\r | |
97 | }\r | |
98 | }\r | |
99 | \r | |
100 | /* fix point fmod */\r | |
101 | n = ix - iy;\r | |
102 | while(n--) {\r | |
103 | hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;\r | |
104 | if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}\r | |
105 | else {\r | |
106 | if((hz|lz)==0) /* return sign(x)*0 */\r | |
107 | return Zero[(u_int32_t)sx>>31];\r | |
108 | hx = hz+hz+(lz>>31); lx = lz+lz;\r | |
109 | }\r | |
110 | }\r | |
111 | hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;\r | |
112 | if(hz>=0) {hx=hz;lx=lz;}\r | |
113 | \r | |
114 | /* convert back to floating value and restore the sign */\r | |
115 | if((hx|lx)==0) /* return sign(x)*0 */\r | |
116 | return Zero[(u_int32_t)sx>>31];\r | |
117 | while(hx<0x00100000) { /* normalize x */\r | |
118 | hx = hx+hx+(lx>>31); lx = lx+lx;\r | |
119 | iy -= 1;\r | |
120 | }\r | |
121 | if(iy>= -1022) { /* normalize output */\r | |
122 | hx = ((hx-0x00100000)|((iy+1023)<<20));\r | |
123 | INSERT_WORDS(x,hx|sx,lx);\r | |
124 | } else { /* subnormal output */\r | |
125 | n = -1022 - iy;\r | |
126 | if(n<=20) {\r | |
127 | lx = (lx>>n)|((u_int32_t)hx<<(32-n));\r | |
128 | hx >>= n;\r | |
129 | } else if (n<=31) {\r | |
130 | lx = (hx<<(32-n))|(lx>>n); hx = sx;\r | |
131 | } else {\r | |
132 | lx = hx>>(n-32); hx = sx;\r | |
133 | }\r | |
134 | INSERT_WORDS(x,hx|sx,lx);\r | |
135 | x *= one; /* create necessary signal */\r | |
136 | }\r | |
137 | return x; /* exact output */\r | |
138 | }\r |