Ignore:
File:
1 edited

Legend:

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

    r750636a rc67aff2  
    11/*
    22 * Copyright (c) 2005 Josef Cejka
     3 * Copyright (c) 2011 Petr Koupy
    34 * All rights reserved.
    45 *
     
    3031 * @{
    3132 */
    32 /** @file
     33/** @file Substraction functions.
    3334 */
    3435
     
    3637#include <sub.h>
    3738#include <comparison.h>
    38 
    39 /** Subtract two float32 numbers with same signs
     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.
    4047 */
    4148float32 subFloat32(float32 a, float32 b)
     
    5259                        /* TODO: fix SigNaN */
    5360                        if (isFloat32SigNaN(b)) {
    54                         };
    55                         return b;
    56                 };
     61                        }
     62                        return b;
     63                }
    5764               
    5865                if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
     
    7279                        /* TODO: fix SigNaN */
    7380                        if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
    74                         };
    75                         return a;
    76                 };
     81                        }
     82                        return a;
     83                }
    7784               
    7885                if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
     
    8289                                result.binary = FLOAT32_NAN;
    8390                                return result;
    84                         };
     91                        }
    8592                        return a;
    8693                }
     
    9299                frac2 = b.parts.fraction;
    93100                exp2 = b.parts.exp;     
    94         };
     101        }
    95102       
    96103        if (exp1 == 0) {
    97104                /* both are denormalized */
    98                 result.parts.fraction = frac1-frac2;
     105                result.parts.fraction = frac1 - frac2;
    99106                if (result.parts.fraction > frac1) {
    100107                        /* TODO: underflow exception */
    101108                        return result;
    102                 };
     109                }
    103110                result.parts.exp = 0;
    104111                return result;
    105         };
     112        }
    106113
    107114        /* add hidden bit */
     
    114121                /* normalized */
    115122                frac2 |= FLOAT32_HIDDEN_BIT_MASK;
    116         };
     123        }
    117124       
    118125        /* create some space for rounding */
     
    121128       
    122129        if (expdiff > FLOAT32_FRACTION_SIZE + 1) {
    123              goto done;
    124              };
     130                goto done;
     131        }
    125132       
    126133        frac1 = frac1 - (frac2 >> expdiff);
     134
    127135done:
    128136        /* TODO: find first nonzero digit and shift result and detect possibly underflow */
     
    130138                --exp1;
    131139                frac1 <<= 1;
    132                         /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
    133         };
     140                /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
     141        }
    134142       
    135143        /* rounding - if first bit after fraction is set then round up */
     
    139147                ++exp1;
    140148                frac1 >>= 1;
    141         };
    142        
    143         /*Clear hidden bit and shift */
     149        }
     150       
     151        /* Clear hidden bit and shift */
    144152        result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
    145153        result.parts.exp = exp1;
     
    148156}
    149157
    150 /** Subtract two float64 numbers with same signs
     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.
    151164 */
    152165float64 subFloat64(float64 a, float64 b)
     
    164177                        /* TODO: fix SigNaN */
    165178                        if (isFloat64SigNaN(b)) {
    166                         };
    167                         return b;
    168                 };
     179                        }
     180                        return b;
     181                }
    169182               
    170183                if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
     
    184197                        /* TODO: fix SigNaN */
    185198                        if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
    186                         };
    187                         return a;
    188                 };
     199                        }
     200                        return a;
     201                }
    189202               
    190203                if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
     
    194207                                result.binary = FLOAT64_NAN;
    195208                                return result;
    196                         };
     209                        }
    197210                        return a;
    198211                }
     
    204217                frac2 = b.parts.fraction;
    205218                exp2 = b.parts.exp;     
    206         };
     219        }
    207220       
    208221        if (exp1 == 0) {
     
    212225                        /* TODO: underflow exception */
    213226                        return result;
    214                 };
     227                }
    215228                result.parts.exp = 0;
    216229                return result;
    217         };
     230        }
    218231
    219232        /* add hidden bit */
     
    226239                /* normalized */
    227240                frac2 |= FLOAT64_HIDDEN_BIT_MASK;
    228         };
     241        }
    229242       
    230243        /* create some space for rounding */
     
    233246       
    234247        if (expdiff > FLOAT64_FRACTION_SIZE + 1) {
    235              goto done;
    236              };
     248                goto done;
     249        }
    237250       
    238251        frac1 = frac1 - (frac2 >> expdiff);
     252
    239253done:
    240254        /* TODO: find first nonzero digit and shift result and detect possibly underflow */
     
    242256                --exp1;
    243257                frac1 <<= 1;
    244                         /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
    245         };
     258                /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
     259        }
    246260       
    247261        /* rounding - if first bit after fraction is set then round up */
     
    251265                ++exp1;
    252266                frac1 >>= 1;
    253         };
    254        
    255         /*Clear hidden bit and shift */
     267        }
     268       
     269        /* Clear hidden bit and shift */
    256270        result.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK));
    257271        result.parts.exp = exp1;
     
    260274}
    261275
     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 */
     283float128 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
     385done:
     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
    262424/** @}
    263425 */
Note: See TracChangeset for help on using the changeset viewer.