Subversion Repositories HelenOS-historic

Rev

Rev 829 | Rev 1657 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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