Changeset c67aff2 in mainline for uspace/lib/softfloat/generic/add.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/add.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 Addition functions.
    3334 */
    3435
    … …  
    3637#include <add.h>
    3738#include <comparison.h>
    38 
    39 /** Add two Float32 numbers with same signs
     39#include <common.h>
     40
     41/**
     42 * Add 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 addition.
    4047 */
    4148float32 addFloat32(float32 a, float32 b)
    4249{
    4350        int expdiff;
    44         uint32_t exp1, exp2,frac1, frac2;
     51        uint32_t exp1, exp2, frac1, frac2;
    4552       
    4653        expdiff = a.parts.exp - b.parts.exp;
    … …  
    4956                        /* TODO: fix SigNaN */
    5057                        if (isFloat32SigNaN(b)) {
    51                         };
    52 
    53                         return b;
    54                 };
     58                        }
     59
     60                        return b;
     61                }
    5562               
    5663                if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
    … …  
    6774                        /* TODO: fix SigNaN */
    6875                        if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
    69                         };
    70                         return (isFloat32NaN(a)?a:b);
    71                 };
     76                        }
     77                        return (isFloat32NaN(a) ? a : b);
     78                }
    7279               
    7380                if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
    … …  
    7986                frac2 = b.parts.fraction;
    8087                exp2 = b.parts.exp;
    81         };
     88        }
    8289       
    8390        if (exp1 == 0) {
    … …  
    8794                        /* result is not denormalized */
    8895                        a.parts.exp = 1;
    89                 };
     96                }
    9097                a.parts.fraction = frac1;
    9198                return a;
    92         };
     99        }
    93100       
    94101        frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
    … …  
    100107                /* add hidden bit to second operand */
    101108                frac2 |= FLOAT32_HIDDEN_BIT_MASK;
    102         };
     109        }
    103110       
    104111        /* create some space for rounding */
    … …  
    118125                ++exp1;
    119126                frac1 >>= 1;
    120         };
     127        }
    121128       
    122129        /* rounding - if first bit after fraction is set then round up */
    … …  
    127134                ++exp1;
    128135                frac1 >>= 1;
    129         };
    130        
     136        }
    131137       
    132138        if ((exp1 == FLOAT32_MAX_EXPONENT ) || (exp2 > exp1)) {
    133                         /* overflow - set infinity as result */
    134                         a.parts.exp = FLOAT32_MAX_EXPONENT;
    135                         a.parts.fraction = 0;
    136                         return a;
    137                         }
     139                /* overflow - set infinity as result */
     140                a.parts.exp = FLOAT32_MAX_EXPONENT;
     141                a.parts.fraction = 0;
     142                return a;
     143        }
    138144       
    139145        a.parts.exp = exp1;
    140146       
    141147        /* Clear hidden bit and shift */
    142         a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)) ;
     148        a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
    143149        return a;
    144150}
    145151
    146 /** Add two Float64 numbers with same signs
     152/**
     153 * Add two double-precision floats with the same signs.
     154 *
     155 * @param a First input operand.
     156 * @param b Second input operand.
     157 * @return Result of addition.
    147158 */
    148159float64 addFloat64(float64 a, float64 b)
    … …  
    152163        uint64_t frac1, frac2;
    153164       
    154         expdiff = ((int )a.parts.exp) - b.parts.exp;
     165        expdiff = ((int) a.parts.exp) - b.parts.exp;
    155166        if (expdiff < 0) {
    156167                if (isFloat64NaN(b)) {
    157168                        /* TODO: fix SigNaN */
    158169                        if (isFloat64SigNaN(b)) {
    159                         };
    160 
    161                         return b;
    162                 };
     170                        }
     171
     172                        return b;
     173                }
    163174               
    164175                /* b is infinity and a not */   
    165                 if (b.parts.exp == FLOAT64_MAX_EXPONENT ) {
     176                if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
    166177                        return b;
    167178                }
    … …  
    176187                        /* TODO: fix SigNaN */
    177188                        if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
    178                         };
     189                        }
    179190                        return a;
    180                 };
     191                }
    181192               
    182193                /* a is infinity and b not */
    183                 if (a.parts.exp == FLOAT64_MAX_EXPONENT ) {
     194                if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
    184195                        return a;
    185196                }
    … …  
    189200                frac2 = b.parts.fraction;
    190201                exp2 = b.parts.exp;
    191         };
     202        }
    192203       
    193204        if (exp1 == 0) {
    … …  
    197208                        /* result is not denormalized */
    198209                        a.parts.exp = 1;
    199                 };
     210                }
    200211                a.parts.fraction = frac1;
    201212                return a;
    202         };
     213        }
    203214       
    204215        /* add hidden bit - frac1 is sure not denormalized */
    … …  
    212223                /* is not denormalized */
    213224                frac2 |= FLOAT64_HIDDEN_BIT_MASK;
    214         };
     225        }
    215226       
    216227        /* create some space for rounding */
    … …  
    218229        frac2 <<= 6;
    219230       
    220         if (expdiff < (FLOAT64_FRACTION_SIZE + 2) ) {
     231        if (expdiff < (FLOAT64_FRACTION_SIZE + 2)) {
    221232                frac2 >>= expdiff;
    222233                frac1 += frac2;
    … …  
    227238        }
    228239       
    229         if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7) ) {
     240        if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
    230241                ++exp1;
    231242                frac1 >>= 1;
    232         };
     243        }
    233244       
    234245        /* rounding - if first bit after fraction is set then round up */
    … …  
    239250                ++exp1;
    240251                frac1 >>= 1;
    241         };
     252        }
    242253       
    243254        if ((exp1 == FLOAT64_MAX_EXPONENT ) || (exp2 > exp1)) {
    244                         /* overflow - set infinity as result */
    245                         a.parts.exp = FLOAT64_MAX_EXPONENT;
    246                         a.parts.fraction = 0;
    247                         return a;
    248                         }
     255                /* overflow - set infinity as result */
     256                a.parts.exp = FLOAT64_MAX_EXPONENT;
     257                a.parts.fraction = 0;
     258                return a;
     259        }
    249260       
    250261        a.parts.exp = exp1;
    251262        /* Clear hidden bit and shift */
    252         a.parts.fraction = ( (frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
    253        
     263        a.parts.fraction = ((frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
    254264        return a;
    255265}
    256266
     267/**
     268 * Add two quadruple-precision floats with the same signs.
     269 *
     270 * @param a First input operand.
     271 * @param b Second input operand.
     272 * @return Result of addition.
     273 */
     274float128 addFloat128(float128 a, float128 b)
     275{
     276        int expdiff;
     277        uint32_t exp1, exp2;
     278        uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
     279
     280        expdiff = ((int) a.parts.exp) - b.parts.exp;
     281        if (expdiff < 0) {
     282                if (isFloat128NaN(b)) {
     283                        /* TODO: fix SigNaN */
     284                        if (isFloat128SigNaN(b)) {
     285                        }
     286
     287                        return b;
     288                }
     289
     290                /* b is infinity and a not */
     291                if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
     292                        return b;
     293                }
     294
     295                frac1_hi = b.parts.frac_hi;
     296                frac1_lo = b.parts.frac_lo;
     297                exp1 = b.parts.exp;
     298                frac2_hi = a.parts.frac_hi;
     299                frac2_lo = a.parts.frac_lo;
     300                exp2 = a.parts.exp;
     301                expdiff *= -1;
     302        } else {
     303                if (isFloat128NaN(a)) {
     304                        /* TODO: fix SigNaN */
     305                        if (isFloat128SigNaN(a) || isFloat128SigNaN(b)) {
     306                        }
     307                        return a;
     308                }
     309
     310                /* a is infinity and b not */
     311                if (a.parts.exp == FLOAT128_MAX_EXPONENT) {
     312                        return a;
     313                }
     314
     315                frac1_hi = a.parts.frac_hi;
     316                frac1_lo = a.parts.frac_lo;
     317                exp1 = a.parts.exp;
     318                frac2_hi = b.parts.frac_hi;
     319                frac2_lo = b.parts.frac_lo;
     320                exp2 = b.parts.exp;
     321        }
     322
     323        if (exp1 == 0) {
     324                /* both are denormalized */
     325                add128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &frac1_hi, &frac1_lo);
     326
     327                and128(frac1_hi, frac1_lo,
     328                    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     329                    &tmp_hi, &tmp_lo);
     330                if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
     331                        /* result is not denormalized */
     332                        a.parts.exp = 1;
     333                }
     334
     335                a.parts.frac_hi = frac1_hi;
     336                a.parts.frac_lo = frac1_lo;
     337                return a;
     338        }
     339
     340        /* add hidden bit - frac1 is sure not denormalized */
     341        or128(frac1_hi, frac1_lo,
     342            FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     343            &frac1_hi, &frac1_lo);
     344
     345        /* second operand ... */
     346        if (exp2 == 0) {
     347                /* ... is denormalized */
     348                --expdiff;
     349        } else {
     350                /* is not denormalized */
     351                or128(frac2_hi, frac2_lo,
     352                    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     353                    &frac2_hi, &frac2_lo);
     354        }
     355
     356        /* create some space for rounding */
     357        lshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
     358        lshift128(frac2_hi, frac2_lo, 6, &frac2_hi, &frac2_lo);
     359
     360        if (expdiff < (FLOAT128_FRACTION_SIZE + 2)) {
     361                rshift128(frac2_hi, frac2_lo, expdiff, &frac2_hi, &frac2_lo);
     362                add128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &frac1_hi, &frac1_lo);
     363        } else {
     364                a.parts.exp = exp1;
     365
     366                rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
     367                not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     368                    &tmp_hi, &tmp_lo);
     369                and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
     370
     371                a.parts.frac_hi = tmp_hi;
     372                a.parts.frac_lo = tmp_lo;
     373                return a;
     374        }
     375
     376        lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7,
     377            &tmp_hi, &tmp_lo);
     378        and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
     379        if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
     380                ++exp1;
     381                rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
     382        }
     383
     384        /* rounding - if first bit after fraction is set then round up */
     385        add128(frac1_hi, frac1_lo, 0x0ll, 0x1ll << 5, &frac1_hi, &frac1_lo);
     386
     387        lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7,
     388           &tmp_hi, &tmp_lo);
     389        and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
     390        if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
     391                /* rounding overflow */
     392                ++exp1;
     393                rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
     394        }
     395
     396        if ((exp1 == FLOAT128_MAX_EXPONENT ) || (exp2 > exp1)) {
     397                /* overflow - set infinity as result */
     398                a.parts.exp = FLOAT64_MAX_EXPONENT;
     399                a.parts.frac_hi = 0;
     400                a.parts.frac_lo = 0;
     401                return a;
     402        }
     403
     404        a.parts.exp = exp1;
     405       
     406        /* Clear hidden bit and shift */
     407        rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
     408        not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     409            &tmp_hi, &tmp_lo);
     410        and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
     411
     412        a.parts.frac_hi = tmp_hi;
     413        a.parts.frac_lo = tmp_lo;
     414
     415        return a;
     416}
     417
    257418/** @}
    258419 */
Note: See TracChangeset for help on using the changeset viewer.