Changeset c67aff2 in mainline for uspace/lib/softfloat/generic/sub.c


Ignore:
Timestamp:
2011-08-06T07:04:50Z (13 years ago)
Author:
Petr Koupy <petr.koupy@…>
Branches:
lfn, master, serial, ticket/834-toolchain-update, topic/msim-upgrade, topic/simplify-dev-export
Children:
d3e241a, e0e922d
Parents:
9a6034a
Message:

Quadruple-precision softfloat, coding style improvements. Details below…

Highlights:

  • completed double-precision support
  • added quadruple-precision support
  • added SPARC quadruple-precision wrappers
  • added doxygen comments
  • corrected and unified coding style

Current state of the softfloat library:

Support for single, double and quadruple precision is currently almost complete (apart from power, square root, complex multiplication and complex division) and provides the same set of features (i.e. the support for all three precisions is now aligned). In order to extend softfloat library consistently, addition of quadruple precision was done in the same spirit as already existing single and double precision written by Josef Cejka in 2006 - that is relaxed standard-compliance for corner cases while mission-critical code sections heavily inspired by the widely used softfloat library written by John R. Hauser (although I personally think it would be more appropriate for HelenOS to use something less optimized, shorter and more readable).

Most of the quadruple-precision code is just an adapted double-precision code to work on 128-bit variables. That means if there is TODO, FIXME or some defect in single or double-precision code, it is most likely also in the quadruple-precision code. Please note that quadruple-precision functions are currently not tested - it is challenging task for itself, especially when the ports that use them are either not finished (mips64) or badly supported by simulators (sparc64). To test whole softfloat library, one would probably have to either write very non-trivial native tester, or use some existing one (e.g. TestFloat from J. R. Hauser) and port it to HelenOS (or rip the softfloat library out of HelenOS and test it on a host system). At the time of writing this, the code dependent on quadruple-precision functions (on mips64 and sparc64) is just a libposix strtold() function (and its callers, most notably scanf backend).

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9a6034a 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.