Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • uspace/lib/softfloat/generic/sub.c

    rc67aff2 r750636a  
    11/*
    22 * Copyright (c) 2005 Josef Cejka
    3  * Copyright (c) 2011 Petr Koupy
    43 * All rights reserved.
    54 *
     
    3130 * @{
    3231 */
    33 /** @file Substraction functions.
     32/** @file
    3433 */
    3534
     
    3736#include <sub.h>
    3837#include <comparison.h>
    39 #include <common.h>
    40 
    41 /**
    42  * Subtract two single-precision floats with the same signs.
    43  *
    44  * @param a First input operand.
    45  * @param b Second input operand.
    46  * @return Result of substraction.
     38
     39/** Subtract two float32 numbers with same signs
    4740 */
    4841float32 subFloat32(float32 a, float32 b)
     
    5952                        /* TODO: fix SigNaN */
    6053                        if (isFloat32SigNaN(b)) {
    61                         }
    62                         return b;
    63                 }
     54                        };
     55                        return b;
     56                };
    6457               
    6558                if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
     
    7972                        /* TODO: fix SigNaN */
    8073                        if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
    81                         }
    82                         return a;
    83                 }
     74                        };
     75                        return a;
     76                };
    8477               
    8578                if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
     
    8982                                result.binary = FLOAT32_NAN;
    9083                                return result;
    91                         }
     84                        };
    9285                        return a;
    9386                }
     
    9992                frac2 = b.parts.fraction;
    10093                exp2 = b.parts.exp;     
    101         }
     94        };
     95       
     96        if (exp1 == 0) {
     97                /* both are denormalized */
     98                result.parts.fraction = frac1-frac2;
     99                if (result.parts.fraction > frac1) {
     100                        /* TODO: underflow exception */
     101                        return result;
     102                };
     103                result.parts.exp = 0;
     104                return result;
     105        };
     106
     107        /* add hidden bit */
     108        frac1 |= FLOAT32_HIDDEN_BIT_MASK;
     109       
     110        if (exp2 == 0) {
     111                /* denormalized */
     112                --expdiff;     
     113        } else {
     114                /* normalized */
     115                frac2 |= FLOAT32_HIDDEN_BIT_MASK;
     116        };
     117       
     118        /* create some space for rounding */
     119        frac1 <<= 6;
     120        frac2 <<= 6;
     121       
     122        if (expdiff > FLOAT32_FRACTION_SIZE + 1) {
     123             goto done;
     124             };
     125       
     126        frac1 = frac1 - (frac2 >> expdiff);
     127done:
     128        /* TODO: find first nonzero digit and shift result and detect possibly underflow */
     129        while ((exp1 > 0) && (!(frac1 & (FLOAT32_HIDDEN_BIT_MASK << 6 )))) {
     130                --exp1;
     131                frac1 <<= 1;
     132                        /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
     133        };
     134       
     135        /* rounding - if first bit after fraction is set then round up */
     136        frac1 += 0x20;
     137
     138        if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
     139                ++exp1;
     140                frac1 >>= 1;
     141        };
     142       
     143        /*Clear hidden bit and shift */
     144        result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
     145        result.parts.exp = exp1;
     146       
     147        return result;
     148}
     149
     150/** Subtract two float64 numbers with same signs
     151 */
     152float64 subFloat64(float64 a, float64 b)
     153{
     154        int expdiff;
     155        uint32_t exp1, exp2;
     156        uint64_t frac1, frac2;
     157        float64 result;
     158
     159        result.d = 0;
     160       
     161        expdiff = a.parts.exp - b.parts.exp;
     162        if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.fraction < b.parts.fraction))) {
     163                if (isFloat64NaN(b)) {
     164                        /* TODO: fix SigNaN */
     165                        if (isFloat64SigNaN(b)) {
     166                        };
     167                        return b;
     168                };
     169               
     170                if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
     171                        b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
     172                        return b;
     173                }
     174               
     175                result.parts.sign = !a.parts.sign;
     176               
     177                frac1 = b.parts.fraction;
     178                exp1 = b.parts.exp;
     179                frac2 = a.parts.fraction;
     180                exp2 = a.parts.exp;
     181                expdiff *= -1;
     182        } else {
     183                if (isFloat64NaN(a)) {
     184                        /* TODO: fix SigNaN */
     185                        if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
     186                        };
     187                        return a;
     188                };
     189               
     190                if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
     191                        if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
     192                                /* inf - inf => nan */
     193                                /* TODO: fix exception */
     194                                result.binary = FLOAT64_NAN;
     195                                return result;
     196                        };
     197                        return a;
     198                }
     199               
     200                result.parts.sign = a.parts.sign;
     201               
     202                frac1 = a.parts.fraction;
     203                exp1 = a.parts.exp;
     204                frac2 = b.parts.fraction;
     205                exp2 = b.parts.exp;     
     206        };
    102207       
    103208        if (exp1 == 0) {
     
    107212                        /* TODO: underflow exception */
    108213                        return result;
    109                 }
     214                };
    110215                result.parts.exp = 0;
    111216                return result;
    112         }
     217        };
    113218
    114219        /* add hidden bit */
    115         frac1 |= FLOAT32_HIDDEN_BIT_MASK;
     220        frac1 |= FLOAT64_HIDDEN_BIT_MASK;
    116221       
    117222        if (exp2 == 0) {
     
    120225        } else {
    121226                /* normalized */
    122                 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
    123         }
     227                frac2 |= FLOAT64_HIDDEN_BIT_MASK;
     228        };
    124229       
    125230        /* create some space for rounding */
     
    127232        frac2 <<= 6;
    128233       
    129         if (expdiff > FLOAT32_FRACTION_SIZE + 1) {
    130                 goto done;
    131         }
     234        if (expdiff > FLOAT64_FRACTION_SIZE + 1) {
     235             goto done;
     236             };
    132237       
    133238        frac1 = frac1 - (frac2 >> expdiff);
    134 
    135 done:
    136         /* TODO: find first nonzero digit and shift result and detect possibly underflow */
    137         while ((exp1 > 0) && (!(frac1 & (FLOAT32_HIDDEN_BIT_MASK << 6 )))) {
    138                 --exp1;
    139                 frac1 <<= 1;
    140                 /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
    141         }
    142        
    143         /* rounding - if first bit after fraction is set then round up */
    144         frac1 += 0x20;
    145 
    146         if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
    147                 ++exp1;
    148                 frac1 >>= 1;
    149         }
    150        
    151         /* Clear hidden bit and shift */
    152         result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
    153         result.parts.exp = exp1;
    154        
    155         return result;
    156 }
    157 
    158 /**
    159  * Subtract two double-precision floats with the same signs.
    160  *
    161  * @param a First input operand.
    162  * @param b Second input operand.
    163  * @return Result of substraction.
    164  */
    165 float64 subFloat64(float64 a, float64 b)
    166 {
    167         int expdiff;
    168         uint32_t exp1, exp2;
    169         uint64_t frac1, frac2;
    170         float64 result;
    171 
    172         result.d = 0;
    173        
    174         expdiff = a.parts.exp - b.parts.exp;
    175         if ((expdiff < 0 ) || ((expdiff == 0) && (a.parts.fraction < b.parts.fraction))) {
    176                 if (isFloat64NaN(b)) {
    177                         /* TODO: fix SigNaN */
    178                         if (isFloat64SigNaN(b)) {
    179                         }
    180                         return b;
    181                 }
    182                
    183                 if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
    184                         b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
    185                         return b;
    186                 }
    187                
    188                 result.parts.sign = !a.parts.sign;
    189                
    190                 frac1 = b.parts.fraction;
    191                 exp1 = b.parts.exp;
    192                 frac2 = a.parts.fraction;
    193                 exp2 = a.parts.exp;
    194                 expdiff *= -1;
    195         } else {
    196                 if (isFloat64NaN(a)) {
    197                         /* TODO: fix SigNaN */
    198                         if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
    199                         }
    200                         return a;
    201                 }
    202                
    203                 if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
    204                         if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
    205                                 /* inf - inf => nan */
    206                                 /* TODO: fix exception */
    207                                 result.binary = FLOAT64_NAN;
    208                                 return result;
    209                         }
    210                         return a;
    211                 }
    212                
    213                 result.parts.sign = a.parts.sign;
    214                
    215                 frac1 = a.parts.fraction;
    216                 exp1 = a.parts.exp;
    217                 frac2 = b.parts.fraction;
    218                 exp2 = b.parts.exp;     
    219         }
    220        
    221         if (exp1 == 0) {
    222                 /* both are denormalized */
    223                 result.parts.fraction = frac1 - frac2;
    224                 if (result.parts.fraction > frac1) {
    225                         /* TODO: underflow exception */
    226                         return result;
    227                 }
    228                 result.parts.exp = 0;
    229                 return result;
    230         }
    231 
    232         /* add hidden bit */
    233         frac1 |= FLOAT64_HIDDEN_BIT_MASK;
    234        
    235         if (exp2 == 0) {
    236                 /* denormalized */
    237                 --expdiff;     
    238         } else {
    239                 /* normalized */
    240                 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
    241         }
    242        
    243         /* create some space for rounding */
    244         frac1 <<= 6;
    245         frac2 <<= 6;
    246        
    247         if (expdiff > FLOAT64_FRACTION_SIZE + 1) {
    248                 goto done;
    249         }
    250        
    251         frac1 = frac1 - (frac2 >> expdiff);
    252 
    253239done:
    254240        /* TODO: find first nonzero digit and shift result and detect possibly underflow */
     
    256242                --exp1;
    257243                frac1 <<= 1;
    258                 /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
    259         }
     244                        /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
     245        };
    260246       
    261247        /* rounding - if first bit after fraction is set then round up */
     
    265251                ++exp1;
    266252                frac1 >>= 1;
    267         }
    268        
    269         /* Clear hidden bit and shift */
     253        };
     254       
     255        /*Clear hidden bit and shift */
    270256        result.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK));
    271257        result.parts.exp = exp1;
     
    274260}
    275261
    276 /**
    277  * Subtract two quadruple-precision floats with the same signs.
    278  *
    279  * @param a First input operand.
    280  * @param b Second input operand.
    281  * @return Result of substraction.
    282  */
    283 float128 subFloat128(float128 a, float128 b)
    284 {
    285         int expdiff;
    286         uint32_t exp1, exp2;
    287         uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
    288         float128 result;
    289 
    290         result.binary.hi = 0;
    291         result.binary.lo = 0;
    292 
    293         expdiff = a.parts.exp - b.parts.exp;
    294         if ((expdiff < 0 ) || ((expdiff == 0) &&
    295             lt128(a.parts.frac_hi, a.parts.frac_lo, b.parts.frac_hi, b.parts.frac_lo))) {
    296                 if (isFloat128NaN(b)) {
    297                         /* TODO: fix SigNaN */
    298                         if (isFloat128SigNaN(b)) {
    299                         }
    300                         return b;
    301                 }
    302 
    303                 if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
    304                         b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
    305                         return b;
    306                 }
    307 
    308                 result.parts.sign = !a.parts.sign;
    309 
    310                 frac1_hi = b.parts.frac_hi;
    311                 frac1_lo = b.parts.frac_lo;
    312                 exp1 = b.parts.exp;
    313                 frac2_hi = a.parts.frac_hi;
    314                 frac2_lo = a.parts.frac_lo;
    315                 exp2 = a.parts.exp;
    316                 expdiff *= -1;
    317         } else {
    318                 if (isFloat128NaN(a)) {
    319                         /* TODO: fix SigNaN */
    320                         if (isFloat128SigNaN(a) || isFloat128SigNaN(b)) {
    321                         }
    322                         return a;
    323                 }
    324 
    325                 if (a.parts.exp == FLOAT128_MAX_EXPONENT) {
    326                         if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
    327                                 /* inf - inf => nan */
    328                                 /* TODO: fix exception */
    329                                 result.binary.hi = FLOAT128_NAN_HI;
    330                                 result.binary.lo = FLOAT128_NAN_LO;
    331                                 return result;
    332                         }
    333                         return a;
    334                 }
    335 
    336                 result.parts.sign = a.parts.sign;
    337 
    338                 frac1_hi = a.parts.frac_hi;
    339                 frac1_lo = a.parts.frac_lo;
    340                 exp1 = a.parts.exp;
    341                 frac2_hi = b.parts.frac_hi;
    342                 frac2_lo = b.parts.frac_lo;
    343                 exp2 = b.parts.exp;
    344         }
    345 
    346         if (exp1 == 0) {
    347                 /* both are denormalized */
    348                 sub128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &tmp_hi, &tmp_lo);
    349                 result.parts.frac_hi = tmp_hi;
    350                 result.parts.frac_lo = tmp_lo;
    351                 if (lt128(frac1_hi, frac1_lo, result.parts.frac_hi, result.parts.frac_lo)) {
    352                         /* TODO: underflow exception */
    353                         return result;
    354                 }
    355                 result.parts.exp = 0;
    356                 return result;
    357         }
    358 
    359         /* add hidden bit */
    360         or128(frac1_hi, frac1_lo,
    361             FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
    362             &frac1_hi, &frac1_lo);
    363 
    364         if (exp2 == 0) {
    365                 /* denormalized */
    366                 --expdiff;
    367         } else {
    368                 /* normalized */
    369                 or128(frac2_hi, frac2_lo,
    370                     FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
    371                     &frac2_hi, &frac2_lo);
    372         }
    373 
    374         /* create some space for rounding */
    375         lshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
    376         lshift128(frac2_hi, frac2_lo, 6, &frac2_hi, &frac2_lo);
    377 
    378         if (expdiff > FLOAT128_FRACTION_SIZE + 1) {
    379                 goto done;
    380         }
    381 
    382         rshift128(frac2_hi, frac2_lo, expdiff, &tmp_hi, &tmp_lo);
    383         sub128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo);
    384 
    385 done:
    386         /* TODO: find first nonzero digit and shift result and detect possibly underflow */
    387         lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6,
    388             &tmp_hi, &tmp_lo);
    389         and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
    390         while ((exp1 > 0) && (!lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo))) {
    391                 --exp1;
    392                 lshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
    393                 /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
    394 
    395                 lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6,
    396                     &tmp_hi, &tmp_lo);
    397                 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
    398         }
    399 
    400         /* rounding - if first bit after fraction is set then round up */
    401         add128(frac1_hi, frac1_lo, 0x0ll, 0x20ll, &frac1_hi, &frac1_lo);
    402 
    403         lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7,
    404            &tmp_hi, &tmp_lo);
    405         and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
    406         if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
    407                 ++exp1;
    408                 rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
    409         }
    410 
    411         /* Clear hidden bit and shift */
    412         rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
    413         not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
    414             &tmp_hi, &tmp_lo);
    415         and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
    416         result.parts.frac_hi = tmp_hi;
    417         result.parts.frac_lo = tmp_lo;
    418 
    419         result.parts.exp = exp1;
    420 
    421         return result;
    422 }
    423 
    424262/** @}
    425263 */
Note: See TracChangeset for help on using the changeset viewer.