Rev 737 | Rev 829 | Go to most recent revision | Only display areas with differences | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 737 | Rev 804 | ||
---|---|---|---|
1 | /* |
1 | /* |
2 | * Copyright (C) 2005 Josef Cejka |
2 | * Copyright (C) 2005 Josef Cejka |
3 | * All rights reserved. |
3 | * All rights reserved. |
4 | * |
4 | * |
5 | * Redistribution and use in source and binary forms, with or without |
5 | * Redistribution and use in source and binary forms, with or without |
6 | * modification, are permitted provided that the following conditions |
6 | * modification, are permitted provided that the following conditions |
7 | * are met: |
7 | * are met: |
8 | * |
8 | * |
9 | * - Redistributions of source code must retain the above copyright |
9 | * - Redistributions of source code must retain the above copyright |
10 | * notice, this list of conditions and the following disclaimer. |
10 | * notice, this list of conditions and the following disclaimer. |
11 | * - Redistributions in binary form must reproduce the above copyright |
11 | * - Redistributions in binary form must reproduce the above copyright |
12 | * notice, this list of conditions and the following disclaimer in the |
12 | * notice, this list of conditions and the following disclaimer in the |
13 | * documentation and/or other materials provided with the distribution. |
13 | * documentation and/or other materials provided with the distribution. |
14 | * - The name of the author may not be used to endorse or promote products |
14 | * - The name of the author may not be used to endorse or promote products |
15 | * derived from this software without specific prior written permission. |
15 | * derived from this software without specific prior written permission. |
16 | * |
16 | * |
17 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
17 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
18 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
18 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
19 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. |
19 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. |
20 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, |
20 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, |
21 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT |
21 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT |
22 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
22 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
23 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
23 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
24 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
24 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
25 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF |
25 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF |
26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
27 | */ |
27 | */ |
28 | 28 | ||
29 | #include<sftypes.h> |
29 | #include<sftypes.h> |
30 | #include<mul.h> |
30 | #include<mul.h> |
31 | #include<comparison.h> |
31 | #include<comparison.h> |
32 | 32 | ||
33 | /** Multiply two 32 bit float numbers |
33 | /** Multiply two 32 bit float numbers |
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; |
58 | return result; |
58 | return result; |
59 | }; |
59 | }; |
60 | 60 | ||
61 | if (isFloat32Infinity(a)) { |
61 | if (isFloat32Infinity(a)) { |
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 */ |
84 | exp = a.parts.exp + b.parts.exp; |
84 | exp = a.parts.exp + b.parts.exp; |
85 | exp -= FLOAT32_BIAS; |
85 | exp -= FLOAT32_BIAS; |
86 | 86 | ||
87 | if (exp >= FLOAT32_MAX_EXPONENT) { |
87 | if (exp >= FLOAT32_MAX_EXPONENT) { |
88 | /* FIXME: overflow */ |
88 | /* FIXME: overflow */ |
89 | /* set infinity as result */ |
89 | /* set infinity as result */ |
90 | result.binary = FLOAT32_INF; |
90 | result.binary = FLOAT32_INF; |
91 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
91 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
92 | return result; |
92 | return result; |
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 | ||
169 | /** Multiply two 64 bit float numbers |
169 | /** Multiply two 64 bit float numbers |
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; |
194 | return result; |
194 | return result; |
195 | }; |
195 | }; |
196 | 196 | ||
197 | if (isFloat64Infinity(a)) { |
197 | if (isFloat64Infinity(a)) { |
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 */ |
220 | exp = a.parts.exp + b.parts.exp; |
220 | exp = a.parts.exp + b.parts.exp; |
221 | exp -= FLOAT64_BIAS; |
221 | exp -= FLOAT64_BIAS; |
222 | 222 | ||
223 | if (exp >= FLOAT64_MAX_EXPONENT) { |
223 | if (exp >= FLOAT64_MAX_EXPONENT) { |
224 | /* FIXME: overflow */ |
224 | /* FIXME: overflow */ |
225 | /* set infinity as result */ |
225 | /* set infinity as result */ |
226 | result.binary = FLOAT64_INF; |
226 | result.binary = FLOAT64_INF; |
227 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
227 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
228 | return result; |
228 | return result; |
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 | ||
318 | /** Multiply two 64 bit numbers and return result in two parts |
318 | /** Multiply two 64 bit numbers and return result in two parts |
319 | * @param a first operand |
319 | * @param a first operand |
320 | * @param b second operand |
320 | * @param b second operand |
321 | * @param lo lower part from result |
321 | * @param lo lower part from result |
322 | * @param hi higher part of result |
322 | * @param hi higher part of result |
323 | */ |
323 | */ |
324 | void mul64integers(__u64 a,__u64 b, __u64 *lo, __u64 *hi) |
324 | void mul64integers(__u64 a,__u64 b, __u64 *lo, __u64 *hi) |
325 | { |
325 | { |
326 | __u64 low, high, middle1, middle2; |
326 | __u64 low, high, middle1, middle2; |
327 | __u32 alow, blow; |
327 | __u32 alow, blow; |
328 | 328 | ||
329 | alow = a & 0xFFFFFFFF; |
329 | alow = a & 0xFFFFFFFF; |
330 | blow = b & 0xFFFFFFFF; |
330 | blow = b & 0xFFFFFFFF; |
331 | 331 | ||
332 | a <<= 32; |
332 | a <<= 32; |
333 | b <<= 32; |
333 | b <<= 32; |
334 | 334 | ||
335 | low = (__u64)alow * blow; |
335 | low = (__u64)alow * blow; |
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; |
348 | } |
348 | } |
349 | 349 | ||
350 | 350 | ||
351 | 351 |