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

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since a47d49f was 00acd66, checked in by Jakub Jermar <jakub@…>, 18 years ago

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