Subversion Repositories HelenOS

Rev

Rev 2071 | Details | Compare with Previous | Last modification | View Log | RSS feed

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