Changeset c67aff2 in mainline for uspace/lib/softfloat/generic/mul.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/mul.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 Multiplication functions.
    3334 */
    3435
    … …  
    3839#include <common.h>
    3940
    40 /** Multiply two 32 bit float numbers
    41  *
     41/**
     42 * Multiply two single-precision floats.
     43 *
     44 * @param a First input operand.
     45 * @param b Second input operand.
     46 * @return Result of multiplication.
    4247 */
    4348float32 mulFloat32(float32 a, float32 b)
    … …  
    4954        result.parts.sign = a.parts.sign ^ b.parts.sign;
    5055       
    51         if (isFloat32NaN(a) || isFloat32NaN(b) ) {
     56        if (isFloat32NaN(a) || isFloat32NaN(b)) {
    5257                /* TODO: fix SigNaNs */
    5358                if (isFloat32SigNaN(a)) {
    … …  
    5560                        result.parts.exp = a.parts.exp;
    5661                        return result;
    57                 };
     62                }
    5863                if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
    5964                        result.parts.fraction = b.parts.fraction;
    6065                        result.parts.exp = b.parts.exp;
    6166                        return result;
    62                 };
     67                }
    6368                /* set NaN as result */
    6469                result.binary = FLOAT32_NAN;
    6570                return result;
    66         };
     71        }
    6772               
    6873        if (isFloat32Infinity(a)) {
    … …  
    98103                result.parts.sign = a.parts.sign ^ b.parts.sign;
    99104                return result;
    100         };
     105        }
    101106       
    102107        if (exp < 0) {
    … …  
    106111                result.parts.exp = 0x0;
    107112                return result;
    108         };
     113        }
    109114       
    110115        frac1 = a.parts.fraction;
    … …  
    113118        } else {
    114119                ++exp;
    115         };
     120        }
    116121       
    117122        frac2 = b.parts.fraction;
    … …  
    121126        } else {
    122127                ++exp;
    123         };
     128        }
    124129
    125130        frac1 <<= 1; /* one bit space for rounding */
    126131
    127132        frac1 = frac1 * frac2;
    128 /* round and return */
    129        
    130         while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) {
    131                 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
     133
     134        /* round and return */
     135        while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 2)))) {
     136                /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left) */
    132137                ++exp;
    133138                frac1 >>= 1;
    134         };
     139        }
    135140
    136141        /* rounding */
    … …  
    141146                ++exp;
    142147                frac1 >>= 1;
    143         };
    144 
    145         if (exp >= FLOAT32_MAX_EXPONENT ) {     
     148        }
     149
     150        if (exp >= FLOAT32_MAX_EXPONENT) {     
    146151                /* TODO: fix overflow */
    147152                /* return infinity*/
    … …  
    159164                        frac1 >>= 1;
    160165                        ++exp;
    161                 };
     166                }
    162167                if (frac1 == 0) {
    163168                        /* FIXME : underflow */
    164                 result.parts.exp = 0;
    165                 result.parts.fraction = 0;
    166                 return result;
    167                 };
    168         };
     169                        result.parts.exp = 0;
     170                        result.parts.fraction = 0;
     171                        return result;
     172                }
     173        }
    169174        result.parts.exp = exp;
    170         result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
     175        result.parts.fraction = frac1 & ((1 << FLOAT32_FRACTION_SIZE) - 1);
    171176       
    172177        return result; 
    173        
    174178}
    175179
    176 /** Multiply two 64 bit float numbers
    177  *
     180/**
     181 * Multiply two double-precision floats.
     182 *
     183 * @param a First input operand.
     184 * @param b Second input operand.
     185 * @return Result of multiplication.
    178186 */
    179187float64 mulFloat64(float64 a, float64 b)
    … …  
    185193        result.parts.sign = a.parts.sign ^ b.parts.sign;
    186194       
    187         if (isFloat64NaN(a) || isFloat64NaN(b) ) {
     195        if (isFloat64NaN(a) || isFloat64NaN(b)) {
    188196                /* TODO: fix SigNaNs */
    189197                if (isFloat64SigNaN(a)) {
    … …  
    191199                        result.parts.exp = a.parts.exp;
    192200                        return result;
    193                 };
     201                }
    194202                if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
    195203                        result.parts.fraction = b.parts.fraction;
    196204                        result.parts.exp = b.parts.exp;
    197205                        return result;
    198                 };
     206                }
    199207                /* set NaN as result */
    200208                result.binary = FLOAT64_NAN;
    201209                return result;
    202         };
     210        }
    203211               
    204212        if (isFloat64Infinity(a)) {
    … …  
    233241        } else {
    234242                ++exp;
    235         };
     243        }
    236244       
    237245        frac2 = b.parts.fraction;
    … …  
    241249        } else {
    242250                ++exp;
    243         };
     251        }
    244252
    245253        frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
    246254        frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
    247255
    248         mul64integers(frac1, frac2, &frac1, &frac2);
    249 
    250         frac2 |= (frac1 != 0);
    251         if (frac2 & (0x1ll << 62)) {
    252                 frac2 <<= 1;
     256        mul64(frac1, frac2, &frac1, &frac2);
     257
     258        frac1 |= (frac2 != 0);
     259        if (frac1 & (0x1ll << 62)) {
     260                frac1 <<= 1;
    253261                exp--;
    254262        }
    255263
    256         result = finishFloat64(exp, frac2, result.parts.sign);
     264        result = finishFloat64(exp, frac1, result.parts.sign);
    257265        return result;
    258266}
    259267
    260 /** Multiply two 64 bit numbers and return result in two parts
    261  * @param a first operand
    262  * @param b second operand
    263  * @param lo lower part from result
    264  * @param hi higher part of result
    265  */
    266 void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi)
     268/**
     269 * Multiply two quadruple-precision floats.
     270 *
     271 * @param a First input operand.
     272 * @param b Second input operand.
     273 * @return Result of multiplication.
     274 */
     275float128 mulFloat128(float128 a, float128 b)
    267276{
    268         uint64_t low, high, middle1, middle2;
    269         uint32_t alow, blow;
    270 
    271         alow = a & 0xFFFFFFFF;
    272         blow = b & 0xFFFFFFFF;
    273        
    274         a >>= 32;
    275         b >>= 32;
    276        
    277         low = ((uint64_t)alow) * blow;
    278         middle1 = a * blow;
    279         middle2 = alow * b;
    280         high = a * b;
    281 
    282         middle1 += middle2;
    283         high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
    284         middle1 <<= 32;
    285         low += middle1;
    286         high += (low < middle1);
    287         *lo = low;
    288         *hi = high;
    289        
    290         return;
     277        float128 result;
     278        uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
     279        int32_t exp;
     280
     281        result.parts.sign = a.parts.sign ^ b.parts.sign;
     282
     283        if (isFloat128NaN(a) || isFloat128NaN(b)) {
     284                /* TODO: fix SigNaNs */
     285                if (isFloat128SigNaN(a)) {
     286                        result.parts.frac_hi = a.parts.frac_hi;
     287                        result.parts.frac_lo = a.parts.frac_lo;
     288                        result.parts.exp = a.parts.exp;
     289                        return result;
     290                }
     291                if (isFloat128SigNaN(b)) { /* TODO: fix SigNaN */
     292                        result.parts.frac_hi = b.parts.frac_hi;
     293                        result.parts.frac_lo = b.parts.frac_lo;
     294                        result.parts.exp = b.parts.exp;
     295                        return result;
     296                }
     297                /* set NaN as result */
     298                result.binary.hi = FLOAT128_NAN_HI;
     299                result.binary.lo = FLOAT128_NAN_LO;
     300                return result;
     301        }
     302
     303        if (isFloat128Infinity(a)) {
     304                if (isFloat128Zero(b)) {
     305                        /* FIXME: zero * infinity */
     306                        result.binary.hi = FLOAT128_NAN_HI;
     307                        result.binary.lo = FLOAT128_NAN_LO;
     308                        return result;
     309                }
     310                result.parts.frac_hi = a.parts.frac_hi;
     311                result.parts.frac_lo = a.parts.frac_lo;
     312                result.parts.exp = a.parts.exp;
     313                return result;
     314        }
     315
     316        if (isFloat128Infinity(b)) {
     317                if (isFloat128Zero(a)) {
     318                        /* FIXME: zero * infinity */
     319                        result.binary.hi = FLOAT128_NAN_HI;
     320                        result.binary.lo = FLOAT128_NAN_LO;
     321                        return result;
     322                }
     323                result.parts.frac_hi = b.parts.frac_hi;
     324                result.parts.frac_lo = b.parts.frac_lo;
     325                result.parts.exp = b.parts.exp;
     326                return result;
     327        }
     328
     329        /* exp is signed so we can easy detect underflow */
     330        exp = a.parts.exp + b.parts.exp - FLOAT128_BIAS - 1;
     331
     332        frac1_hi = a.parts.frac_hi;
     333        frac1_lo = a.parts.frac_lo;
     334
     335        if (a.parts.exp > 0) {
     336                or128(frac1_hi, frac1_lo,
     337                FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     338                &frac1_hi, &frac1_lo);
     339        } else {
     340                ++exp;
     341        }
     342
     343        frac2_hi = b.parts.frac_hi;
     344        frac2_lo = b.parts.frac_lo;
     345
     346        if (b.parts.exp > 0) {
     347                or128(frac2_hi, frac2_lo,
     348                    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
     349                    &frac2_hi, &frac2_lo);
     350        } else {
     351                ++exp;
     352        }
     353
     354        lshift128(frac2_hi, frac2_lo,
     355            128 - FLOAT128_FRACTION_SIZE, &frac2_hi, &frac2_lo);
     356
     357        tmp_hi = frac1_hi;
     358        tmp_lo = frac1_lo;
     359        mul128(frac1_hi, frac1_lo, frac2_hi, frac2_lo,
     360            &frac1_hi, &frac1_lo, &frac2_hi, &frac2_lo);
     361        add128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo);
     362        frac2_hi |= (frac2_lo != 0x0ll);
     363
     364        if ((FLOAT128_HIDDEN_BIT_MASK_HI << 1) <= frac1_hi) {
     365                frac2_hi >>= 1;
     366                if (frac1_lo & 0x1ll) {
     367                        frac2_hi |= (0x1ull < 64);
     368                }
     369                rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
     370                ++exp;
     371        }
     372
     373        result = finishFloat128(exp, frac1_hi, frac1_lo, result.parts.sign, frac2_hi);
     374        return result;
    291375}
    292376
Note: See TracChangeset for help on using the changeset viewer.