source: mainline/softfloat/generic/arithmetic.c@ 3af72dc

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

Function for multiplying added to softfloat lib.
  • Property mode set to 100644
File size: 7.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<arithmetic.h>
31#include<comparison.h>
32
33/** Add two Float32 numbers with same signs
34 */
35float32 addFloat32(float32 a, float32 b)
36{
37 int expdiff;
38 __u32 exp1,exp2,mant1,mant2;
39
40 expdiff=a.parts.exp - b.parts.exp;
41 if (expdiff<0) {
42 if (isFloat32NaN(b)) {
43 //TODO: fix SigNaN
44 if (isFloat32SigNaN(b)) {
45 };
46
47 return b;
48 };
49
50 if (b.parts.exp==0xFF) {
51 return b;
52 }
53
54 mant1=b.parts.mantisa;
55 exp1=b.parts.exp;
56 mant2=a.parts.mantisa;
57 exp2=a.parts.exp;
58 expdiff*=-1;
59 } else {
60 if (isFloat32NaN(a)) {
61 //TODO: fix SigNaN
62 if ((isFloat32SigNaN(a))||(isFloat32SigNaN(b))) {
63 };
64 return a;
65 };
66
67 if (a.parts.exp==0xFF) {
68 return a;
69 }
70
71 mant1=a.parts.mantisa;
72 exp1=a.parts.exp;
73 mant2=b.parts.mantisa;
74 exp2=b.parts.exp;
75 };
76
77 if (exp1==0) {
78 //both are denormalized
79 mant1+=mant2;
80 if (mant1&0xF00000) {
81 a.parts.exp=1;
82 };
83 a.parts.mantisa=mant1;
84 return a;
85 };
86
87 // create some space for rounding
88 mant1<<=6;
89 mant2<<=6;
90
91 mant1|=0x20000000; //add hidden bit
92
93
94 if (exp2==0) {
95 --expdiff;
96 } else {
97 mant2|=0x20000000; //hidden bit
98 };
99
100 if (expdiff>24) {
101 goto done;
102 };
103
104 mant2>>=expdiff;
105 mant1+=mant2;
106done:
107 if (mant1&0x40000000) {
108 ++exp1;
109 mant1>>=1;
110 };
111
112 //rounding - if first bit after mantisa is set then round up
113 mant1+=0x20;
114
115 if (mant1&0x40000000) {
116 ++exp1;
117 mant1>>=1;
118 };
119
120 a.parts.exp=exp1;
121 a.parts.mantisa = ((mant1&(~0x20000000))>>6); /*Clear hidden bit and shift */
122 return a;
123};
124
125/** Subtract two float32 numbers with same signs
126 */
127float32 subFloat32(float32 a, float32 b)
128{
129 int expdiff;
130 __u32 exp1,exp2,mant1,mant2;
131 float32 result;
132
133 result.f = 0;
134
135 expdiff=a.parts.exp - b.parts.exp;
136 if ((expdiff<0)||((expdiff==0)&&(a.parts.mantisa<b.parts.mantisa))) {
137 if (isFloat32NaN(b)) {
138 //TODO: fix SigNaN
139 if (isFloat32SigNaN(b)) {
140 };
141 return b;
142 };
143
144 if (b.parts.exp==0xFF) {
145 b.parts.sign=!b.parts.sign; /* num -(+-inf) = -+inf */
146 return b;
147 }
148
149 result.parts.sign = !a.parts.sign;
150
151 mant1=b.parts.mantisa;
152 exp1=b.parts.exp;
153 mant2=a.parts.mantisa;
154 exp2=a.parts.exp;
155 expdiff*=-1;
156 } else {
157 if (isFloat32NaN(a)) {
158 //TODO: fix SigNaN
159 if ((isFloat32SigNaN(a))||(isFloat32SigNaN(b))) {
160 };
161 return a;
162 };
163
164 if (a.parts.exp==0xFF) {
165 if (b.parts.exp==0xFF) {
166 /* inf - inf => nan */
167 //TODO: fix exception
168 result.binary = FLOAT32_NAN;
169 return result;
170 };
171 return a;
172 }
173
174 result.parts.sign = a.parts.sign;
175
176 mant1=a.parts.mantisa;
177 exp1=a.parts.exp;
178 mant2=b.parts.mantisa;
179 exp2=b.parts.exp;
180
181
182
183 };
184
185 if (exp1==0) {
186 //both are denormalized
187 result.parts.mantisa=mant1-mant2;
188 if (result.parts.mantisa>mant1) {
189 //TODO: underflow exception
190 return result;
191 };
192 result.parts.exp=0;
193 return result;
194 };
195
196 // create some space for rounding
197 mant1<<=6;
198 mant2<<=6;
199
200 mant1|=0x20000000; //add hidden bit
201
202
203 if (exp2==0) {
204 --expdiff;
205 } else {
206 mant2|=0x20000000; //hidden bit
207 };
208
209 if (expdiff>24) {
210 goto done;
211 };
212
213 mant1 = mant1-(mant2>>expdiff);
214done:
215
216 //TODO: find first nonzero digit and shift result and detect possibly underflow
217 while ((exp1>0)&&(!(mant1&0x20000000))) {
218 exp1--;
219 mant1 <<= 1;
220 if(mant1 == 0) {
221 /* Realy is it an underflow? ... */
222 /* TODO: fix underflow */
223 };
224 };
225
226 //rounding - if first bit after mantisa is set then round up
227 mant1 += 0x20;
228
229 if (mant1&0x40000000) {
230 ++exp1;
231 mant1>>=1;
232 };
233
234 result.parts.mantisa = ((mant1&(~0x20000000))>>6); /*Clear hidden bit and shift */
235 result.parts.exp = exp1;
236
237 return result;
238};
239
240/** Multiply two 32 bit float numbers
241 *
242 */
243float32 mulFloat32(float32 a, float32 b)
244{
245 float32 result;
246 __u64 mant1, mant2;
247 __s32 exp;
248
249 result.parts.sign = a.parts.sign ^ b.parts.sign;
250
251 if ((isFloat32NaN(a))||(isFloat32NaN(b))) {
252 /* TODO: fix SigNaNs */
253 if (isFloat32SigNaN(a)) {
254 result.parts.mantisa = a.parts.mantisa;
255 result.parts.exp = a.parts.exp;
256 return result;
257 };
258 if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
259 result.parts.mantisa = b.parts.mantisa;
260 result.parts.exp = b.parts.exp;
261 return result;
262 };
263 /* set NaN as result */
264 result.parts.mantisa = 0x1;
265 result.parts.exp = 0xFF;
266 return result;
267 };
268
269 if (isFloat32Infinity(a)) {
270 if (isFloat32Zero(b)) {
271 /* FIXME: zero * infinity */
272 result.parts.mantisa = 0x1;
273 result.parts.exp = 0xFF;
274 return result;
275 }
276 result.parts.mantisa = a.parts.mantisa;
277 result.parts.exp = a.parts.exp;
278 return result;
279 }
280
281 if (isFloat32Infinity(b)) {
282 if (isFloat32Zero(a)) {
283 /* FIXME: zero * infinity */
284 result.parts.mantisa = 0x1;
285 result.parts.exp = 0xFF;
286 return result;
287 }
288 result.parts.mantisa = b.parts.mantisa;
289 result.parts.exp = b.parts.exp;
290 return result;
291 }
292
293 /* exp is signed so we can easy detect underflow */
294 exp = a.parts.exp + b.parts.exp;
295 exp -= FLOAT32_BIAS;
296
297 if (exp >= 0xFF ) {
298 /* FIXME: overflow */
299 /* set infinity as result */
300 result.parts.mantisa = 0x0;
301 result.parts.exp = 0xFF;
302 return result;
303 };
304
305 if (exp < 0) {
306 /* FIXME: underflow */
307 /* return signed zero */
308 result.parts.mantisa = 0x0;
309 result.parts.exp = 0x0;
310 return result;
311 };
312
313 mant1 = a.parts.mantisa;
314 if (a.parts.exp>0) {
315 mant1 |= 0x800000;
316 } else {
317 ++exp;
318 };
319
320 mant2 = b.parts.mantisa;
321 if (b.parts.exp>0) {
322 mant2 |= 0x800000;
323 } else {
324 ++exp;
325 };
326
327 mant1 <<= 1; /* one bit space for rounding */
328
329 mant1 = mant1 * mant2;
330/* round and return */
331
332 while ((exp < 0xFF )&&(mant1 > 0x1FFFFFF )) { /* 0xFFFFFF is 23 bits of mantisa + one more for hidden bit (all shifted 1 bit left)*/
333 ++exp;
334 mant1 >>= 1;
335 };
336
337 /* rounding */
338 //++mant1; /* FIXME: not works - without it is ok */
339 mant1 >>= 1; /* shift off rounding space */
340
341 if ((exp < 0xFF )&&(mant1 > 0xFFFFFF )) {
342 ++exp;
343 mant1 >>= 1;
344 };
345
346 if (exp >= 0xFF ) {
347 /* TODO: fix overflow */
348 /* return infinity*/
349 result.parts.exp = 0xFF;
350 result.parts.mantisa = 0x0;
351 return result;
352 }
353
354 exp -= FLOAT32_MANTISA_SIZE;
355
356 if (exp <= FLOAT32_MANTISA_SIZE) {
357 /* denormalized number */
358 mant1 >>= 1; /* denormalize */
359 while ((mant1 > 0) && (exp < 0)) {
360 mant1 >>= 1;
361 ++exp;
362 };
363 if (mant1 == 0) {
364 /* FIXME : underflow */
365 result.parts.exp = 0;
366 result.parts.mantisa = 0;
367 return result;
368 };
369 };
370 result.parts.exp = exp;
371 result.parts.mantisa = mant1 & 0x7FFFFF;
372
373 return result;
374
375};
376
377
Note: See TracBrowser for help on using the repository browser.