source: mainline/uspace/lib/softfloat/mul.c@ cf13b17

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

software floating point overhaul
use proper type mapping
fix cosine calculation
  • 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.