source: mainline/uspace/lib/softfloat/mul.c@ 8f1dd3b

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since 8f1dd3b was 2416085, checked in by Martin Decky <martin@…>, 13 years ago

softfloat: move sources to a single directory
implement more ARM EABI bindings
  • Property mode set to 100644
File size: 9.2 KB
Line 
1/*
2 * Copyright (c) 2005 Josef Cejka
3 * Copyright (c) 2011 Petr Koupy
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * - Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * - Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * - The name of the author may not be used to endorse or promote products
16 * derived from this software without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
19 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
20 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
22 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
23 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
25 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
26 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
27 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30/** @addtogroup softfloat
31 * @{
32 */
33/** @file Multiplication functions.
34 */
35
36#include "sftypes.h"
37#include "mul.h"
38#include "comparison.h"
39#include "common.h"
40
41/** Multiply two single-precision floats.
42 *
43 * @param a First input operand.
44 * @param b Second input operand.
45 *
46 * @return Result of multiplication.
47 *
48 */
49float32 mul_float32(float32 a, float32 b)
50{
51 float32 result;
52 uint64_t frac1, frac2;
53 int32_t exp;
54
55 result.parts.sign = a.parts.sign ^ b.parts.sign;
56
57 if (is_float32_nan(a) || is_float32_nan(b)) {
58 /* TODO: fix SigNaNs */
59 if (is_float32_signan(a)) {
60 result.parts.fraction = a.parts.fraction;
61 result.parts.exp = a.parts.exp;
62 return result;
63 }
64 if (is_float32_signan(b)) { /* TODO: fix SigNaN */
65 result.parts.fraction = b.parts.fraction;
66 result.parts.exp = b.parts.exp;
67 return result;
68 }
69 /* set NaN as result */
70 result.bin = FLOAT32_NAN;
71 return result;
72 }
73
74 if (is_float32_infinity(a)) {
75 if (is_float32_zero(b)) {
76 /* FIXME: zero * infinity */
77 result.bin = FLOAT32_NAN;
78 return result;
79 }
80 result.parts.fraction = a.parts.fraction;
81 result.parts.exp = a.parts.exp;
82 return result;
83 }
84
85 if (is_float32_infinity(b)) {
86 if (is_float32_zero(a)) {
87 /* FIXME: zero * infinity */
88 result.bin = FLOAT32_NAN;
89 return result;
90 }
91 result.parts.fraction = b.parts.fraction;
92 result.parts.exp = b.parts.exp;
93 return result;
94 }
95
96 /* exp is signed so we can easy detect underflow */
97 exp = a.parts.exp + b.parts.exp;
98 exp -= FLOAT32_BIAS;
99
100 if (exp >= FLOAT32_MAX_EXPONENT) {
101 /* FIXME: overflow */
102 /* set infinity as result */
103 result.bin = FLOAT32_INF;
104 result.parts.sign = a.parts.sign ^ b.parts.sign;
105 return result;
106 }
107
108 if (exp < 0) {
109 /* FIXME: underflow */
110 /* return signed zero */
111 result.parts.fraction = 0x0;
112 result.parts.exp = 0x0;
113 return result;
114 }
115
116 frac1 = a.parts.fraction;
117 if (a.parts.exp > 0) {
118 frac1 |= FLOAT32_HIDDEN_BIT_MASK;
119 } else {
120 ++exp;
121 }
122
123 frac2 = b.parts.fraction;
124
125 if (b.parts.exp > 0) {
126 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
127 } else {
128 ++exp;
129 }
130
131 frac1 <<= 1; /* one bit space for rounding */
132
133 frac1 = frac1 * frac2;
134
135 /* round and return */
136 while ((exp < FLOAT32_MAX_EXPONENT) &&
137 (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 2)))) {
138 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left) */
139 ++exp;
140 frac1 >>= 1;
141 }
142
143 /* rounding */
144 /* ++frac1; FIXME: not works - without it is ok */
145 frac1 >>= 1; /* shift off rounding space */
146
147 if ((exp < FLOAT32_MAX_EXPONENT) &&
148 (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
149 ++exp;
150 frac1 >>= 1;
151 }
152
153 if (exp >= FLOAT32_MAX_EXPONENT) {
154 /* TODO: fix overflow */
155 /* return infinity*/
156 result.parts.exp = FLOAT32_MAX_EXPONENT;
157 result.parts.fraction = 0x0;
158 return result;
159 }
160
161 exp -= FLOAT32_FRACTION_SIZE;
162
163 if (exp <= FLOAT32_FRACTION_SIZE) {
164 /* denormalized number */
165 frac1 >>= 1; /* denormalize */
166 while ((frac1 > 0) && (exp < 0)) {
167 frac1 >>= 1;
168 ++exp;
169 }
170 if (frac1 == 0) {
171 /* FIXME : underflow */
172 result.parts.exp = 0;
173 result.parts.fraction = 0;
174 return result;
175 }
176 }
177 result.parts.exp = exp;
178 result.parts.fraction = frac1 & ((1 << FLOAT32_FRACTION_SIZE) - 1);
179
180 return result;
181}
182
183/** Multiply two double-precision floats.
184 *
185 * @param a First input operand.
186 * @param b Second input operand.
187 *
188 * @return Result of multiplication.
189 *
190 */
191float64 mul_float64(float64 a, float64 b)
192{
193 float64 result;
194 uint64_t frac1, frac2;
195 int32_t exp;
196
197 result.parts.sign = a.parts.sign ^ b.parts.sign;
198
199 if (is_float64_nan(a) || is_float64_nan(b)) {
200 /* TODO: fix SigNaNs */
201 if (is_float64_signan(a)) {
202 result.parts.fraction = a.parts.fraction;
203 result.parts.exp = a.parts.exp;
204 return result;
205 }
206 if (is_float64_signan(b)) { /* TODO: fix SigNaN */
207 result.parts.fraction = b.parts.fraction;
208 result.parts.exp = b.parts.exp;
209 return result;
210 }
211 /* set NaN as result */
212 result.bin = FLOAT64_NAN;
213 return result;
214 }
215
216 if (is_float64_infinity(a)) {
217 if (is_float64_zero(b)) {
218 /* FIXME: zero * infinity */
219 result.bin = FLOAT64_NAN;
220 return result;
221 }
222 result.parts.fraction = a.parts.fraction;
223 result.parts.exp = a.parts.exp;
224 return result;
225 }
226
227 if (is_float64_infinity(b)) {
228 if (is_float64_zero(a)) {
229 /* FIXME: zero * infinity */
230 result.bin = FLOAT64_NAN;
231 return result;
232 }
233 result.parts.fraction = b.parts.fraction;
234 result.parts.exp = b.parts.exp;
235 return result;
236 }
237
238 /* exp is signed so we can easy detect underflow */
239 exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
240
241 frac1 = a.parts.fraction;
242
243 if (a.parts.exp > 0) {
244 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
245 } else {
246 ++exp;
247 }
248
249 frac2 = b.parts.fraction;
250
251 if (b.parts.exp > 0) {
252 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
253 } else {
254 ++exp;
255 }
256
257 frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
258 frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
259
260 mul64(frac1, frac2, &frac1, &frac2);
261
262 frac1 |= (frac2 != 0);
263 if (frac1 & (0x1ll << 62)) {
264 frac1 <<= 1;
265 exp--;
266 }
267
268 result = finish_float64(exp, frac1, result.parts.sign);
269 return result;
270}
271
272/** Multiply two quadruple-precision floats.
273 *
274 * @param a First input operand.
275 * @param b Second input operand.
276 *
277 * @return Result of multiplication.
278 *
279 */
280float128 mul_float128(float128 a, float128 b)
281{
282 float128 result;
283 uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
284 int32_t exp;
285
286 result.parts.sign = a.parts.sign ^ b.parts.sign;
287
288 if (is_float128_nan(a) || is_float128_nan(b)) {
289 /* TODO: fix SigNaNs */
290 if (is_float128_signan(a)) {
291 result.parts.frac_hi = a.parts.frac_hi;
292 result.parts.frac_lo = a.parts.frac_lo;
293 result.parts.exp = a.parts.exp;
294 return result;
295 }
296 if (is_float128_signan(b)) { /* TODO: fix SigNaN */
297 result.parts.frac_hi = b.parts.frac_hi;
298 result.parts.frac_lo = b.parts.frac_lo;
299 result.parts.exp = b.parts.exp;
300 return result;
301 }
302 /* set NaN as result */
303 result.bin.hi = FLOAT128_NAN_HI;
304 result.bin.lo = FLOAT128_NAN_LO;
305 return result;
306 }
307
308 if (is_float128_infinity(a)) {
309 if (is_float128_zero(b)) {
310 /* FIXME: zero * infinity */
311 result.bin.hi = FLOAT128_NAN_HI;
312 result.bin.lo = FLOAT128_NAN_LO;
313 return result;
314 }
315 result.parts.frac_hi = a.parts.frac_hi;
316 result.parts.frac_lo = a.parts.frac_lo;
317 result.parts.exp = a.parts.exp;
318 return result;
319 }
320
321 if (is_float128_infinity(b)) {
322 if (is_float128_zero(a)) {
323 /* FIXME: zero * infinity */
324 result.bin.hi = FLOAT128_NAN_HI;
325 result.bin.lo = FLOAT128_NAN_LO;
326 return result;
327 }
328 result.parts.frac_hi = b.parts.frac_hi;
329 result.parts.frac_lo = b.parts.frac_lo;
330 result.parts.exp = b.parts.exp;
331 return result;
332 }
333
334 /* exp is signed so we can easy detect underflow */
335 exp = a.parts.exp + b.parts.exp - FLOAT128_BIAS - 1;
336
337 frac1_hi = a.parts.frac_hi;
338 frac1_lo = a.parts.frac_lo;
339
340 if (a.parts.exp > 0) {
341 or128(frac1_hi, frac1_lo,
342 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
343 &frac1_hi, &frac1_lo);
344 } else {
345 ++exp;
346 }
347
348 frac2_hi = b.parts.frac_hi;
349 frac2_lo = b.parts.frac_lo;
350
351 if (b.parts.exp > 0) {
352 or128(frac2_hi, frac2_lo,
353 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
354 &frac2_hi, &frac2_lo);
355 } else {
356 ++exp;
357 }
358
359 lshift128(frac2_hi, frac2_lo,
360 128 - FLOAT128_FRACTION_SIZE, &frac2_hi, &frac2_lo);
361
362 tmp_hi = frac1_hi;
363 tmp_lo = frac1_lo;
364 mul128(frac1_hi, frac1_lo, frac2_hi, frac2_lo,
365 &frac1_hi, &frac1_lo, &frac2_hi, &frac2_lo);
366 add128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo);
367 frac2_hi |= (frac2_lo != 0x0ll);
368
369 if ((FLOAT128_HIDDEN_BIT_MASK_HI << 1) <= frac1_hi) {
370 frac2_hi >>= 1;
371 if (frac1_lo & 0x1ll) {
372 frac2_hi |= (0x1ull < 64);
373 }
374 rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
375 ++exp;
376 }
377
378 result = finish_float128(exp, frac1_hi, frac1_lo, result.parts.sign, frac2_hi);
379 return result;
380}
381
382/** @}
383 */
Note: See TracBrowser for help on using the repository browser.