source: mainline/uspace/lib/softfloat/mul.c@ 22fb7ab

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since 22fb7ab 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
RevLine 
[b5440cf]1/*
[df4ed85]2 * Copyright (c) 2005 Josef Cejka
[c67aff2]3 * Copyright (c) 2011 Petr Koupy
[b5440cf]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
[750636a]30/** @addtogroup softfloat
[846848a6]31 * @{
32 */
[c67aff2]33/** @file Multiplication functions.
[846848a6]34 */
35
[2416085]36#include "mul.h"
37#include "comparison.h"
38#include "common.h"
[b5440cf]39
[88d5c1e]40/** Multiply two single-precision floats.
[3af72dc]41 *
[c67aff2]42 * @param a First input operand.
43 * @param b Second input operand.
[88d5c1e]44 *
[c67aff2]45 * @return Result of multiplication.
[88d5c1e]46 *
[3af72dc]47 */
[88d5c1e]48float32 mul_float32(float32 a, float32 b)
[3af72dc]49{
50 float32 result;
[aa59fa0]51 uint64_t frac1, frac2;
52 int32_t exp;
[88d5c1e]53
[3af72dc]54 result.parts.sign = a.parts.sign ^ b.parts.sign;
55
[88d5c1e]56 if (is_float32_nan(a) || is_float32_nan(b)) {
[3af72dc]57 /* TODO: fix SigNaNs */
[88d5c1e]58 if (is_float32_signan(a)) {
[1266543]59 result.parts.fraction = a.parts.fraction;
[3af72dc]60 result.parts.exp = a.parts.exp;
61 return result;
[c67aff2]62 }
[c0c38c7c]63
[88d5c1e]64 if (is_float32_signan(b)) { /* TODO: fix SigNaN */
[1266543]65 result.parts.fraction = b.parts.fraction;
[3af72dc]66 result.parts.exp = b.parts.exp;
67 return result;
[c67aff2]68 }
[c0c38c7c]69
[3af72dc]70 /* set NaN as result */
[88d5c1e]71 result.bin = FLOAT32_NAN;
[3af72dc]72 return result;
[c67aff2]73 }
[88d5c1e]74
75 if (is_float32_infinity(a)) {
76 if (is_float32_zero(b)) {
[3af72dc]77 /* FIXME: zero * infinity */
[88d5c1e]78 result.bin = FLOAT32_NAN;
[3af72dc]79 return result;
80 }
[c0c38c7c]81
[1266543]82 result.parts.fraction = a.parts.fraction;
[3af72dc]83 result.parts.exp = a.parts.exp;
84 return result;
85 }
[88d5c1e]86
87 if (is_float32_infinity(b)) {
88 if (is_float32_zero(a)) {
[3af72dc]89 /* FIXME: zero * infinity */
[88d5c1e]90 result.bin = FLOAT32_NAN;
[3af72dc]91 return result;
92 }
[c0c38c7c]93
[1266543]94 result.parts.fraction = b.parts.fraction;
[3af72dc]95 result.parts.exp = b.parts.exp;
96 return result;
97 }
[88d5c1e]98
[3af72dc]99 /* exp is signed so we can easy detect underflow */
100 exp = a.parts.exp + b.parts.exp;
101 exp -= FLOAT32_BIAS;
102
[bff16dd]103 if (exp >= FLOAT32_MAX_EXPONENT) {
[3af72dc]104 /* FIXME: overflow */
105 /* set infinity as result */
[88d5c1e]106 result.bin = FLOAT32_INF;
[bff16dd]107 result.parts.sign = a.parts.sign ^ b.parts.sign;
[3af72dc]108 return result;
[c67aff2]109 }
[3af72dc]110
[c0c38c7c]111 if (exp < 0) {
[3af72dc]112 /* FIXME: underflow */
113 /* return signed zero */
[1266543]114 result.parts.fraction = 0x0;
[3af72dc]115 result.parts.exp = 0x0;
116 return result;
[c67aff2]117 }
[3af72dc]118
[1266543]119 frac1 = a.parts.fraction;
[bff16dd]120 if (a.parts.exp > 0) {
[1266543]121 frac1 |= FLOAT32_HIDDEN_BIT_MASK;
[3af72dc]122 } else {
123 ++exp;
[c67aff2]124 }
[3af72dc]125
[1266543]126 frac2 = b.parts.fraction;
[88d5c1e]127
[bff16dd]128 if (b.parts.exp > 0) {
[1266543]129 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
[3af72dc]130 } else {
131 ++exp;
[c67aff2]132 }
[88d5c1e]133
[1266543]134 frac1 <<= 1; /* one bit space for rounding */
[88d5c1e]135
[1266543]136 frac1 = frac1 * frac2;
[88d5c1e]137
[c67aff2]138 /* round and return */
[88d5c1e]139 while ((exp < FLOAT32_MAX_EXPONENT) &&
140 (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 2)))) {
[c67aff2]141 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left) */
[3af72dc]142 ++exp;
[1266543]143 frac1 >>= 1;
[c67aff2]144 }
[88d5c1e]145
[3af72dc]146 /* rounding */
[1266543]147 /* ++frac1; FIXME: not works - without it is ok */
148 frac1 >>= 1; /* shift off rounding space */
[3af72dc]149
[88d5c1e]150 if ((exp < FLOAT32_MAX_EXPONENT) &&
151 (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
[3af72dc]152 ++exp;
[1266543]153 frac1 >>= 1;
[c67aff2]154 }
[88d5c1e]155
156 if (exp >= FLOAT32_MAX_EXPONENT) {
[3af72dc]157 /* TODO: fix overflow */
158 /* return infinity*/
[bff16dd]159 result.parts.exp = FLOAT32_MAX_EXPONENT;
[1266543]160 result.parts.fraction = 0x0;
[3af72dc]161 return result;
162 }
163
[1266543]164 exp -= FLOAT32_FRACTION_SIZE;
[88d5c1e]165
166 if (exp <= FLOAT32_FRACTION_SIZE) {
[3af72dc]167 /* denormalized number */
[1266543]168 frac1 >>= 1; /* denormalize */
[c0c38c7c]169
[1266543]170 while ((frac1 > 0) && (exp < 0)) {
171 frac1 >>= 1;
[3af72dc]172 ++exp;
[c67aff2]173 }
[c0c38c7c]174
[1266543]175 if (frac1 == 0) {
[3af72dc]176 /* FIXME : underflow */
[c67aff2]177 result.parts.exp = 0;
178 result.parts.fraction = 0;
179 return result;
180 }
181 }
[c0c38c7c]182
[3af72dc]183 result.parts.exp = exp;
[c67aff2]184 result.parts.fraction = frac1 & ((1 << FLOAT32_FRACTION_SIZE) - 1);
[bff16dd]185
[88d5c1e]186 return result;
[bff16dd]187}
188
[88d5c1e]189/** Multiply two double-precision floats.
[bff16dd]190 *
[c67aff2]191 * @param a First input operand.
192 * @param b Second input operand.
[88d5c1e]193 *
[c67aff2]194 * @return Result of multiplication.
[88d5c1e]195 *
[bff16dd]196 */
[88d5c1e]197float64 mul_float64(float64 a, float64 b)
[bff16dd]198{
199 float64 result;
[aa59fa0]200 uint64_t frac1, frac2;
201 int32_t exp;
[88d5c1e]202
[bff16dd]203 result.parts.sign = a.parts.sign ^ b.parts.sign;
204
[88d5c1e]205 if (is_float64_nan(a) || is_float64_nan(b)) {
[bff16dd]206 /* TODO: fix SigNaNs */
[88d5c1e]207 if (is_float64_signan(a)) {
[1266543]208 result.parts.fraction = a.parts.fraction;
[bff16dd]209 result.parts.exp = a.parts.exp;
210 return result;
[c67aff2]211 }
[88d5c1e]212 if (is_float64_signan(b)) { /* TODO: fix SigNaN */
[1266543]213 result.parts.fraction = b.parts.fraction;
[bff16dd]214 result.parts.exp = b.parts.exp;
215 return result;
[c67aff2]216 }
[bff16dd]217 /* set NaN as result */
[88d5c1e]218 result.bin = FLOAT64_NAN;
[bff16dd]219 return result;
[c67aff2]220 }
[88d5c1e]221
222 if (is_float64_infinity(a)) {
223 if (is_float64_zero(b)) {
[bff16dd]224 /* FIXME: zero * infinity */
[88d5c1e]225 result.bin = FLOAT64_NAN;
[bff16dd]226 return result;
227 }
[1266543]228 result.parts.fraction = a.parts.fraction;
[bff16dd]229 result.parts.exp = a.parts.exp;
230 return result;
231 }
[88d5c1e]232
233 if (is_float64_infinity(b)) {
234 if (is_float64_zero(a)) {
[bff16dd]235 /* FIXME: zero * infinity */
[88d5c1e]236 result.bin = FLOAT64_NAN;
[bff16dd]237 return result;
238 }
[1266543]239 result.parts.fraction = b.parts.fraction;
[bff16dd]240 result.parts.exp = b.parts.exp;
241 return result;
242 }
[88d5c1e]243
[bff16dd]244 /* exp is signed so we can easy detect underflow */
[e979fea]245 exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
[bff16dd]246
[1266543]247 frac1 = a.parts.fraction;
[88d5c1e]248
[bff16dd]249 if (a.parts.exp > 0) {
[1266543]250 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
[bff16dd]251 } else {
252 ++exp;
[c67aff2]253 }
[bff16dd]254
[1266543]255 frac2 = b.parts.fraction;
[88d5c1e]256
[bff16dd]257 if (b.parts.exp > 0) {
[1266543]258 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
[bff16dd]259 } else {
260 ++exp;
[c67aff2]261 }
[88d5c1e]262
[e979fea]263 frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
264 frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
[88d5c1e]265
[c67aff2]266 mul64(frac1, frac2, &frac1, &frac2);
[88d5c1e]267
[c67aff2]268 frac1 |= (frac2 != 0);
269 if (frac1 & (0x1ll << 62)) {
270 frac1 <<= 1;
[e979fea]271 exp--;
[bff16dd]272 }
[88d5c1e]273
274 result = finish_float64(exp, frac1, result.parts.sign);
[e979fea]275 return result;
[12c6f2d]276}
277
[88d5c1e]278/** Multiply two quadruple-precision floats.
[c67aff2]279 *
280 * @param a First input operand.
281 * @param b Second input operand.
[88d5c1e]282 *
[c67aff2]283 * @return Result of multiplication.
[88d5c1e]284 *
[bff16dd]285 */
[88d5c1e]286float128 mul_float128(float128 a, float128 b)
[bff16dd]287{
[c67aff2]288 float128 result;
289 uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
290 int32_t exp;
[88d5c1e]291
[c67aff2]292 result.parts.sign = a.parts.sign ^ b.parts.sign;
[88d5c1e]293
294 if (is_float128_nan(a) || is_float128_nan(b)) {
[c67aff2]295 /* TODO: fix SigNaNs */
[88d5c1e]296 if (is_float128_signan(a)) {
[c67aff2]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 }
[88d5c1e]302 if (is_float128_signan(b)) { /* TODO: fix SigNaN */
[c67aff2]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 */
[88d5c1e]309 result.bin.hi = FLOAT128_NAN_HI;
310 result.bin.lo = FLOAT128_NAN_LO;
[c67aff2]311 return result;
312 }
[88d5c1e]313
314 if (is_float128_infinity(a)) {
315 if (is_float128_zero(b)) {
[c67aff2]316 /* FIXME: zero * infinity */
[88d5c1e]317 result.bin.hi = FLOAT128_NAN_HI;
318 result.bin.lo = FLOAT128_NAN_LO;
[c67aff2]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 }
[88d5c1e]326
327 if (is_float128_infinity(b)) {
328 if (is_float128_zero(a)) {
[c67aff2]329 /* FIXME: zero * infinity */
[88d5c1e]330 result.bin.hi = FLOAT128_NAN_HI;
331 result.bin.lo = FLOAT128_NAN_LO;
[c67aff2]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 }
[88d5c1e]339
[c67aff2]340 /* exp is signed so we can easy detect underflow */
341 exp = a.parts.exp + b.parts.exp - FLOAT128_BIAS - 1;
[88d5c1e]342
[c67aff2]343 frac1_hi = a.parts.frac_hi;
344 frac1_lo = a.parts.frac_lo;
[88d5c1e]345
[c67aff2]346 if (a.parts.exp > 0) {
347 or128(frac1_hi, frac1_lo,
[88d5c1e]348 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
349 &frac1_hi, &frac1_lo);
[c67aff2]350 } else {
351 ++exp;
352 }
[88d5c1e]353
[c67aff2]354 frac2_hi = b.parts.frac_hi;
355 frac2_lo = b.parts.frac_lo;
[88d5c1e]356
[c67aff2]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 }
[88d5c1e]364
[c67aff2]365 lshift128(frac2_hi, frac2_lo,
366 128 - FLOAT128_FRACTION_SIZE, &frac2_hi, &frac2_lo);
[88d5c1e]367
[c67aff2]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);
[88d5c1e]374
[c67aff2]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 }
[88d5c1e]383
384 result = finish_float128(exp, frac1_hi, frac1_lo, result.parts.sign, frac2_hi);
[c67aff2]385 return result;
[bff16dd]386}
[3af72dc]387
[c0c38c7c]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
[231a60a]475/** @}
[846848a6]476 */
Note: See TracBrowser for help on using the repository browser.