Subversion Repositories HelenOS

Rev

Rev 804 | Rev 829 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 804 Rev 828
1
/*
1
/*
2
 * Copyright (C) 2005 Josef Cejka
2
 * Copyright (C) 2005 Josef Cejka
3
 * All rights reserved.
3
 * All rights reserved.
4
 *
4
 *
5
 * Redistribution and use in source and binary forms, with or without
5
 * Redistribution and use in source and binary forms, with or without
6
 * modification, are permitted provided that the following conditions
6
 * modification, are permitted provided that the following conditions
7
 * are met:
7
 * are met:
8
 *
8
 *
9
 * - Redistributions of source code must retain the above copyright
9
 * - Redistributions of source code must retain the above copyright
10
 *   notice, this list of conditions and the following disclaimer.
10
 *   notice, this list of conditions and the following disclaimer.
11
 * - Redistributions in binary form must reproduce the above copyright
11
 * - Redistributions in binary form must reproduce the above copyright
12
 *   notice, this list of conditions and the following disclaimer in the
12
 *   notice, this list of conditions and the following disclaimer in the
13
 *   documentation and/or other materials provided with the distribution.
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
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.
15
 *   derived from this software without specific prior written permission.
16
 *
16
 *
17
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
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
18
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
19
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
19
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
20
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
21
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
22
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
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
23
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
24
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
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
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.
26
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27
 */
27
 */
28
 
28
 
29
#include<sftypes.h>
29
#include<sftypes.h>
30
#include<add.h>
30
#include<add.h>
-
 
31
#include<div.h>
31
#include<comparison.h>
32
#include<comparison.h>
-
 
33
#include<mul.h>
32
 
34
 
33
float32 divFloat32(float32 a, float32 b)
35
float32 divFloat32(float32 a, float32 b)
34
{
36
{
35
    float32 result;
37
    float32 result;
36
    __s32 aexp, bexp, cexp;
38
    __s32 aexp, bexp, cexp;
37
    __u64 afrac, bfrac, cfrac;
39
    __u64 afrac, bfrac, cfrac;
38
   
40
   
39
    result.parts.sign = a.parts.sign ^ b.parts.sign;
41
    result.parts.sign = a.parts.sign ^ b.parts.sign;
40
   
42
   
41
    if (isFloat32NaN(a)) {
43
    if (isFloat32NaN(a)) {
42
        if (isFloat32SigNaN(a)) {
44
        if (isFloat32SigNaN(a)) {
43
            /*FIXME: SigNaN*/
45
            /*FIXME: SigNaN*/
44
        }
46
        }
45
        /*NaN*/
47
        /*NaN*/
46
        return a;
48
        return a;
47
    }
49
    }
48
   
50
   
49
    if (isFloat32NaN(b)) {
51
    if (isFloat32NaN(b)) {
50
        if (isFloat32SigNaN(b)) {
52
        if (isFloat32SigNaN(b)) {
51
            /*FIXME: SigNaN*/
53
            /*FIXME: SigNaN*/
52
        }
54
        }
53
        /*NaN*/
55
        /*NaN*/
54
        return b;
56
        return b;
55
    }
57
    }
56
   
58
   
57
    if (isFloat32Infinity(a)) {
59
    if (isFloat32Infinity(a)) {
58
        if (isFloat32Infinity(b)) {
60
        if (isFloat32Infinity(b)) {
59
            /*FIXME: inf / inf */
61
            /*FIXME: inf / inf */
60
            result.binary = FLOAT32_NAN;
62
            result.binary = FLOAT32_NAN;
61
            return result;
63
            return result;
62
        }
64
        }
63
        /* inf / num */
65
        /* inf / num */
64
        result.parts.exp = a.parts.exp;
66
        result.parts.exp = a.parts.exp;
65
        result.parts.fraction = a.parts.fraction;
67
        result.parts.fraction = a.parts.fraction;
66
        return result;
68
        return result;
67
    }
69
    }
68
 
70
 
69
    if (isFloat32Infinity(b)) {
71
    if (isFloat32Infinity(b)) {
70
        if (isFloat32Zero(a)) {
72
        if (isFloat32Zero(a)) {
71
            /* FIXME 0 / inf */
73
            /* FIXME 0 / inf */
72
            result.parts.exp = 0;
74
            result.parts.exp = 0;
73
            result.parts.fraction = 0;
75
            result.parts.fraction = 0;
74
            return result;
76
            return result;
75
        }
77
        }
76
        /* FIXME: num / inf*/
78
        /* FIXME: num / inf*/
77
        result.parts.exp = 0;
79
        result.parts.exp = 0;
78
        result.parts.fraction = 0;
80
        result.parts.fraction = 0;
79
        return result;
81
        return result;
80
    }
82
    }
81
   
83
   
82
    if (isFloat32Zero(b)) {
84
    if (isFloat32Zero(b)) {
83
        if (isFloat32Zero(a)) {
85
        if (isFloat32Zero(a)) {
84
            /*FIXME: 0 / 0*/
86
            /*FIXME: 0 / 0*/
85
            result.binary = FLOAT32_NAN;
87
            result.binary = FLOAT32_NAN;
86
            return result;
88
            return result;
87
        }
89
        }
88
        /* FIXME: division by zero */
90
        /* FIXME: division by zero */
89
        result.parts.exp = 0;
91
        result.parts.exp = 0;
90
        result.parts.fraction = 0;
92
        result.parts.fraction = 0;
91
        return result;
93
        return result;
92
    }
94
    }
93
 
95
 
94
   
96
   
95
    afrac = a.parts.fraction;
97
    afrac = a.parts.fraction;
96
    aexp = a.parts.exp;
98
    aexp = a.parts.exp;
97
    bfrac = b.parts.fraction;
99
    bfrac = b.parts.fraction;
98
    bexp = b.parts.exp;
100
    bexp = b.parts.exp;
99
   
101
   
100
    /* denormalized numbers */
102
    /* denormalized numbers */
101
    if (aexp == 0) {
103
    if (aexp == 0) {
102
        if (afrac == 0) {
104
        if (afrac == 0) {
103
        result.parts.exp = 0;
105
        result.parts.exp = 0;
104
        result.parts.fraction = 0;
106
        result.parts.fraction = 0;
105
        return result;
107
        return result;
106
        }
108
        }
107
        /* normalize it*/
109
        /* normalize it*/
108
       
110
       
109
        afrac <<= 1;
111
        afrac <<= 1;
110
            /* afrac is nonzero => it must stop */ 
112
            /* afrac is nonzero => it must stop */ 
111
        while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
113
        while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
112
            afrac <<= 1;
114
            afrac <<= 1;
113
            aexp--;
115
            aexp--;
114
        }
116
        }
115
    }
117
    }
116
 
118
 
117
    if (bexp == 0) {
119
    if (bexp == 0) {
118
        bfrac <<= 1;
120
        bfrac <<= 1;
119
            /* bfrac is nonzero => it must stop */ 
121
            /* bfrac is nonzero => it must stop */ 
120
        while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
122
        while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
121
            bfrac <<= 1;
123
            bfrac <<= 1;
122
            bexp--;
124
            bexp--;
123
        }
125
        }
124
    }
126
    }
125
 
127
 
126
    afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
128
    afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
127
    bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
129
    bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
128
 
130
 
129
    if ( bfrac <= (afrac << 1) ) {
131
    if ( bfrac <= (afrac << 1) ) {
130
        afrac >>= 1;
132
        afrac >>= 1;
131
        aexp++;
133
        aexp++;
132
    }
134
    }
133
   
135
   
134
    cexp = aexp - bexp + FLOAT32_BIAS - 2;
136
    cexp = aexp - bexp + FLOAT32_BIAS - 2;
135
   
137
   
136
    cfrac = (afrac << 32) / bfrac;
138
    cfrac = (afrac << 32) / bfrac;
137
    if ((  cfrac & 0x3F ) == 0) {
139
    if ((  cfrac & 0x3F ) == 0) {
138
        cfrac |= ( bfrac * cfrac != afrac << 32 );
140
        cfrac |= ( bfrac * cfrac != afrac << 32 );
139
    }
141
    }
140
   
142
   
141
    /* pack and round */
143
    /* pack and round */
142
   
144
   
143
    /* TODO: find first nonzero digit and shift result and detect possibly underflow */
145
    /* find first nonzero digit and shift result and detect possibly underflow */
144
    while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
146
    while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
145
        cexp--;
147
        cexp--;
146
        cfrac <<= 1;
148
        cfrac <<= 1;
147
            /* TODO: fix underflow */
149
            /* TODO: fix underflow */
148
    };
150
    };
149
   
151
   
150
   
-
 
151
    cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
152
    cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
152
   
153
   
153
    if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
154
    if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
154
        ++cexp;
155
        ++cexp;
155
        cfrac >>= 1;
156
        cfrac >>= 1;
156
        }  
157
        }  
157
 
158
 
158
    /* check overflow */
159
    /* check overflow */
159
    if (cexp >= FLOAT32_MAX_EXPONENT ) {
160
    if (cexp >= FLOAT32_MAX_EXPONENT ) {
160
        /* FIXME: overflow, return infinity */
161
        /* FIXME: overflow, return infinity */
161
        result.parts.exp = FLOAT32_MAX_EXPONENT;
162
        result.parts.exp = FLOAT32_MAX_EXPONENT;
162
        result.parts.fraction = 0;
163
        result.parts.fraction = 0;
163
        return result;
164
        return result;
164
    }
165
    }
165
 
166
 
166
    if (cexp < 0) {
167
    if (cexp < 0) {
167
        /* FIXME: underflow */
168
        /* FIXME: underflow */
168
        result.parts.exp = 0;
169
        result.parts.exp = 0;
169
        if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
170
        if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
170
            result.parts.fraction = 0;
171
            result.parts.fraction = 0;
171
            return result;
172
            return result;
172
        }
173
        }
173
        cfrac >>= 1;
174
        cfrac >>= 1;
174
        while (cexp < 0) {
175
        while (cexp < 0) {
175
            cexp ++;
176
            cexp ++;
176
            cfrac >>= 1;
177
            cfrac >>= 1;
177
        }
178
        }
178
       
179
       
179
    } else {
180
    } else {
180
        result.parts.exp = (__u32)cexp;
181
        result.parts.exp = (__u32)cexp;
181
    }
182
    }
182
   
183
   
183
    result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
184
    result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
184
   
185
   
185
    return result; 
186
    return result; 
186
}
187
}
-
 
188
 
-
 
189
float64 divFloat64(float64 a, float64 b)
-
 
190
{
-
 
191
    float64 result;
-
 
192
    __s32 aexp, bexp, cexp;
-
 
193
    __u64 afrac, bfrac, cfrac;
-
 
194
    __u64 remlo, remhi;
-
 
195
   
-
 
196
    result.parts.sign = a.parts.sign ^ b.parts.sign;
-
 
197
   
-
 
198
    if (isFloat64NaN(a)) {
-
 
199
        if (isFloat64SigNaN(a)) {
-
 
200
            /*FIXME: SigNaN*/
-
 
201
        }
-
 
202
        /*NaN*/
-
 
203
        return a;
-
 
204
    }
-
 
205
   
-
 
206
    if (isFloat64NaN(b)) {
-
 
207
        if (isFloat64SigNaN(b)) {
-
 
208
            /*FIXME: SigNaN*/
-
 
209
        }
-
 
210
        /*NaN*/
-
 
211
        return b;
-
 
212
    }
-
 
213
   
-
 
214
    if (isFloat64Infinity(a)) {
-
 
215
        if (isFloat64Infinity(b)) {
-
 
216
            /*FIXME: inf / inf */
-
 
217
            result.binary = FLOAT64_NAN;
-
 
218
            return result;
-
 
219
        }
-
 
220
        /* inf / num */
-
 
221
        result.parts.exp = a.parts.exp;
-
 
222
        result.parts.fraction = a.parts.fraction;
-
 
223
        return result;
-
 
224
    }
-
 
225
 
-
 
226
    if (isFloat64Infinity(b)) {
-
 
227
        if (isFloat64Zero(a)) {
-
 
228
            /* FIXME 0 / inf */
-
 
229
            result.parts.exp = 0;
-
 
230
            result.parts.fraction = 0;
-
 
231
            return result;
-
 
232
        }
-
 
233
        /* FIXME: num / inf*/
-
 
234
        result.parts.exp = 0;
-
 
235
        result.parts.fraction = 0;
-
 
236
        return result;
-
 
237
    }
-
 
238
   
-
 
239
    if (isFloat64Zero(b)) {
-
 
240
        if (isFloat64Zero(a)) {
-
 
241
            /*FIXME: 0 / 0*/
-
 
242
            result.binary = FLOAT64_NAN;
-
 
243
            return result;
-
 
244
        }
-
 
245
        /* FIXME: division by zero */
-
 
246
        result.parts.exp = 0;
-
 
247
        result.parts.fraction = 0;
-
 
248
        return result;
-
 
249
    }
-
 
250
 
-
 
251
   
-
 
252
    afrac = a.parts.fraction;
-
 
253
    aexp = a.parts.exp;
-
 
254
    bfrac = b.parts.fraction;
-
 
255
    bexp = b.parts.exp;
-
 
256
   
-
 
257
    /* denormalized numbers */
-
 
258
    if (aexp == 0) {
-
 
259
        if (afrac == 0) {
-
 
260
        result.parts.exp = 0;
-
 
261
        result.parts.fraction = 0;
-
 
262
        return result;
-
 
263
        }
-
 
264
        /* normalize it*/
-
 
265
       
-
 
266
        afrac <<= 1;
-
 
267
            /* afrac is nonzero => it must stop */ 
-
 
268
        while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
-
 
269
            afrac <<= 1;
-
 
270
            aexp--;
-
 
271
        }
-
 
272
    }
-
 
273
 
-
 
274
    if (bexp == 0) {
-
 
275
        bfrac <<= 1;
-
 
276
            /* bfrac is nonzero => it must stop */ 
-
 
277
        while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
-
 
278
            bfrac <<= 1;
-
 
279
            bexp--;
-
 
280
        }
-
 
281
    }
-
 
282
 
-
 
283
    afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
-
 
284
    bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
-
 
285
 
-
 
286
    if ( bfrac <= (afrac << 1) ) {
-
 
287
        afrac >>= 1;
-
 
288
        aexp++;
-
 
289
    }
-
 
290
   
-
 
291
    cexp = aexp - bexp + FLOAT64_BIAS - 2;
-
 
292
   
-
 
293
    cfrac = divFloat64estim(afrac, bfrac);
-
 
294
   
-
 
295
    if ((  cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
-
 
296
        mul64integers( bfrac, cfrac, &remlo, &remhi);
-
 
297
        /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ 
-
 
298
        remhi = afrac - remhi - ( remlo > 0);
-
 
299
        remlo = - remlo;
-
 
300
       
-
 
301
        while ((__s64) remhi < 0) {
-
 
302
            cfrac--;
-
 
303
            remlo += bfrac;
-
 
304
            remhi += ( remlo < bfrac );
-
 
305
        }
-
 
306
        cfrac |= ( remlo != 0 );
-
 
307
    }
-
 
308
   
-
 
309
    /* pack and round */
-
 
310
   
-
 
311
    /* find first nonzero digit and shift result and detect possibly underflow */
-
 
312
    while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ) )))) {
-
 
313
        cexp--;
-
 
314
        cfrac <<= 1;
-
 
315
            /* TODO: fix underflow */
-
 
316
    };
-
 
317
   
-
 
318
   
-
 
319
    cfrac >>= 1;
-
 
320
    ++cexp;
-
 
321
    cfrac += (0x1 << (64 - FLOAT64_FRACTION_SIZE - 3));
-
 
322
 
-
 
323
    if (cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ))) {
-
 
324
        ++cexp;
-
 
325
        cfrac >>= 1;
-
 
326
        }  
-
 
327
 
-
 
328
    /* check overflow */
-
 
329
    if (cexp >= FLOAT64_MAX_EXPONENT ) {
-
 
330
        /* FIXME: overflow, return infinity */
-
 
331
        result.parts.exp = FLOAT64_MAX_EXPONENT;
-
 
332
        result.parts.fraction = 0;
-
 
333
        return result;
-
 
334
    }
-
 
335
 
-
 
336
    if (cexp < 0) {
-
 
337
        /* FIXME: underflow */
-
 
338
        result.parts.exp = 0;
-
 
339
        if ((cexp + FLOAT64_FRACTION_SIZE) < 0) {
-
 
340
            result.parts.fraction = 0;
-
 
341
            return result;
-
 
342
        }
-
 
343
        cfrac >>= 1;
-
 
344
        while (cexp < 0) {
-
 
345
            cexp ++;
-
 
346
            cfrac >>= 1;
-
 
347
        }
-
 
348
        return result;
-
 
349
       
-
 
350
    } else {
-
 
351
        cexp ++; /*normalized*/
-
 
352
        result.parts.exp = (__u32)cexp;
-
 
353
    }
-
 
354
   
-
 
355
    result.parts.fraction = ((cfrac >>(64 - FLOAT64_FRACTION_SIZE - 2 ) ) & (~FLOAT64_HIDDEN_BIT_MASK));
-
 
356
   
-
 
357
    return result; 
-
 
358
}
-
 
359
 
-
 
360
__u64 divFloat64estim(__u64 a, __u64 b)
-
 
361
{
-
 
362
    __u64 bhi;
-
 
363
    __u64 remhi, remlo;
-
 
364
    __u64 result;
-
 
365
   
-
 
366
    if ( b <= a ) {
-
 
367
        return 0xFFFFFFFFFFFFFFFFull;
-
 
368
    }
-
 
369
   
-
 
370
    bhi = b >> 32;
-
 
371
    result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
-
 
372
    mul64integers(b, result, &remlo, &remhi);
-
 
373
   
-
 
374
    remhi = a - remhi - (remlo > 0);
-
 
375
    remlo = - remlo;
-
 
376
 
-
 
377
    b <<= 32;
-
 
378
    while ( (__s64) remhi < 0 ) {
-
 
379
            result -= 0x1ll << 32; 
-
 
380
            remlo += b;
-
 
381
            remhi += bhi + ( remlo < b );
-
 
382
        }
-
 
383
    remhi = (remhi << 32) | (remlo >> 32);
-
 
384
    if (( bhi << 32) <= remhi) {
-
 
385
        result |= 0xFFFFFFFF;
-
 
386
    } else {
-
 
387
        result |= remhi / bhi;
-
 
388
    }
-
 
389
   
-
 
390
   
-
 
391
    return result;
-
 
392
}
187
 
393
 
188
 
394