source: mainline/softfloat/generic/add.c@ d3ca210

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

Fixed most problems with 64bit arithmetic (but division is still buggy).
  • Property mode set to 100644
File size: 5.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#include<sftypes.h>
30#include<add.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,frac1, frac2;
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 == FLOAT32_MAX_EXPONENT) {
51 return b;
52 }
53
54 frac1 = b.parts.fraction;
55 exp1 = b.parts.exp;
56 frac2 = a.parts.fraction;
57 exp2 = a.parts.exp;
58 expdiff *= -1;
59 } else {
60 if ((isFloat32NaN(a)) || (isFloat32NaN(b))) {
61 /* TODO: fix SigNaN */
62 if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
63 };
64 return (isFloat32NaN(a)?a:b);
65 };
66
67 if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
68 return a;
69 }
70
71 frac1 = a.parts.fraction;
72 exp1 = a.parts.exp;
73 frac2 = b.parts.fraction;
74 exp2 = b.parts.exp;
75 };
76
77 if (exp1 == 0) {
78 /* both are denormalized */
79 frac1 += frac2;
80 if (frac1 & FLOAT32_HIDDEN_BIT_MASK ) {
81 /* result is not denormalized */
82 a.parts.exp = 1;
83 };
84 a.parts.fraction = frac1;
85 return a;
86 };
87
88 frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
89
90 if (exp2 == 0) {
91 /* second operand is denormalized */
92 --expdiff;
93 } else {
94 /* add hidden bit to second operand */
95 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
96 };
97
98 /* create some space for rounding */
99 frac1 <<= 6;
100 frac2 <<= 6;
101
102 if (expdiff < (FLOAT32_FRACTION_SIZE + 2) ) {
103 frac2 >>= expdiff;
104 frac1 += frac2;
105 } else {
106 a.parts.exp = exp1;
107 a.parts.fraction = (frac1 >> 6) & (~(FLOAT32_HIDDEN_BIT_MASK));
108 return a;
109 }
110
111 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
112 ++exp1;
113 frac1 >>= 1;
114 };
115
116 /* rounding - if first bit after fraction is set then round up */
117 frac1 += (0x1 << 5);
118
119 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
120 /* rounding overflow */
121 ++exp1;
122 frac1 >>= 1;
123 };
124
125
126 if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
127 /* overflow - set infinity as result */
128 a.parts.exp = FLOAT32_MAX_EXPONENT;
129 a.parts.fraction = 0;
130 return a;
131 }
132
133 a.parts.exp = exp1;
134
135 /*Clear hidden bit and shift */
136 a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
137 return a;
138}
139
140
141/** Add two Float64 numbers with same signs
142 */
143float64 addFloat64(float64 a, float64 b)
144{
145 int expdiff;
146 __u32 exp1, exp2;
147 __u64 frac1, frac2;
148
149 expdiff = ((int )a.parts.exp) - b.parts.exp;
150 if (expdiff < 0) {
151 if (isFloat64NaN(b)) {
152 /* TODO: fix SigNaN */
153 if (isFloat64SigNaN(b)) {
154 };
155
156 return b;
157 };
158
159 /* b is infinity and a not */
160 if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
161 return b;
162 }
163
164 frac1 = b.parts.fraction;
165 exp1 = b.parts.exp;
166 frac2 = a.parts.fraction;
167 exp2 = a.parts.exp;
168 expdiff *= -1;
169 } else {
170 if (isFloat64NaN(a)) {
171 /* TODO: fix SigNaN */
172 if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
173 };
174 return a;
175 };
176
177 /* a is infinity and b not */
178 if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
179 return a;
180 }
181
182 frac1 = a.parts.fraction;
183 exp1 = a.parts.exp;
184 frac2 = b.parts.fraction;
185 exp2 = b.parts.exp;
186 };
187
188 if (exp1 == 0) {
189 /* both are denormalized */
190 frac1 += frac2;
191 if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
192 /* result is not denormalized */
193 a.parts.exp = 1;
194 };
195 a.parts.fraction = frac1;
196 return a;
197 };
198
199 /* add hidden bit - frac1 is sure not denormalized */
200 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
201
202 /* second operand ... */
203 if (exp2 == 0) {
204 /* ... is denormalized */
205 --expdiff;
206 } else {
207 /* is not denormalized */
208 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
209 };
210
211 /* create some space for rounding */
212 frac1 <<= 6;
213 frac2 <<= 6;
214
215 if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
216 frac2 >>= expdiff;
217 frac1 += frac2;
218 } else {
219 a.parts.exp = exp1;
220 a.parts.fraction = (frac1 >> 6) & (~(FLOAT64_HIDDEN_BIT_MASK));
221 return a;
222 }
223
224 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
225 ++exp1;
226 frac1 >>= 1;
227 };
228
229 /* rounding - if first bit after fraction is set then round up */
230 frac1 += (0x1 << 5);
231
232 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
233 /* rounding overflow */
234 ++exp1;
235 frac1 >>= 1;
236 };
237
238 if ((exp1 == FLOAT64_MAX_EXPONENT ) || (exp2 > exp1)) {
239 /* overflow - set infinity as result */
240 a.parts.exp = FLOAT64_MAX_EXPONENT;
241 a.parts.fraction = 0;
242 return a;
243 }
244
245 a.parts.exp = exp1;
246 /*Clear hidden bit and shift */
247 a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
248
249 return a;
250}
251
252
Note: See TracBrowser for help on using the repository browser.