source: mainline/uspace/lib/softfloat/generic/mul.c@ d61c3a7

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

softfloat: slightly improve coding style (no change in functionality)
  • Property mode set to 100644
File size: 6.8 KB
Line 
1/*
2 * Copyright (c) 2005 Josef Cejka
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 *
9 * - Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * - Redistributions in binary form must reproduce the above copyright
12 * notice, this list of conditions and the following disclaimer in the
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
15 * derived from this software without specific prior written permission.
16 *
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
19 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
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
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
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29/** @addtogroup softfloat
30 * @{
31 */
32/** @file
33 */
34
35#include <sftypes.h>
36#include <mul.h>
37#include <comparison.h>
38#include <common.h>
39
40/** Multiply two 32 bit float numbers
41 *
42 */
43float32 mulFloat32(float32 a, float32 b)
44{
45 float32 result;
46 uint64_t frac1, frac2;
47 int32_t exp;
48
49 result.parts.sign = a.parts.sign ^ b.parts.sign;
50
51 if (isFloat32NaN(a) || isFloat32NaN(b) ) {
52 /* TODO: fix SigNaNs */
53 if (isFloat32SigNaN(a)) {
54 result.parts.fraction = a.parts.fraction;
55 result.parts.exp = a.parts.exp;
56 return result;
57 };
58 if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
59 result.parts.fraction = b.parts.fraction;
60 result.parts.exp = b.parts.exp;
61 return result;
62 };
63 /* set NaN as result */
64 result.binary = FLOAT32_NAN;
65 return result;
66 };
67
68 if (isFloat32Infinity(a)) {
69 if (isFloat32Zero(b)) {
70 /* FIXME: zero * infinity */
71 result.binary = FLOAT32_NAN;
72 return result;
73 }
74 result.parts.fraction = a.parts.fraction;
75 result.parts.exp = a.parts.exp;
76 return result;
77 }
78
79 if (isFloat32Infinity(b)) {
80 if (isFloat32Zero(a)) {
81 /* FIXME: zero * infinity */
82 result.binary = FLOAT32_NAN;
83 return result;
84 }
85 result.parts.fraction = b.parts.fraction;
86 result.parts.exp = b.parts.exp;
87 return result;
88 }
89
90 /* exp is signed so we can easy detect underflow */
91 exp = a.parts.exp + b.parts.exp;
92 exp -= FLOAT32_BIAS;
93
94 if (exp >= FLOAT32_MAX_EXPONENT) {
95 /* FIXME: overflow */
96 /* set infinity as result */
97 result.binary = FLOAT32_INF;
98 result.parts.sign = a.parts.sign ^ b.parts.sign;
99 return result;
100 };
101
102 if (exp < 0) {
103 /* FIXME: underflow */
104 /* return signed zero */
105 result.parts.fraction = 0x0;
106 result.parts.exp = 0x0;
107 return result;
108 };
109
110 frac1 = a.parts.fraction;
111 if (a.parts.exp > 0) {
112 frac1 |= FLOAT32_HIDDEN_BIT_MASK;
113 } else {
114 ++exp;
115 };
116
117 frac2 = b.parts.fraction;
118
119 if (b.parts.exp > 0) {
120 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
121 } else {
122 ++exp;
123 };
124
125 frac1 <<= 1; /* one bit space for rounding */
126
127 frac1 = frac1 * frac2;
128/* round and return */
129
130 while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
131 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
132 ++exp;
133 frac1 >>= 1;
134 };
135
136 /* rounding */
137 /* ++frac1; FIXME: not works - without it is ok */
138 frac1 >>= 1; /* shift off rounding space */
139
140 if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
141 ++exp;
142 frac1 >>= 1;
143 };
144
145 if (exp >= FLOAT32_MAX_EXPONENT ) {
146 /* TODO: fix overflow */
147 /* return infinity*/
148 result.parts.exp = FLOAT32_MAX_EXPONENT;
149 result.parts.fraction = 0x0;
150 return result;
151 }
152
153 exp -= FLOAT32_FRACTION_SIZE;
154
155 if (exp <= FLOAT32_FRACTION_SIZE) {
156 /* denormalized number */
157 frac1 >>= 1; /* denormalize */
158 while ((frac1 > 0) && (exp < 0)) {
159 frac1 >>= 1;
160 ++exp;
161 };
162 if (frac1 == 0) {
163 /* FIXME : underflow */
164 result.parts.exp = 0;
165 result.parts.fraction = 0;
166 return result;
167 };
168 };
169 result.parts.exp = exp;
170 result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
171
172 return result;
173
174}
175
176/** Multiply two 64 bit float numbers
177 *
178 */
179float64 mulFloat64(float64 a, float64 b)
180{
181 float64 result;
182 uint64_t frac1, frac2;
183 int32_t exp;
184
185 result.parts.sign = a.parts.sign ^ b.parts.sign;
186
187 if (isFloat64NaN(a) || isFloat64NaN(b) ) {
188 /* TODO: fix SigNaNs */
189 if (isFloat64SigNaN(a)) {
190 result.parts.fraction = a.parts.fraction;
191 result.parts.exp = a.parts.exp;
192 return result;
193 };
194 if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
195 result.parts.fraction = b.parts.fraction;
196 result.parts.exp = b.parts.exp;
197 return result;
198 };
199 /* set NaN as result */
200 result.binary = FLOAT64_NAN;
201 return result;
202 };
203
204 if (isFloat64Infinity(a)) {
205 if (isFloat64Zero(b)) {
206 /* FIXME: zero * infinity */
207 result.binary = FLOAT64_NAN;
208 return result;
209 }
210 result.parts.fraction = a.parts.fraction;
211 result.parts.exp = a.parts.exp;
212 return result;
213 }
214
215 if (isFloat64Infinity(b)) {
216 if (isFloat64Zero(a)) {
217 /* FIXME: zero * infinity */
218 result.binary = FLOAT64_NAN;
219 return result;
220 }
221 result.parts.fraction = b.parts.fraction;
222 result.parts.exp = b.parts.exp;
223 return result;
224 }
225
226 /* exp is signed so we can easy detect underflow */
227 exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
228
229 frac1 = a.parts.fraction;
230
231 if (a.parts.exp > 0) {
232 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
233 } else {
234 ++exp;
235 };
236
237 frac2 = b.parts.fraction;
238
239 if (b.parts.exp > 0) {
240 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
241 } else {
242 ++exp;
243 };
244
245 frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
246 frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
247
248 mul64integers(frac1, frac2, &frac1, &frac2);
249
250 frac2 |= (frac1 != 0);
251 if (frac2 & (0x1ll << 62)) {
252 frac2 <<= 1;
253 exp--;
254 }
255
256 result = finishFloat64(exp, frac2, result.parts.sign);
257 return result;
258}
259
260/** Multiply two 64 bit numbers and return result in two parts
261 * @param a first operand
262 * @param b second operand
263 * @param lo lower part from result
264 * @param hi higher part of result
265 */
266void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi)
267{
268 uint64_t low, high, middle1, middle2;
269 uint32_t alow, blow;
270
271 alow = a & 0xFFFFFFFF;
272 blow = b & 0xFFFFFFFF;
273
274 a >>= 32;
275 b >>= 32;
276
277 low = ((uint64_t)alow) * blow;
278 middle1 = a * blow;
279 middle2 = alow * b;
280 high = a * b;
281
282 middle1 += middle2;
283 high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
284 middle1 <<= 32;
285 low += middle1;
286 high += (low < middle1);
287 *lo = low;
288 *hi = high;
289
290 return;
291}
292
293/** @}
294 */
Note: See TracBrowser for help on using the repository browser.