source: mainline/softfloat/generic/mul.c@ ce5bcb4

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since ce5bcb4 was aa59fa0, checked in by Josef Cejka <malyzelenyhnus@…>, 19 years ago

SoftFloat integrated into HelenOS uspace.
  • Property mode set to 100644
File size: 6.7 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#include<sftypes.h>
30#include<mul.h>
31#include<comparison.h>
32#include<common.h>
33
34/** Multiply two 32 bit float numbers
35 *
36 */
37float32 mulFloat32(float32 a, float32 b)
38{
39 float32 result;
40 uint64_t frac1, frac2;
41 int32_t exp;
42
43 result.parts.sign = a.parts.sign ^ b.parts.sign;
44
45 if (isFloat32NaN(a) || isFloat32NaN(b) ) {
46 /* TODO: fix SigNaNs */
47 if (isFloat32SigNaN(a)) {
48 result.parts.fraction = a.parts.fraction;
49 result.parts.exp = a.parts.exp;
50 return result;
51 };
52 if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
53 result.parts.fraction = b.parts.fraction;
54 result.parts.exp = b.parts.exp;
55 return result;
56 };
57 /* set NaN as result */
58 result.binary = FLOAT32_NAN;
59 return result;
60 };
61
62 if (isFloat32Infinity(a)) {
63 if (isFloat32Zero(b)) {
64 /* FIXME: zero * infinity */
65 result.binary = FLOAT32_NAN;
66 return result;
67 }
68 result.parts.fraction = a.parts.fraction;
69 result.parts.exp = a.parts.exp;
70 return result;
71 }
72
73 if (isFloat32Infinity(b)) {
74 if (isFloat32Zero(a)) {
75 /* FIXME: zero * infinity */
76 result.binary = FLOAT32_NAN;
77 return result;
78 }
79 result.parts.fraction = b.parts.fraction;
80 result.parts.exp = b.parts.exp;
81 return result;
82 }
83
84 /* exp is signed so we can easy detect underflow */
85 exp = a.parts.exp + b.parts.exp;
86 exp -= FLOAT32_BIAS;
87
88 if (exp >= FLOAT32_MAX_EXPONENT) {
89 /* FIXME: overflow */
90 /* set infinity as result */
91 result.binary = FLOAT32_INF;
92 result.parts.sign = a.parts.sign ^ b.parts.sign;
93 return result;
94 };
95
96 if (exp < 0) {
97 /* FIXME: underflow */
98 /* return signed zero */
99 result.parts.fraction = 0x0;
100 result.parts.exp = 0x0;
101 return result;
102 };
103
104 frac1 = a.parts.fraction;
105 if (a.parts.exp > 0) {
106 frac1 |= FLOAT32_HIDDEN_BIT_MASK;
107 } else {
108 ++exp;
109 };
110
111 frac2 = b.parts.fraction;
112
113 if (b.parts.exp > 0) {
114 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
115 } else {
116 ++exp;
117 };
118
119 frac1 <<= 1; /* one bit space for rounding */
120
121 frac1 = frac1 * frac2;
122/* round and return */
123
124 while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
125 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
126 ++exp;
127 frac1 >>= 1;
128 };
129
130 /* rounding */
131 /* ++frac1; FIXME: not works - without it is ok */
132 frac1 >>= 1; /* shift off rounding space */
133
134 if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
135 ++exp;
136 frac1 >>= 1;
137 };
138
139 if (exp >= FLOAT32_MAX_EXPONENT ) {
140 /* TODO: fix overflow */
141 /* return infinity*/
142 result.parts.exp = FLOAT32_MAX_EXPONENT;
143 result.parts.fraction = 0x0;
144 return result;
145 }
146
147 exp -= FLOAT32_FRACTION_SIZE;
148
149 if (exp <= FLOAT32_FRACTION_SIZE) {
150 /* denormalized number */
151 frac1 >>= 1; /* denormalize */
152 while ((frac1 > 0) && (exp < 0)) {
153 frac1 >>= 1;
154 ++exp;
155 };
156 if (frac1 == 0) {
157 /* FIXME : underflow */
158 result.parts.exp = 0;
159 result.parts.fraction = 0;
160 return result;
161 };
162 };
163 result.parts.exp = exp;
164 result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
165
166 return result;
167
168}
169
170/** Multiply two 64 bit float numbers
171 *
172 */
173float64 mulFloat64(float64 a, float64 b)
174{
175 float64 result;
176 uint64_t frac1, frac2;
177 int32_t exp;
178
179 result.parts.sign = a.parts.sign ^ b.parts.sign;
180
181 if (isFloat64NaN(a) || isFloat64NaN(b) ) {
182 /* TODO: fix SigNaNs */
183 if (isFloat64SigNaN(a)) {
184 result.parts.fraction = a.parts.fraction;
185 result.parts.exp = a.parts.exp;
186 return result;
187 };
188 if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
189 result.parts.fraction = b.parts.fraction;
190 result.parts.exp = b.parts.exp;
191 return result;
192 };
193 /* set NaN as result */
194 result.binary = FLOAT64_NAN;
195 return result;
196 };
197
198 if (isFloat64Infinity(a)) {
199 if (isFloat64Zero(b)) {
200 /* FIXME: zero * infinity */
201 result.binary = FLOAT64_NAN;
202 return result;
203 }
204 result.parts.fraction = a.parts.fraction;
205 result.parts.exp = a.parts.exp;
206 return result;
207 }
208
209 if (isFloat64Infinity(b)) {
210 if (isFloat64Zero(a)) {
211 /* FIXME: zero * infinity */
212 result.binary = FLOAT64_NAN;
213 return result;
214 }
215 result.parts.fraction = b.parts.fraction;
216 result.parts.exp = b.parts.exp;
217 return result;
218 }
219
220 /* exp is signed so we can easy detect underflow */
221 exp = a.parts.exp + b.parts.exp - FLOAT64_BIAS;
222
223 frac1 = a.parts.fraction;
224
225 if (a.parts.exp > 0) {
226 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
227 } else {
228 ++exp;
229 };
230
231 frac2 = b.parts.fraction;
232
233 if (b.parts.exp > 0) {
234 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
235 } else {
236 ++exp;
237 };
238
239 frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
240 frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
241
242 mul64integers(frac1, frac2, &frac1, &frac2);
243
244 frac2 |= (frac1 != 0);
245 if (frac2 & (0x1ll << 62)) {
246 frac2 <<= 1;
247 exp--;
248 }
249
250 result = finishFloat64(exp, frac2, result.parts.sign);
251 return result;
252}
253
254/** Multiply two 64 bit numbers and return result in two parts
255 * @param a first operand
256 * @param b second operand
257 * @param lo lower part from result
258 * @param hi higher part of result
259 */
260void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi)
261{
262 uint64_t low, high, middle1, middle2;
263 uint32_t alow, blow;
264
265 alow = a & 0xFFFFFFFF;
266 blow = b & 0xFFFFFFFF;
267
268 a >>= 32;
269 b >>= 32;
270
271 low = ((uint64_t)alow) * blow;
272 middle1 = a * blow;
273 middle2 = alow * b;
274 high = a * b;
275
276 middle1 += middle2;
277 high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
278 middle1 <<= 32;
279 low += middle1;
280 high += (low < middle1);
281 *lo = low;
282 *hi = high;
283
284 return;
285}
286
287
Note: See TracBrowser for help on using the repository browser.