source: mainline/uspace/lib/softfloat/generic/div.c@ 98e4507

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

softfloat: slightly improve coding style (no change in functionality)
  • Property mode set to 100644
File size: 7.7 KB
RevLine 
[b5440cf]1/*
[df4ed85]2 * Copyright (c) 2005 Josef Cejka
[b5440cf]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
[750636a]29/** @addtogroup softfloat
[846848a6]30 * @{
31 */
32/** @file
33 */
34
[750636a]35#include <sftypes.h>
36#include <add.h>
37#include <div.h>
38#include <comparison.h>
39#include <mul.h>
40#include <common.h>
[b5440cf]41
[12c6f2d]42float32 divFloat32(float32 a, float32 b)
43{
[1266543]44 float32 result;
[aa59fa0]45 int32_t aexp, bexp, cexp;
46 uint64_t afrac, bfrac, cfrac;
[12c6f2d]47
[1266543]48 result.parts.sign = a.parts.sign ^ b.parts.sign;
49
50 if (isFloat32NaN(a)) {
51 if (isFloat32SigNaN(a)) {
52 /*FIXME: SigNaN*/
53 }
54 /*NaN*/
55 return a;
56 }
57
58 if (isFloat32NaN(b)) {
59 if (isFloat32SigNaN(b)) {
60 /*FIXME: SigNaN*/
61 }
62 /*NaN*/
63 return b;
64 }
65
66 if (isFloat32Infinity(a)) {
67 if (isFloat32Infinity(b)) {
68 /*FIXME: inf / inf */
69 result.binary = FLOAT32_NAN;
70 return result;
71 }
72 /* inf / num */
73 result.parts.exp = a.parts.exp;
74 result.parts.fraction = a.parts.fraction;
75 return result;
76 }
77
78 if (isFloat32Infinity(b)) {
79 if (isFloat32Zero(a)) {
80 /* FIXME 0 / inf */
81 result.parts.exp = 0;
82 result.parts.fraction = 0;
83 return result;
84 }
85 /* FIXME: num / inf*/
86 result.parts.exp = 0;
87 result.parts.fraction = 0;
88 return result;
89 }
90
91 if (isFloat32Zero(b)) {
92 if (isFloat32Zero(a)) {
93 /*FIXME: 0 / 0*/
94 result.binary = FLOAT32_NAN;
95 return result;
96 }
97 /* FIXME: division by zero */
98 result.parts.exp = 0;
99 result.parts.fraction = 0;
100 return result;
101 }
102
103
104 afrac = a.parts.fraction;
105 aexp = a.parts.exp;
106 bfrac = b.parts.fraction;
107 bexp = b.parts.exp;
108
109 /* denormalized numbers */
110 if (aexp == 0) {
111 if (afrac == 0) {
112 result.parts.exp = 0;
113 result.parts.fraction = 0;
114 return result;
115 }
116 /* normalize it*/
117
118 afrac <<= 1;
119 /* afrac is nonzero => it must stop */
120 while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
121 afrac <<= 1;
122 aexp--;
123 }
124 }
125
126 if (bexp == 0) {
127 bfrac <<= 1;
128 /* bfrac is nonzero => it must stop */
129 while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
130 bfrac <<= 1;
131 bexp--;
132 }
133 }
134
135 afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
136 bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
137
138 if ( bfrac <= (afrac << 1) ) {
139 afrac >>= 1;
140 aexp++;
141 }
142
143 cexp = aexp - bexp + FLOAT32_BIAS - 2;
144
145 cfrac = (afrac << 32) / bfrac;
146 if (( cfrac & 0x3F ) == 0) {
147 cfrac |= ( bfrac * cfrac != afrac << 32 );
148 }
149
150 /* pack and round */
151
[e6a40ac]152 /* find first nonzero digit and shift result and detect possibly underflow */
[1266543]153 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
154 cexp--;
155 cfrac <<= 1;
156 /* TODO: fix underflow */
157 };
158
159 cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
160
161 if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
162 ++cexp;
163 cfrac >>= 1;
164 }
165
166 /* check overflow */
167 if (cexp >= FLOAT32_MAX_EXPONENT ) {
168 /* FIXME: overflow, return infinity */
169 result.parts.exp = FLOAT32_MAX_EXPONENT;
170 result.parts.fraction = 0;
171 return result;
172 }
173
174 if (cexp < 0) {
175 /* FIXME: underflow */
176 result.parts.exp = 0;
177 if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
178 result.parts.fraction = 0;
179 return result;
180 }
181 cfrac >>= 1;
182 while (cexp < 0) {
183 cexp ++;
184 cfrac >>= 1;
185 }
186
187 } else {
[aa59fa0]188 result.parts.exp = (uint32_t)cexp;
[1266543]189 }
190
191 result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
192
193 return result;
[12c6f2d]194}
[b5440cf]195
[e6a40ac]196float64 divFloat64(float64 a, float64 b)
197{
198 float64 result;
[aa59fa0]199 int64_t aexp, bexp, cexp;
200 uint64_t afrac, bfrac, cfrac;
201 uint64_t remlo, remhi;
[e6a40ac]202
203 result.parts.sign = a.parts.sign ^ b.parts.sign;
204
205 if (isFloat64NaN(a)) {
[f1f95f2]206
207 if (isFloat64SigNaN(b)) {
208 /*FIXME: SigNaN*/
209 return b;
210 }
211
[e6a40ac]212 if (isFloat64SigNaN(a)) {
213 /*FIXME: SigNaN*/
214 }
215 /*NaN*/
216 return a;
217 }
218
219 if (isFloat64NaN(b)) {
220 if (isFloat64SigNaN(b)) {
221 /*FIXME: SigNaN*/
222 }
223 /*NaN*/
224 return b;
225 }
226
227 if (isFloat64Infinity(a)) {
[f1f95f2]228 if (isFloat64Infinity(b) || isFloat64Zero(b)) {
[e6a40ac]229 /*FIXME: inf / inf */
230 result.binary = FLOAT64_NAN;
231 return result;
232 }
233 /* inf / num */
234 result.parts.exp = a.parts.exp;
235 result.parts.fraction = a.parts.fraction;
236 return result;
237 }
238
239 if (isFloat64Infinity(b)) {
240 if (isFloat64Zero(a)) {
241 /* FIXME 0 / inf */
242 result.parts.exp = 0;
243 result.parts.fraction = 0;
244 return result;
245 }
246 /* FIXME: num / inf*/
247 result.parts.exp = 0;
248 result.parts.fraction = 0;
249 return result;
250 }
251
252 if (isFloat64Zero(b)) {
253 if (isFloat64Zero(a)) {
254 /*FIXME: 0 / 0*/
255 result.binary = FLOAT64_NAN;
256 return result;
257 }
258 /* FIXME: division by zero */
259 result.parts.exp = 0;
260 result.parts.fraction = 0;
261 return result;
262 }
263
264
265 afrac = a.parts.fraction;
266 aexp = a.parts.exp;
267 bfrac = b.parts.fraction;
268 bexp = b.parts.exp;
269
270 /* denormalized numbers */
271 if (aexp == 0) {
272 if (afrac == 0) {
[f1f95f2]273 result.parts.exp = 0;
274 result.parts.fraction = 0;
275 return result;
[e6a40ac]276 }
277 /* normalize it*/
278
[f1f95f2]279 aexp++;
[e6a40ac]280 /* afrac is nonzero => it must stop */
281 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
282 afrac <<= 1;
283 aexp--;
284 }
285 }
286
287 if (bexp == 0) {
[f1f95f2]288 bexp++;
[e6a40ac]289 /* bfrac is nonzero => it must stop */
290 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
291 bfrac <<= 1;
292 bexp--;
293 }
294 }
295
296 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
297 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
298
299 if ( bfrac <= (afrac << 1) ) {
300 afrac >>= 1;
301 aexp++;
302 }
303
304 cexp = aexp - bexp + FLOAT64_BIAS - 2;
305
306 cfrac = divFloat64estim(afrac, bfrac);
307
308 if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
309 mul64integers( bfrac, cfrac, &remlo, &remhi);
310 /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/
311 remhi = afrac - remhi - ( remlo > 0);
312 remlo = - remlo;
313
[aa59fa0]314 while ((int64_t) remhi < 0) {
[e6a40ac]315 cfrac--;
316 remlo += bfrac;
317 remhi += ( remlo < bfrac );
318 }
319 cfrac |= ( remlo != 0 );
320 }
321
[e979fea]322 /* round and shift */
323 result = finishFloat64(cexp, cfrac, result.parts.sign);
324 return result;
[e6a40ac]325
326}
327
[aa59fa0]328uint64_t divFloat64estim(uint64_t a, uint64_t b)
[e6a40ac]329{
[aa59fa0]330 uint64_t bhi;
331 uint64_t remhi, remlo;
332 uint64_t result;
[e6a40ac]333
334 if ( b <= a ) {
335 return 0xFFFFFFFFFFFFFFFFull;
336 }
337
338 bhi = b >> 32;
339 result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
340 mul64integers(b, result, &remlo, &remhi);
341
342 remhi = a - remhi - (remlo > 0);
343 remlo = - remlo;
344
345 b <<= 32;
[aa59fa0]346 while ( (int64_t) remhi < 0 ) {
[e6a40ac]347 result -= 0x1ll << 32;
348 remlo += b;
349 remhi += bhi + ( remlo < b );
350 }
351 remhi = (remhi << 32) | (remlo >> 32);
352 if (( bhi << 32) <= remhi) {
353 result |= 0xFFFFFFFF;
354 } else {
355 result |= remhi / bhi;
356 }
357
358
359 return result;
360}
361
[231a60a]362/** @}
[846848a6]363 */
Note: See TracBrowser for help on using the repository browser.