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

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

32 bit float division added.
Some small bugs fixed.
Code cleanup.
  • Property mode set to 100644
File size: 8.3 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
33/** Multiply two 32 bit float numbers
34 *
35 */
36float32 mulFloat32(float32 a, float32 b)
37{
38 float32 result;
39 __u64 frac1, frac2;
40 __s32 exp;
41
42 result.parts.sign = a.parts.sign ^ b.parts.sign;
43
44 if (isFloat32NaN(a) || isFloat32NaN(b) ) {
45 /* TODO: fix SigNaNs */
46 if (isFloat32SigNaN(a)) {
47 result.parts.fraction = a.parts.fraction;
48 result.parts.exp = a.parts.exp;
49 return result;
50 };
51 if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
52 result.parts.fraction = b.parts.fraction;
53 result.parts.exp = b.parts.exp;
54 return result;
55 };
56 /* set NaN as result */
57 result.binary = FLOAT32_NAN;
58 return result;
59 };
60
61 if (isFloat32Infinity(a)) {
62 if (isFloat32Zero(b)) {
63 /* FIXME: zero * infinity */
64 result.binary = FLOAT32_NAN;
65 return result;
66 }
67 result.parts.fraction = a.parts.fraction;
68 result.parts.exp = a.parts.exp;
69 return result;
70 }
71
72 if (isFloat32Infinity(b)) {
73 if (isFloat32Zero(a)) {
74 /* FIXME: zero * infinity */
75 result.binary = FLOAT32_NAN;
76 return result;
77 }
78 result.parts.fraction = b.parts.fraction;
79 result.parts.exp = b.parts.exp;
80 return result;
81 }
82
83 /* exp is signed so we can easy detect underflow */
84 exp = a.parts.exp + b.parts.exp;
85 exp -= FLOAT32_BIAS;
86
87 if (exp >= FLOAT32_MAX_EXPONENT) {
88 /* FIXME: overflow */
89 /* set infinity as result */
90 result.binary = FLOAT32_INF;
91 result.parts.sign = a.parts.sign ^ b.parts.sign;
92 return result;
93 };
94
95 if (exp < 0) {
96 /* FIXME: underflow */
97 /* return signed zero */
98 result.parts.fraction = 0x0;
99 result.parts.exp = 0x0;
100 return result;
101 };
102
103 frac1 = a.parts.fraction;
104 if (a.parts.exp > 0) {
105 frac1 |= FLOAT32_HIDDEN_BIT_MASK;
106 } else {
107 ++exp;
108 };
109
110 frac2 = b.parts.fraction;
111
112 if (b.parts.exp > 0) {
113 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
114 } else {
115 ++exp;
116 };
117
118 frac1 <<= 1; /* one bit space for rounding */
119
120 frac1 = frac1 * frac2;
121/* round and return */
122
123 while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
124 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
125 ++exp;
126 frac1 >>= 1;
127 };
128
129 /* rounding */
130 /* ++frac1; FIXME: not works - without it is ok */
131 frac1 >>= 1; /* shift off rounding space */
132
133 if ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 1)))) {
134 ++exp;
135 frac1 >>= 1;
136 };
137
138 if (exp >= FLOAT32_MAX_EXPONENT ) {
139 /* TODO: fix overflow */
140 /* return infinity*/
141 result.parts.exp = FLOAT32_MAX_EXPONENT;
142 result.parts.fraction = 0x0;
143 return result;
144 }
145
146 exp -= FLOAT32_FRACTION_SIZE;
147
148 if (exp <= FLOAT32_FRACTION_SIZE) {
149 /* denormalized number */
150 frac1 >>= 1; /* denormalize */
151 while ((frac1 > 0) && (exp < 0)) {
152 frac1 >>= 1;
153 ++exp;
154 };
155 if (frac1 == 0) {
156 /* FIXME : underflow */
157 result.parts.exp = 0;
158 result.parts.fraction = 0;
159 return result;
160 };
161 };
162 result.parts.exp = exp;
163 result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
164
165 return result;
166
167}
168
169/** Multiply two 64 bit float numbers
170 *
171 */
172float64 mulFloat64(float64 a, float64 b)
173{
174 float64 result;
175 __u64 frac1, frac2;
176 __s32 exp;
177
178 result.parts.sign = a.parts.sign ^ b.parts.sign;
179
180 if (isFloat64NaN(a) || isFloat64NaN(b) ) {
181 /* TODO: fix SigNaNs */
182 if (isFloat64SigNaN(a)) {
183 result.parts.fraction = a.parts.fraction;
184 result.parts.exp = a.parts.exp;
185 return result;
186 };
187 if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
188 result.parts.fraction = b.parts.fraction;
189 result.parts.exp = b.parts.exp;
190 return result;
191 };
192 /* set NaN as result */
193 result.binary = FLOAT64_NAN;
194 return result;
195 };
196
197 if (isFloat64Infinity(a)) {
198 if (isFloat64Zero(b)) {
199 /* FIXME: zero * infinity */
200 result.binary = FLOAT64_NAN;
201 return result;
202 }
203 result.parts.fraction = a.parts.fraction;
204 result.parts.exp = a.parts.exp;
205 return result;
206 }
207
208 if (isFloat64Infinity(b)) {
209 if (isFloat64Zero(a)) {
210 /* FIXME: zero * infinity */
211 result.binary = FLOAT64_NAN;
212 return result;
213 }
214 result.parts.fraction = b.parts.fraction;
215 result.parts.exp = b.parts.exp;
216 return result;
217 }
218
219 /* exp is signed so we can easy detect underflow */
220 exp = a.parts.exp + b.parts.exp;
221 exp -= FLOAT64_BIAS;
222
223 if (exp >= FLOAT64_MAX_EXPONENT) {
224 /* FIXME: overflow */
225 /* set infinity as result */
226 result.binary = FLOAT64_INF;
227 result.parts.sign = a.parts.sign ^ b.parts.sign;
228 return result;
229 };
230
231 if (exp < 0) {
232 /* FIXME: underflow */
233 /* return signed zero */
234 result.parts.fraction = 0x0;
235 result.parts.exp = 0x0;
236 return result;
237 };
238
239 frac1 = a.parts.fraction;
240 if (a.parts.exp > 0) {
241 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
242 } else {
243 ++exp;
244 };
245
246 frac2 = b.parts.fraction;
247
248 if (b.parts.exp > 0) {
249 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
250 } else {
251 ++exp;
252 };
253
254 frac1 <<= 1; /* one bit space for rounding */
255
256 mul64integers(frac1, frac2, &frac1, &frac2);
257
258/* round and return */
259 /* FIXME: ugly soulution is to shift whole frac2 >> as in 32bit version
260 * Here is is more slower because we have to shift two numbers with carry
261 * Better is find first nonzero bit and make only one shift
262 * Third version is to shift both numbers a bit to right and result will be then
263 * placed in higher part of result. Then lower part will be good only for rounding.
264 */
265
266 while ((exp < FLOAT64_MAX_EXPONENT) && (frac2 > 0 )) {
267 frac1 >>= 1;
268 frac1 &= ((frac2 & 0x1) << 63);
269 frac2 >>= 1;
270 ++exp;
271 }
272
273 while ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ( (__u64)1 << (FLOAT64_FRACTION_SIZE + 2)))) {
274 ++exp;
275 frac1 >>= 1;
276 };
277
278 /* rounding */
279 /* ++frac1; FIXME: not works - without it is ok */
280 frac1 >>= 1; /* shift off rounding space */
281
282 if ((exp < FLOAT64_MAX_EXPONENT) && (frac1 >= ((__u64)1 << (FLOAT64_FRACTION_SIZE + 1)))) {
283 ++exp;
284 frac1 >>= 1;
285 };
286
287 if (exp >= FLOAT64_MAX_EXPONENT ) {
288 /* TODO: fix overflow */
289 /* return infinity*/
290 result.parts.exp = FLOAT64_MAX_EXPONENT;
291 result.parts.fraction = 0x0;
292 return result;
293 }
294
295 exp -= FLOAT64_FRACTION_SIZE;
296
297 if (exp <= FLOAT64_FRACTION_SIZE) {
298 /* denormalized number */
299 frac1 >>= 1; /* denormalize */
300 while ((frac1 > 0) && (exp < 0)) {
301 frac1 >>= 1;
302 ++exp;
303 };
304 if (frac1 == 0) {
305 /* FIXME : underflow */
306 result.parts.exp = 0;
307 result.parts.fraction = 0;
308 return result;
309 };
310 };
311 result.parts.exp = exp;
312 result.parts.fraction = frac1 & ( ((__u64)1 << FLOAT64_FRACTION_SIZE) - 1);
313
314 return result;
315
316}
317
318/** Multiply two 64 bit numbers and return result in two parts
319 * @param a first operand
320 * @param b second operand
321 * @param lo lower part from result
322 * @param hi higher part of result
323 */
324void mul64integers(__u64 a,__u64 b, __u64 *lo, __u64 *hi)
325{
326 __u64 low, high, middle1, middle2;
327 __u32 alow, blow;
328
329 alow = a & 0xFFFFFFFF;
330 blow = b & 0xFFFFFFFF;
331
332 a <<= 32;
333 b <<= 32;
334
335 low = (__u64)alow * blow;
336 middle1 = a * blow;
337 middle2 = alow * b;
338 high = a * b;
339
340 middle1 += middle2;
341 high += ((__u64)(middle1 < middle2) << 32) + (middle1 >> 32);
342 middle1 <<= 32;
343 low += middle1;
344 high += (low < middle1);
345 *lo = low;
346 *hi = high;
347 return;
348}
349
350
Note: See TracBrowser for help on using the repository browser.