Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • uspace/lib/softfloat/generic/mul.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 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.