Rev 737 | Rev 828 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
| Rev 737 | Rev 804 | ||
|---|---|---|---|
| Line 34... | Line 34... | ||
| 34 | * |
34 | * |
| 35 | */ |
35 | */ |
| 36 | float32 mulFloat32(float32 a, float32 b) |
36 | float32 mulFloat32(float32 a, float32 b) |
| 37 | { |
37 | { |
| 38 | float32 result; |
38 | float32 result; |
| 39 | __u64 mant1, mant2; |
39 | __u64 frac1, frac2; |
| 40 | __s32 exp; |
40 | __s32 exp; |
| 41 | 41 | ||
| 42 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
42 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
| 43 | 43 | ||
| 44 | if (isFloat32NaN(a) || isFloat32NaN(b) ) { |
44 | if (isFloat32NaN(a) || isFloat32NaN(b) ) { |
| 45 | /* TODO: fix SigNaNs */ |
45 | /* TODO: fix SigNaNs */ |
| 46 | if (isFloat32SigNaN(a)) { |
46 | if (isFloat32SigNaN(a)) { |
| 47 | result.parts.mantisa = a.parts.mantisa; |
47 | result.parts.fraction = a.parts.fraction; |
| 48 | result.parts.exp = a.parts.exp; |
48 | result.parts.exp = a.parts.exp; |
| 49 | return result; |
49 | return result; |
| 50 | }; |
50 | }; |
| 51 | if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */ |
51 | if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */ |
| 52 | result.parts.mantisa = b.parts.mantisa; |
52 | result.parts.fraction = b.parts.fraction; |
| 53 | result.parts.exp = b.parts.exp; |
53 | result.parts.exp = b.parts.exp; |
| 54 | return result; |
54 | return result; |
| 55 | }; |
55 | }; |
| 56 | /* set NaN as result */ |
56 | /* set NaN as result */ |
| 57 | result.binary = FLOAT32_NAN; |
57 | result.binary = FLOAT32_NAN; |
| Line 62... | Line 62... | ||
| 62 | if (isFloat32Zero(b)) { |
62 | if (isFloat32Zero(b)) { |
| 63 | /* FIXME: zero * infinity */ |
63 | /* FIXME: zero * infinity */ |
| 64 | result.binary = FLOAT32_NAN; |
64 | result.binary = FLOAT32_NAN; |
| 65 | return result; |
65 | return result; |
| 66 | } |
66 | } |
| 67 | result.parts.mantisa = a.parts.mantisa; |
67 | result.parts.fraction = a.parts.fraction; |
| 68 | result.parts.exp = a.parts.exp; |
68 | result.parts.exp = a.parts.exp; |
| 69 | return result; |
69 | return result; |
| 70 | } |
70 | } |
| 71 | 71 | ||
| 72 | if (isFloat32Infinity(b)) { |
72 | if (isFloat32Infinity(b)) { |
| 73 | if (isFloat32Zero(a)) { |
73 | if (isFloat32Zero(a)) { |
| 74 | /* FIXME: zero * infinity */ |
74 | /* FIXME: zero * infinity */ |
| 75 | result.binary = FLOAT32_NAN; |
75 | result.binary = FLOAT32_NAN; |
| 76 | return result; |
76 | return result; |
| 77 | } |
77 | } |
| 78 | result.parts.mantisa = b.parts.mantisa; |
78 | result.parts.fraction = b.parts.fraction; |
| 79 | result.parts.exp = b.parts.exp; |
79 | result.parts.exp = b.parts.exp; |
| 80 | return result; |
80 | return result; |
| 81 | } |
81 | } |
| 82 | 82 | ||
| 83 | /* exp is signed so we can easy detect underflow */ |
83 | /* exp is signed so we can easy detect underflow */ |
| Line 93... | Line 93... | ||
| 93 | }; |
93 | }; |
| 94 | 94 | ||
| 95 | if (exp < 0) { |
95 | if (exp < 0) { |
| 96 | /* FIXME: underflow */ |
96 | /* FIXME: underflow */ |
| 97 | /* return signed zero */ |
97 | /* return signed zero */ |
| 98 | result.parts.mantisa = 0x0; |
98 | result.parts.fraction = 0x0; |
| 99 | result.parts.exp = 0x0; |
99 | result.parts.exp = 0x0; |
| 100 | return result; |
100 | return result; |
| 101 | }; |
101 | }; |
| 102 | 102 | ||
| 103 | mant1 = a.parts.mantisa; |
103 | frac1 = a.parts.fraction; |
| 104 | if (a.parts.exp > 0) { |
104 | if (a.parts.exp > 0) { |
| 105 | mant1 |= FLOAT32_HIDDEN_BIT_MASK; |
105 | frac1 |= FLOAT32_HIDDEN_BIT_MASK; |
| 106 | } else { |
106 | } else { |
| 107 | ++exp; |
107 | ++exp; |
| 108 | }; |
108 | }; |
| 109 | 109 | ||
| 110 | mant2 = b.parts.mantisa; |
110 | frac2 = b.parts.fraction; |
| 111 | 111 | ||
| 112 | if (b.parts.exp > 0) { |
112 | if (b.parts.exp > 0) { |
| 113 | mant2 |= FLOAT32_HIDDEN_BIT_MASK; |
113 | frac2 |= FLOAT32_HIDDEN_BIT_MASK; |
| 114 | } else { |
114 | } else { |
| 115 | ++exp; |
115 | ++exp; |
| 116 | }; |
116 | }; |
| 117 | 117 | ||
| 118 | mant1 <<= 1; /* one bit space for rounding */ |
118 | frac1 <<= 1; /* one bit space for rounding */ |
| 119 | 119 | ||
| 120 | mant1 = mant1 * mant2; |
120 | frac1 = frac1 * frac2; |
| 121 | /* round and return */ |
121 | /* round and return */ |
| 122 | 122 | ||
| 123 | while ((exp < FLOAT32_MAX_EXPONENT) && (mant1 >= ( 1 << (FLOAT32_MANTISA_SIZE + 2)))) { |
123 | while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) { |
| 124 | /* 23 bits of mantisa + one more for hidden bit (all shifted 1 bit left)*/ |
124 | /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/ |
| 125 | ++exp; |
125 | ++exp; |
| 126 | mant1 >>= 1; |
126 | frac1 >>= 1; |
| 127 | }; |
127 | }; |
| 128 | 128 | ||
| 129 | /* rounding */ |
129 | /* rounding */ |
| 130 | //++mant1; /* FIXME: not works - without it is ok */ |
130 | /* ++frac1; FIXME: not works - without it is ok */ |
| 131 | mant1 >>= 1; /* shift off rounding space */ |
131 | frac1 >>= 1; /* shift off rounding space */ |
| 132 | 132 | ||
| 133 | if ((exp < FLOAT32_MAX_EXPONENT) && (mant1 >= (1 << (FLOAT32_MANTISA_SIZE + 1)))) { |
133 | if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) { |
| 134 | ++exp; |
134 | ++exp; |
| 135 | mant1 >>= 1; |
135 | frac1 >>= 1; |
| 136 | }; |
136 | }; |
| 137 | 137 | ||
| 138 | if (exp >= FLOAT32_MAX_EXPONENT ) { |
138 | if (exp >= FLOAT32_MAX_EXPONENT ) { |
| 139 | /* TODO: fix overflow */ |
139 | /* TODO: fix overflow */ |
| 140 | /* return infinity*/ |
140 | /* return infinity*/ |
| 141 | result.parts.exp = FLOAT32_MAX_EXPONENT; |
141 | result.parts.exp = FLOAT32_MAX_EXPONENT; |
| 142 | result.parts.mantisa = 0x0; |
142 | result.parts.fraction = 0x0; |
| 143 | return result; |
143 | return result; |
| 144 | } |
144 | } |
| 145 | 145 | ||
| 146 | exp -= FLOAT32_MANTISA_SIZE; |
146 | exp -= FLOAT32_FRACTION_SIZE; |
| 147 | 147 | ||
| 148 | if (exp <= FLOAT32_MANTISA_SIZE) { |
148 | if (exp <= FLOAT32_FRACTION_SIZE) { |
| 149 | /* denormalized number */ |
149 | /* denormalized number */ |
| 150 | mant1 >>= 1; /* denormalize */ |
150 | frac1 >>= 1; /* denormalize */ |
| 151 | while ((mant1 > 0) && (exp < 0)) { |
151 | while ((frac1 > 0) && (exp < 0)) { |
| 152 | mant1 >>= 1; |
152 | frac1 >>= 1; |
| 153 | ++exp; |
153 | ++exp; |
| 154 | }; |
154 | }; |
| 155 | if (mant1 == 0) { |
155 | if (frac1 == 0) { |
| 156 | /* FIXME : underflow */ |
156 | /* FIXME : underflow */ |
| 157 | result.parts.exp = 0; |
157 | result.parts.exp = 0; |
| 158 | result.parts.mantisa = 0; |
158 | result.parts.fraction = 0; |
| 159 | return result; |
159 | return result; |
| 160 | }; |
160 | }; |
| 161 | }; |
161 | }; |
| 162 | result.parts.exp = exp; |
162 | result.parts.exp = exp; |
| 163 | result.parts.mantisa = mant1 & ( (1 << FLOAT32_MANTISA_SIZE) - 1); |
163 | result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1); |
| 164 | 164 | ||
| 165 | return result; |
165 | return result; |
| 166 | 166 | ||
| 167 | } |
167 | } |
| 168 | 168 | ||
| Line 170... | Line 170... | ||
| 170 | * |
170 | * |
| 171 | */ |
171 | */ |
| 172 | float64 mulFloat64(float64 a, float64 b) |
172 | float64 mulFloat64(float64 a, float64 b) |
| 173 | { |
173 | { |
| 174 | float64 result; |
174 | float64 result; |
| 175 | __u64 mant1, mant2; |
175 | __u64 frac1, frac2; |
| 176 | __s32 exp; |
176 | __s32 exp; |
| 177 | 177 | ||
| 178 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
178 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
| 179 | 179 | ||
| 180 | if (isFloat64NaN(a) || isFloat64NaN(b) ) { |
180 | if (isFloat64NaN(a) || isFloat64NaN(b) ) { |
| 181 | /* TODO: fix SigNaNs */ |
181 | /* TODO: fix SigNaNs */ |
| 182 | if (isFloat64SigNaN(a)) { |
182 | if (isFloat64SigNaN(a)) { |
| 183 | result.parts.mantisa = a.parts.mantisa; |
183 | result.parts.fraction = a.parts.fraction; |
| 184 | result.parts.exp = a.parts.exp; |
184 | result.parts.exp = a.parts.exp; |
| 185 | return result; |
185 | return result; |
| 186 | }; |
186 | }; |
| 187 | if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */ |
187 | if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */ |
| 188 | result.parts.mantisa = b.parts.mantisa; |
188 | result.parts.fraction = b.parts.fraction; |
| 189 | result.parts.exp = b.parts.exp; |
189 | result.parts.exp = b.parts.exp; |
| 190 | return result; |
190 | return result; |
| 191 | }; |
191 | }; |
| 192 | /* set NaN as result */ |
192 | /* set NaN as result */ |
| 193 | result.binary = FLOAT64_NAN; |
193 | result.binary = FLOAT64_NAN; |
| Line 198... | Line 198... | ||
| 198 | if (isFloat64Zero(b)) { |
198 | if (isFloat64Zero(b)) { |
| 199 | /* FIXME: zero * infinity */ |
199 | /* FIXME: zero * infinity */ |
| 200 | result.binary = FLOAT64_NAN; |
200 | result.binary = FLOAT64_NAN; |
| 201 | return result; |
201 | return result; |
| 202 | } |
202 | } |
| 203 | result.parts.mantisa = a.parts.mantisa; |
203 | result.parts.fraction = a.parts.fraction; |
| 204 | result.parts.exp = a.parts.exp; |
204 | result.parts.exp = a.parts.exp; |
| 205 | return result; |
205 | return result; |
| 206 | } |
206 | } |
| 207 | 207 | ||
| 208 | if (isFloat64Infinity(b)) { |
208 | if (isFloat64Infinity(b)) { |
| 209 | if (isFloat64Zero(a)) { |
209 | if (isFloat64Zero(a)) { |
| 210 | /* FIXME: zero * infinity */ |
210 | /* FIXME: zero * infinity */ |
| 211 | result.binary = FLOAT64_NAN; |
211 | result.binary = FLOAT64_NAN; |
| 212 | return result; |
212 | return result; |
| 213 | } |
213 | } |
| 214 | result.parts.mantisa = b.parts.mantisa; |
214 | result.parts.fraction = b.parts.fraction; |
| 215 | result.parts.exp = b.parts.exp; |
215 | result.parts.exp = b.parts.exp; |
| 216 | return result; |
216 | return result; |
| 217 | } |
217 | } |
| 218 | 218 | ||
| 219 | /* exp is signed so we can easy detect underflow */ |
219 | /* exp is signed so we can easy detect underflow */ |
| Line 229... | Line 229... | ||
| 229 | }; |
229 | }; |
| 230 | 230 | ||
| 231 | if (exp < 0) { |
231 | if (exp < 0) { |
| 232 | /* FIXME: underflow */ |
232 | /* FIXME: underflow */ |
| 233 | /* return signed zero */ |
233 | /* return signed zero */ |
| 234 | result.parts.mantisa = 0x0; |
234 | result.parts.fraction = 0x0; |
| 235 | result.parts.exp = 0x0; |
235 | result.parts.exp = 0x0; |
| 236 | return result; |
236 | return result; |
| 237 | }; |
237 | }; |
| 238 | 238 | ||
| 239 | mant1 = a.parts.mantisa; |
239 | frac1 = a.parts.fraction; |
| 240 | if (a.parts.exp > 0) { |
240 | if (a.parts.exp > 0) { |
| 241 | mant1 |= FLOAT64_HIDDEN_BIT_MASK; |
241 | frac1 |= FLOAT64_HIDDEN_BIT_MASK; |
| 242 | } else { |
242 | } else { |
| 243 | ++exp; |
243 | ++exp; |
| 244 | }; |
244 | }; |
| 245 | 245 | ||
| 246 | mant2 = b.parts.mantisa; |
246 | frac2 = b.parts.fraction; |
| 247 | 247 | ||
| 248 | if (b.parts.exp > 0) { |
248 | if (b.parts.exp > 0) { |
| 249 | mant2 |= FLOAT64_HIDDEN_BIT_MASK; |
249 | frac2 |= FLOAT64_HIDDEN_BIT_MASK; |
| 250 | } else { |
250 | } else { |
| 251 | ++exp; |
251 | ++exp; |
| 252 | }; |
252 | }; |
| 253 | 253 | ||
| 254 | mant1 <<= 1; /* one bit space for rounding */ |
254 | frac1 <<= 1; /* one bit space for rounding */ |
| 255 | 255 | ||
| 256 | mul64integers(mant1, mant2, &mant1, &mant2); |
256 | mul64integers(frac1, frac2, &frac1, &frac2); |
| 257 | 257 | ||
| 258 | /* round and return */ |
258 | /* round and return */ |
| 259 | /* FIXME: ugly soulution is to shift whole mant2 >> as in 32bit version |
259 | /* FIXME: ugly soulution is to shift whole frac2 >> as in 32bit version |
| 260 | * Here is is more slower because we have to shift two numbers with carry |
260 | * Here is is more slower because we have to shift two numbers with carry |
| 261 | * Better is find first nonzero bit and make only one shift |
261 | * Better is find first nonzero bit and make only one shift |
| 262 | * Third version is to shift both numbers a bit to right and result will be then |
262 | * Third version is to shift both numbers a bit to right and result will be then |
| 263 | * placed in higher part of result. Then lower part will be good only for rounding. |
263 | * placed in higher part of result. Then lower part will be good only for rounding. |
| 264 | */ |
264 | */ |
| 265 | 265 | ||
| 266 | while ((exp < FLOAT64_MAX_EXPONENT) && (mant2 > 0 )) { |
266 | while ((exp < FLOAT64_MAX_EXPONENT) && (frac2 > 0 )) { |
| 267 | mant1 >>= 1; |
267 | frac1 >>= 1; |
| 268 | mant1 &= ((mant2 & 0x1) << 63); |
268 | frac1 &= ((frac2 & 0x1) << 63); |
| 269 | mant2 >>= 1; |
269 | frac2 >>= 1; |
| 270 | ++exp; |
270 | ++exp; |
| 271 | } |
271 | } |
| 272 | 272 | ||
| 273 | while ((exp < FLOAT64_MAX_EXPONENT) && (mant1 >= ( (__u64)1 << (FLOAT64_MANTISA_SIZE + 2)))) { |
273 | while ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ( (__u64)1 << (FLOAT64_FRACTION_SIZE + 2)))) { |
| 274 | ++exp; |
274 | ++exp; |
| 275 | mant1 >>= 1; |
275 | frac1 >>= 1; |
| 276 | }; |
276 | }; |
| 277 | 277 | ||
| 278 | /* rounding */ |
278 | /* rounding */ |
| 279 | //++mant1; /* FIXME: not works - without it is ok */ |
279 | /* ++frac1; FIXME: not works - without it is ok */ |
| 280 | mant1 >>= 1; /* shift off rounding space */ |
280 | frac1 >>= 1; /* shift off rounding space */ |
| 281 | 281 | ||
| 282 | if ((exp < FLOAT64_MAX_EXPONENT) && (mant1 >= ((__u64)1 << (FLOAT64_MANTISA_SIZE + 1)))) { |
282 | if ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ((__u64)1 << (FLOAT64_FRACTION_SIZE + 1)))) { |
| 283 | ++exp; |
283 | ++exp; |
| 284 | mant1 >>= 1; |
284 | frac1 >>= 1; |
| 285 | }; |
285 | }; |
| 286 | 286 | ||
| 287 | if (exp >= FLOAT64_MAX_EXPONENT ) { |
287 | if (exp >= FLOAT64_MAX_EXPONENT ) { |
| 288 | /* TODO: fix overflow */ |
288 | /* TODO: fix overflow */ |
| 289 | /* return infinity*/ |
289 | /* return infinity*/ |
| 290 | result.parts.exp = FLOAT64_MAX_EXPONENT; |
290 | result.parts.exp = FLOAT64_MAX_EXPONENT; |
| 291 | result.parts.mantisa = 0x0; |
291 | result.parts.fraction = 0x0; |
| 292 | return result; |
292 | return result; |
| 293 | } |
293 | } |
| 294 | 294 | ||
| 295 | exp -= FLOAT64_MANTISA_SIZE; |
295 | exp -= FLOAT64_FRACTION_SIZE; |
| 296 | 296 | ||
| 297 | if (exp <= FLOAT64_MANTISA_SIZE) { |
297 | if (exp <= FLOAT64_FRACTION_SIZE) { |
| 298 | /* denormalized number */ |
298 | /* denormalized number */ |
| 299 | mant1 >>= 1; /* denormalize */ |
299 | frac1 >>= 1; /* denormalize */ |
| 300 | while ((mant1 > 0) && (exp < 0)) { |
300 | while ((frac1 > 0) && (exp < 0)) { |
| 301 | mant1 >>= 1; |
301 | frac1 >>= 1; |
| 302 | ++exp; |
302 | ++exp; |
| 303 | }; |
303 | }; |
| 304 | if (mant1 == 0) { |
304 | if (frac1 == 0) { |
| 305 | /* FIXME : underflow */ |
305 | /* FIXME : underflow */ |
| 306 | result.parts.exp = 0; |
306 | result.parts.exp = 0; |
| 307 | result.parts.mantisa = 0; |
307 | result.parts.fraction = 0; |
| 308 | return result; |
308 | return result; |
| 309 | }; |
309 | }; |
| 310 | }; |
310 | }; |
| 311 | result.parts.exp = exp; |
311 | result.parts.exp = exp; |
| 312 | result.parts.mantisa = mant1 & ( ((__u64)1 << FLOAT64_MANTISA_SIZE) - 1); |
312 | result.parts.fraction = frac1 & ( ((__u64)1 << FLOAT64_FRACTION_SIZE) - 1); |
| 313 | 313 | ||
| 314 | return result; |
314 | return result; |
| 315 | 315 | ||
| 316 | } |
316 | } |
| 317 | 317 | ||
| Line 336... | Line 336... | ||
| 336 | middle1 = a * blow; |
336 | middle1 = a * blow; |
| 337 | middle2 = alow * b; |
337 | middle2 = alow * b; |
| 338 | high = a * b; |
338 | high = a * b; |
| 339 | 339 | ||
| 340 | middle1 += middle2; |
340 | middle1 += middle2; |
| 341 | high += ((__u64)(middle1 < middle2) << 32) + middle1>>32; |
341 | high += ((__u64)(middle1 < middle2) << 32) + (middle1 >> 32); |
| 342 | middle1 << 32; |
342 | middle1 <<= 32; |
| 343 | low += middle1; |
343 | low += middle1; |
| 344 | high += (low < middle1); |
344 | high += (low < middle1); |
| 345 | *lo = low; |
345 | *lo = low; |
| 346 | *hi = high; |
346 | *hi = high; |
| 347 | return; |
347 | return; |