Changeset 888207c9 in mainline for uspace/lib/softfloat/generic/div.c


Ignore:
Timestamp:
2011-09-07T22:19:24Z (13 years ago)
Author:
Jakub Jermar <jakub@…>
Branches:
lfn, master, serial, ticket/834-toolchain-update, topic/msim-upgrade, topic/simplify-dev-export
Children:
2aeec8ef
Parents:
038b289 (diff), 5081276 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge from lp:~helenos-posix/helenos/devel.

  • Petr Koupy's port of binutils.
  • Jiri Zarevucky's port of pcc.
  • Petr and Jiri's libposix library.
  • Changes needed to make it all work together.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • uspace/lib/softfloat/generic/div.c

    r038b289 r888207c9  
    11/*
    22 * Copyright (c) 2005 Josef Cejka
     3 * Copyright (c) 2011 Petr Koupy
    34 * All rights reserved.
    45 *
    … …  
    3031 * @{
    3132 */
    32 /** @file
     33/** @file Division functions.
    3334 */
    3435
    … …  
    4041#include <common.h>
    4142
     43/**
     44 * Divide two single-precision floats.
     45 *
     46 * @param a Nominator.
     47 * @param b Denominator.
     48 * @return Result of division.
     49 */
    4250float32 divFloat32(float32 a, float32 b)
    4351{
    … …  
    100108                return result;
    101109        }
    102 
    103110       
    104111        afrac = a.parts.fraction;
    … …  
    110117        if (aexp == 0) {
    111118                if (afrac == 0) {
    112                 result.parts.exp = 0;
    113                 result.parts.fraction = 0;
    114                 return result;
    115                 }
     119                        result.parts.exp = 0;
     120                        result.parts.fraction = 0;
     121                        return result;
     122                }
     123
    116124                /* normalize it*/
    117                
    118125                afrac <<= 1;
    119                         /* afrac is nonzero => it must stop */ 
    120                 while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
     126                /* afrac is nonzero => it must stop */ 
     127                while (!(afrac & FLOAT32_HIDDEN_BIT_MASK)) {
    121128                        afrac <<= 1;
    122129                        aexp--;
    … …  
    126133        if (bexp == 0) {
    127134                bfrac <<= 1;
    128                         /* bfrac is nonzero => it must stop */ 
    129                 while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
     135                /* bfrac is nonzero => it must stop */ 
     136                while (!(bfrac & FLOAT32_HIDDEN_BIT_MASK)) {
    130137                        bfrac <<= 1;
    131138                        bexp--;
    … …  
    133140        }
    134141
    135         afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
    136         bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
    137 
    138         if ( bfrac <= (afrac << 1) ) {
     142        afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE - 1);
     143        bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE);
     144
     145        if (bfrac <= (afrac << 1)) {
    139146                afrac >>= 1;
    140147                aexp++;
    … …  
    144151       
    145152        cfrac = (afrac << 32) / bfrac;
    146         if ((  cfrac & 0x3F ) == 0) {
    147                 cfrac |= ( bfrac * cfrac != afrac << 32 );
     153        if ((cfrac & 0x3F) == 0) {
     154                cfrac |= (bfrac * cfrac != afrac << 32);
    148155        }
    149156       
    … …  
    151158       
    152159        /* find first nonzero digit and shift result and detect possibly underflow */
    153         while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
     160        while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)))) {
    154161                cexp--;
    155162                cfrac <<= 1;
    156                         /* TODO: fix underflow */
    157         };
     163                /* TODO: fix underflow */
     164        }
    158165       
    159166        cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
    … …  
    162169                ++cexp;
    163170                cfrac >>= 1;
    164                 }       
     171        }       
    165172
    166173        /* check overflow */
    167         if (cexp >= FLOAT32_MAX_EXPONENT ) {
     174        if (cexp >= FLOAT32_MAX_EXPONENT) {
    168175                /* FIXME: overflow, return infinity */
    169176                result.parts.exp = FLOAT32_MAX_EXPONENT;
    … …  
    181188                cfrac >>= 1;
    182189                while (cexp < 0) {
    183                         cexp ++;
     190                        cexp++;
    184191                        cfrac >>= 1;
    185                 }
    186                
     192                }       
    187193        } else {
    188                 result.parts.exp = (uint32_t)cexp;
     194                result.parts.exp = (uint32_t) cexp;
    189195        }
    190196       
    … …  
    194200}
    195201
     202/**
     203 * Divide two double-precision floats.
     204 *
     205 * @param a Nominator.
     206 * @param b Denominator.
     207 * @return Result of division.
     208 */
    196209float64 divFloat64(float64 a, float64 b)
    197210{
    … …  
    200213        uint64_t afrac, bfrac, cfrac;
    201214        uint64_t remlo, remhi;
     215        uint64_t tmplo, tmphi;
    202216       
    203217        result.parts.sign = a.parts.sign ^ b.parts.sign;
    204218       
    205219        if (isFloat64NaN(a)) {
    206                
    207220                if (isFloat64SigNaN(b)) {
    208221                        /*FIXME: SigNaN*/
    … …  
    262275        }
    263276
    264        
    265277        afrac = a.parts.fraction;
    266278        aexp = a.parts.exp;
    … …  
    275287                        return result;
    276288                }
     289
    277290                /* normalize it*/
    278                
    279291                aexp++;
    280                         /* afrac is nonzero => it must stop */ 
    281                 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
     292                /* afrac is nonzero => it must stop */ 
     293                while (!(afrac & FLOAT64_HIDDEN_BIT_MASK)) {
    282294                        afrac <<= 1;
    283295                        aexp--;
    … …  
    287299        if (bexp == 0) {
    288300                bexp++;
    289                         /* bfrac is nonzero => it must stop */ 
    290                 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
     301                /* bfrac is nonzero => it must stop */ 
     302                while (!(bfrac & FLOAT64_HIDDEN_BIT_MASK)) {
    291303                        bfrac <<= 1;
    292304                        bexp--;
    … …  
    294306        }
    295307
    296         afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
    297         bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
    298 
    299         if ( bfrac <= (afrac << 1) ) {
     308        afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 2);
     309        bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 1);
     310
     311        if (bfrac <= (afrac << 1)) {
    300312                afrac >>= 1;
    301313                aexp++;
    … …  
    304316        cexp = aexp - bexp + FLOAT64_BIAS - 2;
    305317       
    306         cfrac = divFloat64estim(afrac, bfrac);
    307        
    308         if ((  cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
    309                 mul64integers( bfrac, cfrac, &remlo, &remhi);
    310                 /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/     
    311                 remhi = afrac - remhi - ( remlo > 0);
    312                 remlo = - remlo;
     318        cfrac = div128est(afrac, 0x0ll, bfrac);
     319       
     320        if ((cfrac & 0x1FF) <= 2) {
     321                mul64(bfrac, cfrac, &tmphi, &tmplo);
     322                sub128(afrac, 0x0ll, tmphi, tmplo, &remhi, &remlo);
    313323               
    314324                while ((int64_t) remhi < 0) {
    315325                        cfrac--;
    316                         remlo += bfrac;
    317                         remhi += ( remlo < bfrac );
    318                 }
    319                 cfrac |= ( remlo != 0 );
     326                        add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo);
     327                }
     328                cfrac |= (remlo != 0);
    320329        }
    321330       
    … …  
    323332        result = finishFloat64(cexp, cfrac, result.parts.sign);
    324333        return result;
    325 
    326334}
    327335
    328 uint64_t divFloat64estim(uint64_t a, uint64_t b)
     336/**
     337 * Divide two quadruple-precision floats.
     338 *
     339 * @param a Nominator.
     340 * @param b Denominator.
     341 * @return Result of division.
     342 */
     343float128 divFloat128(float128 a, float128 b)
    329344{
    330         uint64_t bhi;
    331         uint64_t remhi, remlo;
    332         uint64_t result;
    333        
    334         if ( b <= a ) {
    335                 return 0xFFFFFFFFFFFFFFFFull;
    336         }
    337        
    338         bhi = b >> 32;
    339         result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
    340         mul64integers(b, result, &remlo, &remhi);
    341        
    342         remhi = a - remhi - (remlo > 0);
    343         remlo = - remlo;
    344 
    345         b <<= 32;
    346         while ( (int64_t) remhi < 0 ) {
    347                         result -= 0x1ll << 32; 
    348                         remlo += b;
    349                         remhi += bhi + ( remlo < b );
    350                 }
    351         remhi = (remhi << 32) | (remlo >> 32);
    352         if (( bhi << 32) <= remhi) {
    353                 result |= 0xFFFFFFFF;
    354         } else {
    355                 result |= remhi / bhi;
    356         }
    357        
    358        
     345        float128 result;
     346        int64_t aexp, bexp, cexp;
     347        uint64_t afrac_hi, afrac_lo, bfrac_hi, bfrac_lo, cfrac_hi, cfrac_lo;
     348        uint64_t shift_out;
     349        uint64_t rem_hihi, rem_hilo, rem_lohi, rem_lolo;
     350        uint64_t tmp_hihi, tmp_hilo, tmp_lohi, tmp_lolo;
     351
     352        result.parts.sign = a.parts.sign ^ b.parts.sign;
     353
     354        if (isFloat128NaN(a)) {
     355                if (isFloat128SigNaN(b)) {
     356                        /*FIXME: SigNaN*/
     357                        return b;
     358                }
     359
     360                if (isFloat128SigNaN(a)) {
     361                        /*FIXME: SigNaN*/
     362                }
     363                /*NaN*/
     364                return a;
     365        }
     366
     367        if (isFloat128NaN(b)) {
     368                if (isFloat128SigNaN(b)) {
     369                        /*FIXME: SigNaN*/
     370                }
     371                /*NaN*/
     372                return b;
     373        }
     374
     375        if (isFloat128Infinity(a)) {
     376                if (isFloat128Infinity(b) || isFloat128Zero(b)) {
     377                        /*FIXME: inf / inf */
     378                        result.binary.hi = FLOAT128_NAN_HI;
     379                        result.binary.lo = FLOAT128_NAN_LO;
     380                        return result;
     381                }
     382                /* inf / num */
     383                result.parts.exp = a.parts.exp;
     384                result.parts.frac_hi = a.parts.frac_hi;
     385                result.parts.frac_lo = a.parts.frac_lo;
     386                return result;
     387        }
     388
     389        if (isFloat128Infinity(b)) {
     390                if (isFloat128Zero(a)) {
     391                        /* FIXME 0 / inf */
     392                        result.parts.exp = 0;
     393                        result.parts.frac_hi = 0;
     394                        result.parts.frac_lo = 0;
     395                        return result;
     396                }
     397                /* FIXME: num / inf*/
     398                result.parts.exp = 0;
     399                result.parts.frac_hi = 0;
     400                result.parts.frac_lo = 0;
     401                return result;
     402        }
     403
     404        if (isFloat128Zero(b)) {
     405                if (isFloat128Zero(a)) {
     406                        /*FIXME: 0 / 0*/
     407                        result.binary.hi = FLOAT128_NAN_HI;
     408                        result.binary.lo = FLOAT128_NAN_LO;
     409                        return result;
     410                }
     411                /* FIXME: division by zero */
     412                result.parts.exp = 0;
     413                result.parts.frac_hi = 0;
     414                result.parts.frac_lo = 0;
     415                return result;
     416        }
     417
     418        afrac_hi = a.parts.frac_hi;
     419        afrac_lo = a.parts.frac_lo;
     420        aexp = a.parts.exp;
     421        bfrac_hi = b.parts.frac_hi;
     422        bfrac_lo = b.parts.frac_lo;
     423        bexp = b.parts.exp;
     424
     425        /* denormalized numbers */
     426        if (aexp == 0) {
     427                if (eq128(afrac_hi, afrac_lo, 0x0ll, 0x0ll)) {
     428                        result.parts.exp = 0;
     429                        result.parts.frac_hi = 0;
     430                        result.parts.frac_lo = 0;
     431                        return result;
     432                }
     433
     434                /* normalize it*/
     435                aexp++;
     436                /* afrac is nonzero => it must stop */
     437                and128(afrac_hi, afrac_lo,
     438                    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     439                    &tmp_hihi, &tmp_lolo);
     440                while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
     441                        lshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
     442                        aexp--;
     443                }
     444        }
     445
     446        if (bexp == 0) {
     447                bexp++;
     448                /* bfrac is nonzero => it must stop */
     449                and128(bfrac_hi, bfrac_lo,
     450                    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     451                    &tmp_hihi, &tmp_lolo);
     452                while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
     453                        lshift128(bfrac_hi, bfrac_lo, 1, &bfrac_hi, &bfrac_lo);
     454                        bexp--;
     455                }
     456        }
     457
     458        or128(afrac_hi, afrac_lo,
     459            FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     460            &afrac_hi, &afrac_lo);
     461        lshift128(afrac_hi, afrac_lo,
     462            (128 - FLOAT128_FRACTION_SIZE - 1), &afrac_hi, &afrac_lo);
     463        or128(bfrac_hi, bfrac_lo,
     464            FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     465            &bfrac_hi, &bfrac_lo);
     466        lshift128(bfrac_hi, bfrac_lo,
     467            (128 - FLOAT128_FRACTION_SIZE - 1), &bfrac_hi, &bfrac_lo);
     468
     469        if (le128(bfrac_hi, bfrac_lo, afrac_hi, afrac_lo)) {
     470                rshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
     471                aexp++;
     472        }
     473
     474        cexp = aexp - bexp + FLOAT128_BIAS - 2;
     475
     476        cfrac_hi = div128est(afrac_hi, afrac_lo, bfrac_hi);
     477
     478        mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_hi,
     479            &tmp_lolo /* dummy */, &tmp_hihi, &tmp_hilo, &tmp_lohi);
     480
     481        /* sub192(afrac_hi, afrac_lo, 0,
     482         *     tmp_hihi, tmp_hilo, tmp_lohi
     483         *     &rem_hihi, &rem_hilo, &rem_lohi); */
     484        sub128(afrac_hi, afrac_lo, tmp_hihi, tmp_hilo, &rem_hihi, &rem_hilo);
     485        if (tmp_lohi > 0) {
     486                sub128(rem_hihi, rem_hilo, 0x0ll, 0x1ll, &rem_hihi, &rem_hilo);
     487        }
     488        rem_lohi = -tmp_lohi;
     489
     490        while ((int64_t) rem_hihi < 0) {
     491                --cfrac_hi;
     492                /* add192(rem_hihi, rem_hilo, rem_lohi,
     493                 *     0, bfrac_hi, bfrac_lo,
     494                 *     &rem_hihi, &rem_hilo, &rem_lohi); */
     495                add128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo, &rem_hilo, &rem_lohi);
     496                if (lt128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo)) {
     497                        ++rem_hihi;
     498                }
     499        }
     500
     501        cfrac_lo = div128est(rem_hilo, rem_lohi, bfrac_lo);
     502
     503        if ((cfrac_lo & 0x3FFF) <= 4) {
     504                mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_lo,
     505                &tmp_hihi /* dummy */, &tmp_hilo, &tmp_lohi, &tmp_lolo);
     506
     507                /* sub192(rem_hilo, rem_lohi, 0,
     508                 *     tmp_hilo, tmp_lohi, tmp_lolo,
     509                 *     &rem_hilo, &rem_lohi, &rem_lolo); */
     510                sub128(rem_hilo, rem_lohi, tmp_hilo, tmp_lohi, &rem_hilo, &rem_lohi);
     511                if (tmp_lolo > 0) {
     512                        sub128(rem_hilo, rem_lohi, 0x0ll, 0x1ll, &rem_hilo, &rem_lohi);
     513                }
     514                rem_lolo = -tmp_lolo;
     515
     516                while ((int64_t) rem_hilo < 0) {
     517                        --cfrac_lo;
     518                        /* add192(rem_hilo, rem_lohi, rem_lolo,
     519                         *     0, bfrac_hi, bfrac_lo,
     520                         *     &rem_hilo, &rem_lohi, &rem_lolo); */
     521                        add128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo, &rem_lohi, &rem_lolo);
     522                        if (lt128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo)) {
     523                                ++rem_hilo;
     524                        }
     525                }
     526
     527                cfrac_lo |= ((rem_hilo | rem_lohi | rem_lolo) != 0 );
     528        }
     529
     530        shift_out = cfrac_lo << (64 - (128 - FLOAT128_FRACTION_SIZE - 1));
     531        rshift128(cfrac_hi, cfrac_lo, (128 - FLOAT128_FRACTION_SIZE - 1),
     532            &cfrac_hi, &cfrac_lo);
     533
     534        result = finishFloat128(cexp, cfrac_hi, cfrac_lo, result.parts.sign, shift_out);
    359535        return result;
    360536}
Note: See TracChangeset for help on using the changeset viewer.