Subversion Repositories HelenOS

Rev

Rev 804 | Rev 829 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 804 Rev 828
Line 26... Line 26...
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;
Line 138... Line 140...
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;
Line 183... Line 184...
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
}
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
}
-
 
393