source: mainline/uspace/lib/softfloat/mul.c@ 349b0f7

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since 349b0f7 was a35b458, checked in by Jiří Zárevúcky <zarevucky.jiri@…>, 8 years ago

style: Remove trailing whitespace on _all_ lines, including empty ones, for particular file types.

Command used: tools/srepl '\s\+$' '' -- *.c *.h *.py *.sh *.s *.S *.ag

Currently, whitespace on empty lines is very inconsistent.
There are two basic choices: Either remove the whitespace, or keep empty lines
indented to the level of surrounding code. The former is AFAICT more common,
and also much easier to do automatically.

Alternatively, we could write script for automatic indentation, and use that
instead. However, if such a script exists, it's possible to use the indented
style locally, by having the editor apply relevant conversions on load/save,
without affecting remote repository. IMO, it makes more sense to adopt
the simpler rule.
  • Property mode set to 100644
File size: 10.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 "mul.h"
37#include "comparison.h"
38#include "common.h"
39
40/** Multiply two single-precision floats.
41 *
42 * @param a First input operand.
43 * @param b Second input operand.
44 *
45 * @return Result of multiplication.
46 *
47 */
48float32 mul_float32(float32 a, float32 b)
49{
50 float32 result;
51 uint64_t frac1, frac2;
52 int32_t exp;
53
54 result.parts.sign = a.parts.sign ^ b.parts.sign;
55
56 if (is_float32_nan(a) || is_float32_nan(b)) {
57 /* TODO: fix SigNaNs */
58 if (is_float32_signan(a)) {
59 result.parts.fraction = a.parts.fraction;
60 result.parts.exp = a.parts.exp;
61 return result;
62 }
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
70 /* set NaN as result */
71 result.bin = FLOAT32_NAN;
72 return result;
73 }
74
75 if (is_float32_infinity(a)) {
76 if (is_float32_zero(b)) {
77 /* FIXME: zero * infinity */
78 result.bin = FLOAT32_NAN;
79 return result;
80 }
81
82 result.parts.fraction = a.parts.fraction;
83 result.parts.exp = a.parts.exp;
84 return result;
85 }
86
87 if (is_float32_infinity(b)) {
88 if (is_float32_zero(a)) {
89 /* FIXME: zero * infinity */
90 result.bin = FLOAT32_NAN;
91 return result;
92 }
93
94 result.parts.fraction = b.parts.fraction;
95 result.parts.exp = b.parts.exp;
96 return result;
97 }
98
99 /* exp is signed so we can easy detect underflow */
100 exp = a.parts.exp + b.parts.exp;
101 exp -= FLOAT32_BIAS;
102
103 if (exp >= FLOAT32_MAX_EXPONENT) {
104 /* FIXME: overflow */
105 /* set infinity as result */
106 result.bin = FLOAT32_INF;
107 result.parts.sign = a.parts.sign ^ b.parts.sign;
108 return result;
109 }
110
111 if (exp < 0) {
112 /* FIXME: underflow */
113 /* return signed zero */
114 result.parts.fraction = 0x0;
115 result.parts.exp = 0x0;
116 return result;
117 }
118
119 frac1 = a.parts.fraction;
120 if (a.parts.exp > 0) {
121 frac1 |= FLOAT32_HIDDEN_BIT_MASK;
122 } else {
123 ++exp;
124 }
125
126 frac2 = b.parts.fraction;
127
128 if (b.parts.exp > 0) {
129 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
130 } else {
131 ++exp;
132 }
133
134 frac1 <<= 1; /* one bit space for rounding */
135
136 frac1 = frac1 * frac2;
137
138 /* round and return */
139 while ((exp < FLOAT32_MAX_EXPONENT) &&
140 (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 2)))) {
141 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left) */
142 ++exp;
143 frac1 >>= 1;
144 }
145
146 /* rounding */
147 /* ++frac1; FIXME: not works - without it is ok */
148 frac1 >>= 1; /* shift off rounding space */
149
150 if ((exp < FLOAT32_MAX_EXPONENT) &&
151 (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
152 ++exp;
153 frac1 >>= 1;
154 }
155
156 if (exp >= FLOAT32_MAX_EXPONENT) {
157 /* TODO: fix overflow */
158 /* return infinity*/
159 result.parts.exp = FLOAT32_MAX_EXPONENT;
160 result.parts.fraction = 0x0;
161 return result;
162 }
163
164 exp -= FLOAT32_FRACTION_SIZE;
165
166 if (exp <= FLOAT32_FRACTION_SIZE) {
167 /* denormalized number */
168 frac1 >>= 1; /* denormalize */
169
170 while ((frac1 > 0) && (exp < 0)) {
171 frac1 >>= 1;
172 ++exp;
173 }
174
175 if (frac1 == 0) {
176 /* FIXME : underflow */
177 result.parts.exp = 0;
178 result.parts.fraction = 0;
179 return result;
180 }
181 }
182
183 result.parts.exp = exp;
184 result.parts.fraction = frac1 & ((1 << FLOAT32_FRACTION_SIZE) - 1);
185
186 return result;
187}
188
189/** Multiply two double-precision floats.
190 *
191 * @param a First input operand.
192 * @param b Second input operand.
193 *
194 * @return Result of multiplication.
195 *
196 */
197float64 mul_float64(float64 a, float64 b)
198{
199 float64 result;
200 uint64_t frac1, frac2;
201 int32_t exp;
202
203 result.parts.sign = a.parts.sign ^ b.parts.sign;
204
205 if (is_float64_nan(a) || is_float64_nan(b)) {
206 /* TODO: fix SigNaNs */
207 if (is_float64_signan(a)) {
208 result.parts.fraction = a.parts.fraction;
209 result.parts.exp = a.parts.exp;
210 return result;
211 }
212 if (is_float64_signan(b)) { /* TODO: fix SigNaN */
213 result.parts.fraction = b.parts.fraction;
214 result.parts.exp = b.parts.exp;
215 return result;
216 }
217 /* set NaN as result */
218 result.bin = FLOAT64_NAN;
219 return result;
220 }
221
222 if (is_float64_infinity(a)) {
223 if (is_float64_zero(b)) {
224 /* FIXME: zero * infinity */
225 result.bin = FLOAT64_NAN;
226 return result;
227 }
228 result.parts.fraction = a.parts.fraction;
229 result.parts.exp = a.parts.exp;
230 return result;
231 }
232
233 if (is_float64_infinity(b)) {
234 if (is_float64_zero(a)) {
235 /* FIXME: zero * infinity */
236 result.bin = FLOAT64_NAN;
237 return result;
238 }
239 result.parts.fraction = b.parts.fraction;
240 result.parts.exp = b.parts.exp;
241 return result;
242 }
243
244 /* exp is signed so we can easy detect underflow */
245 exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
246
247 frac1 = a.parts.fraction;
248
249 if (a.parts.exp > 0) {
250 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
251 } else {
252 ++exp;
253 }
254
255 frac2 = b.parts.fraction;
256
257 if (b.parts.exp > 0) {
258 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
259 } else {
260 ++exp;
261 }
262
263 frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
264 frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
265
266 mul64(frac1, frac2, &frac1, &frac2);
267
268 frac1 |= (frac2 != 0);
269 if (frac1 & (0x1ll << 62)) {
270 frac1 <<= 1;
271 exp--;
272 }
273
274 result = finish_float64(exp, frac1, result.parts.sign);
275 return result;
276}
277
278/** Multiply two quadruple-precision floats.
279 *
280 * @param a First input operand.
281 * @param b Second input operand.
282 *
283 * @return Result of multiplication.
284 *
285 */
286float128 mul_float128(float128 a, float128 b)
287{
288 float128 result;
289 uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
290 int32_t exp;
291
292 result.parts.sign = a.parts.sign ^ b.parts.sign;
293
294 if (is_float128_nan(a) || is_float128_nan(b)) {
295 /* TODO: fix SigNaNs */
296 if (is_float128_signan(a)) {
297 result.parts.frac_hi = a.parts.frac_hi;
298 result.parts.frac_lo = a.parts.frac_lo;
299 result.parts.exp = a.parts.exp;
300 return result;
301 }
302 if (is_float128_signan(b)) { /* TODO: fix SigNaN */
303 result.parts.frac_hi = b.parts.frac_hi;
304 result.parts.frac_lo = b.parts.frac_lo;
305 result.parts.exp = b.parts.exp;
306 return result;
307 }
308 /* set NaN as result */
309 result.bin.hi = FLOAT128_NAN_HI;
310 result.bin.lo = FLOAT128_NAN_LO;
311 return result;
312 }
313
314 if (is_float128_infinity(a)) {
315 if (is_float128_zero(b)) {
316 /* FIXME: zero * infinity */
317 result.bin.hi = FLOAT128_NAN_HI;
318 result.bin.lo = FLOAT128_NAN_LO;
319 return result;
320 }
321 result.parts.frac_hi = a.parts.frac_hi;
322 result.parts.frac_lo = a.parts.frac_lo;
323 result.parts.exp = a.parts.exp;
324 return result;
325 }
326
327 if (is_float128_infinity(b)) {
328 if (is_float128_zero(a)) {
329 /* FIXME: zero * infinity */
330 result.bin.hi = FLOAT128_NAN_HI;
331 result.bin.lo = FLOAT128_NAN_LO;
332 return result;
333 }
334 result.parts.frac_hi = b.parts.frac_hi;
335 result.parts.frac_lo = b.parts.frac_lo;
336 result.parts.exp = b.parts.exp;
337 return result;
338 }
339
340 /* exp is signed so we can easy detect underflow */
341 exp = a.parts.exp + b.parts.exp - FLOAT128_BIAS - 1;
342
343 frac1_hi = a.parts.frac_hi;
344 frac1_lo = a.parts.frac_lo;
345
346 if (a.parts.exp > 0) {
347 or128(frac1_hi, frac1_lo,
348 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
349 &frac1_hi, &frac1_lo);
350 } else {
351 ++exp;
352 }
353
354 frac2_hi = b.parts.frac_hi;
355 frac2_lo = b.parts.frac_lo;
356
357 if (b.parts.exp > 0) {
358 or128(frac2_hi, frac2_lo,
359 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
360 &frac2_hi, &frac2_lo);
361 } else {
362 ++exp;
363 }
364
365 lshift128(frac2_hi, frac2_lo,
366 128 - FLOAT128_FRACTION_SIZE, &frac2_hi, &frac2_lo);
367
368 tmp_hi = frac1_hi;
369 tmp_lo = frac1_lo;
370 mul128(frac1_hi, frac1_lo, frac2_hi, frac2_lo,
371 &frac1_hi, &frac1_lo, &frac2_hi, &frac2_lo);
372 add128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo);
373 frac2_hi |= (frac2_lo != 0x0ll);
374
375 if ((FLOAT128_HIDDEN_BIT_MASK_HI << 1) <= frac1_hi) {
376 frac2_hi >>= 1;
377 if (frac1_lo & 0x1ll) {
378 frac2_hi |= (0x1ull < 64);
379 }
380 rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
381 ++exp;
382 }
383
384 result = finish_float128(exp, frac1_hi, frac1_lo, result.parts.sign, frac2_hi);
385 return result;
386}
387
388#ifdef float32_t
389
390float32_t __mulsf3(float32_t a, float32_t b)
391{
392 float32_u ua;
393 ua.val = a;
394
395 float32_u ub;
396 ub.val = b;
397
398 float32_u res;
399 res.data = mul_float32(ua.data, ub.data);
400
401 return res.val;
402}
403
404float32_t __aeabi_fmul(float32_t a, float32_t b)
405{
406 float32_u ua;
407 ua.val = a;
408
409 float32_u ub;
410 ub.val = b;
411
412 float32_u res;
413 res.data = mul_float32(ua.data, ub.data);
414
415 return res.val;
416}
417
418#endif
419
420#ifdef float64_t
421
422float64_t __muldf3(float64_t a, float64_t b)
423{
424 float64_u ua;
425 ua.val = a;
426
427 float64_u ub;
428 ub.val = b;
429
430 float64_u res;
431 res.data = mul_float64(ua.data, ub.data);
432
433 return res.val;
434}
435
436float64_t __aeabi_dmul(float64_t a, float64_t b)
437{
438 float64_u ua;
439 ua.val = a;
440
441 float64_u ub;
442 ub.val = b;
443
444 float64_u res;
445 res.data = mul_float64(ua.data, ub.data);
446
447 return res.val;
448}
449
450#endif
451
452#ifdef float128_t
453
454float128_t __multf3(float128_t a, float128_t b)
455{
456 float128_u ua;
457 ua.val = a;
458
459 float128_u ub;
460 ub.val = b;
461
462 float128_u res;
463 res.data = mul_float128(ua.data, ub.data);
464
465 return res.val;
466}
467
468void _Qp_mul(float128_t *c, float128_t *a, float128_t *b)
469{
470 *c = __multf3(*a, *b);
471}
472
473#endif
474
475/** @}
476 */
Note: See TracBrowser for help on using the repository browser.