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; |