← Back to branch summary

~galfy/helenos/bird-port-mainline

« back to all changes in this revision

Viewing changes to uspace/lib/softfloat/generic/div.c

  • Committer: Martin Decky
  • Date: 2009-08-04 11:19:19 UTC
  • Revision ID: martin@uranus.dsrg.hide.ms.mff.cuni.cz-20090804111919-evyclddlr3v5lhmp
Initial import

Show diffs side-by-side

added added

removed removed

Lines of Context:
uspace/lib/softfloat/generic/div.c
 
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
/** @addtogroup softfloat       
 
30
 * @{
 
31
 */
 
32
/** @file
 
33
 */
 
34
 
 
35
#include<sftypes.h>
 
36
#include<add.h>
 
37
#include<div.h>
 
38
#include<comparison.h>
 
39
#include<mul.h>
 
40
#include<common.h>
 
41
 
 
42
 
 
43
float32 divFloat32(float32 a, float32 b) 
 
44
{
 
45
        float32 result;
 
46
        int32_t aexp, bexp, cexp;
 
47
        uint64_t afrac, bfrac, cfrac;
 
48
        
 
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
        
 
153
        /* find first nonzero digit and shift result and detect possibly underflow */
 
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 {
 
189
                result.parts.exp = (uint32_t)cexp;
 
190
        }
 
191
        
 
192
        result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); 
 
193
        
 
194
        return result;  
 
195
}
 
196
 
 
197
float64 divFloat64(float64 a, float64 b) 
 
198
{
 
199
        float64 result;
 
200
        int64_t aexp, bexp, cexp;
 
201
        uint64_t afrac, bfrac, cfrac; 
 
202
        uint64_t remlo, remhi;
 
203
        
 
204
        result.parts.sign = a.parts.sign ^ b.parts.sign;
 
205
        
 
206
        if (isFloat64NaN(a)) {
 
207
                
 
208
                if (isFloat64SigNaN(b)) {
 
209
                        /*FIXME: SigNaN*/
 
210
                        return b;
 
211
                }
 
212
                
 
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)) {
 
229
                if (isFloat64Infinity(b) || isFloat64Zero(b)) {
 
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) {
 
274
                        result.parts.exp = 0;
 
275
                        result.parts.fraction = 0;
 
276
                        return result;
 
277
                }
 
278
                /* normalize it*/
 
279
                
 
280
                aexp++;
 
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) {
 
289
                bexp++;
 
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
                
 
315
                while ((int64_t) remhi < 0) {
 
316
                        cfrac--;
 
317
                        remlo += bfrac;
 
318
                        remhi += ( remlo < bfrac );
 
319
                }
 
320
                cfrac |= ( remlo != 0 );
 
321
        }
 
322
        
 
323
        /* round and shift */
 
324
        result = finishFloat64(cexp, cfrac, result.parts.sign);
 
325
        return result;
 
326
 
 
327
}
 
328
 
 
329
uint64_t divFloat64estim(uint64_t a, uint64_t b)
 
330
{
 
331
        uint64_t bhi;
 
332
        uint64_t remhi, remlo;
 
333
        uint64_t result;
 
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;
 
347
        while ( (int64_t) remhi < 0 ) {
 
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
 
 
363
/** @}
 
364
 */