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

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:

• uspace/lib/softfloat/generic/div.c (modified) (19 diffs)

uspace/lib/softfloat/generic/div.c

@@ -1 +1 @@
 /*
  * Copyright (c) 2005 Josef Cejka
+ * Copyright (c) 2011 Petr Koupy
  * All rights reserved.
  *
@@ -30 +31 @@
  * @{
  */
-/** @file
+/** @file Division functions.
  */
 
@@ -40 +41 @@
 #include <common.h>
 
+/**
+ * Divide two single-precision floats.
+ * 
+ * @param a Nominator.
+ * @param b Denominator.
+ * @return Result of division.
+ */
 float32 divFloat32(float32 a, float32 b) 
 {
@@ -100 +108 @@
                 return result;
         }
-
         
         afrac = a.parts.fraction;
@@ -110 +117 @@
         if (aexp == 0) {
                 if (afrac == 0) {
-                result.parts.exp = 0;
-                result.parts.fraction = 0;
-                return result;
-                }
+                        result.parts.exp = 0;
+                        result.parts.fraction = 0;
+                        return result;
+                }
+
                 /* normalize it*/
-                
                 afrac <<= 1;
-                        /* afrac is nonzero => it must stop */  
-                while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
+                /* afrac is nonzero => it must stop */  
+                while (!(afrac & FLOAT32_HIDDEN_BIT_MASK)) {
                         afrac <<= 1;
                         aexp--;
@@ -126 +133 @@
         if (bexp == 0) {
                 bfrac <<= 1;
-                        /* bfrac is nonzero => it must stop */  
-                while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
+                /* bfrac is nonzero => it must stop */  
+                while (!(bfrac & FLOAT32_HIDDEN_BIT_MASK)) {
                         bfrac <<= 1;
                         bexp--;
@@ -133 +140 @@
         }
 
-        afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
-        bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
-
-        if ( bfrac <= (afrac << 1) ) {
+        afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE - 1);
+        bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE);
+
+        if (bfrac <= (afrac << 1)) {
                 afrac >>= 1;
                 aexp++;
@@ -144 +151 @@
         
         cfrac = (afrac << 32) / bfrac;
-        if ((  cfrac & 0x3F ) == 0) { 
-                cfrac |= ( bfrac * cfrac != afrac << 32 );
+        if ((cfrac & 0x3F) == 0) { 
+                cfrac |= (bfrac * cfrac != afrac << 32);
         }
         
@@ -151 +158 @@
         
         /* find first nonzero digit and shift result and detect possibly underflow */
-        while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
+        while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)))) {
                 cexp--;
                 cfrac <<= 1;
-                        /* TODO: fix underflow */
-        };
+                /* TODO: fix underflow */
+        }
         
         cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
@@ -162 +169 @@
                 ++cexp;
                 cfrac >>= 1;
-                }       
+        }       
 
         /* check overflow */
-        if (cexp >= FLOAT32_MAX_EXPONENT ) {
+        if (cexp >= FLOAT32_MAX_EXPONENT) {
                 /* FIXME: overflow, return infinity */
                 result.parts.exp = FLOAT32_MAX_EXPONENT;
@@ -181 +188 @@
                 cfrac >>= 1;
                 while (cexp < 0) {
-                        cexp ++;
+                        cexp++;
                         cfrac >>= 1;
-                }
-                
+                }       
         } else {
-                result.parts.exp = (uint32_t)cexp;
+                result.parts.exp = (uint32_t) cexp;
         }
         
@@ -194 +200 @@
 }
 
+/**
+ * Divide two double-precision floats.
+ *
+ * @param a Nominator.
+ * @param b Denominator.
+ * @return Result of division.
+ */
 float64 divFloat64(float64 a, float64 b) 
 {
@@ -200 +213 @@
         uint64_t afrac, bfrac, cfrac; 
         uint64_t remlo, remhi;
+        uint64_t tmplo, tmphi;
         
         result.parts.sign = a.parts.sign ^ b.parts.sign;
         
         if (isFloat64NaN(a)) {
-                
                 if (isFloat64SigNaN(b)) {
                         /*FIXME: SigNaN*/
@@ -262 +275 @@
         }
 
-        
         afrac = a.parts.fraction;
         aexp = a.parts.exp;
@@ -275 +287 @@
                         return result;
                 }
+
                 /* normalize it*/
-                
                 aexp++;
-                        /* afrac is nonzero => it must stop */  
-                while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
+                /* afrac is nonzero => it must stop */  
+                while (!(afrac & FLOAT64_HIDDEN_BIT_MASK)) {
                         afrac <<= 1;
                         aexp--;
@@ -287 +299 @@
         if (bexp == 0) {
                 bexp++;
-                        /* bfrac is nonzero => it must stop */  
-                while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
+                /* bfrac is nonzero => it must stop */  
+                while (!(bfrac & FLOAT64_HIDDEN_BIT_MASK)) {
                         bfrac <<= 1;
                         bexp--;
@@ -294 +306 @@
         }
 
-        afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
-        bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
-
-        if ( bfrac <= (afrac << 1) ) {
+        afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 2);
+        bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 1);
+
+        if (bfrac <= (afrac << 1)) {
                 afrac >>= 1;
                 aexp++;
@@ -304 +316 @@
         cexp = aexp - bexp + FLOAT64_BIAS - 2; 
         
-        cfrac = divFloat64estim(afrac, bfrac);
-        
-        if ((  cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
-                mul64integers( bfrac, cfrac, &remlo, &remhi);
-                /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/      
-                remhi = afrac - remhi - ( remlo > 0);
-                remlo = - remlo;
+        cfrac = div128est(afrac, 0x0ll, bfrac);
+        
+        if ((cfrac & 0x1FF) <= 2) {
+                mul64(bfrac, cfrac, &tmphi, &tmplo);
+                sub128(afrac, 0x0ll, tmphi, tmplo, &remhi, &remlo);
                 
                 while ((int64_t) remhi < 0) {
                         cfrac--;
-                        remlo += bfrac;
-                        remhi += ( remlo < bfrac );
-                }
-                cfrac |= ( remlo != 0 );
+                        add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo);
+                }
+                cfrac |= (remlo != 0);
         }
         
@@ -323 +332 @@
         result = finishFloat64(cexp, cfrac, result.parts.sign);
         return result;
-
 }
 
-uint64_t divFloat64estim(uint64_t a, uint64_t b)
+/**
+ * Divide two quadruple-precision floats.
+ *
+ * @param a Nominator.
+ * @param b Denominator.
+ * @return Result of division.
+ */
+float128 divFloat128(float128 a, float128 b)
 {
-        uint64_t bhi;
-        uint64_t remhi, remlo;
-        uint64_t result;
-        
-        if ( b <= a ) {
-                return 0xFFFFFFFFFFFFFFFFull;
-        }
-        
-        bhi = b >> 32;
-        result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
-        mul64integers(b, result, &remlo, &remhi);
-        
-        remhi = a - remhi - (remlo > 0);
-        remlo = - remlo;
-
-        b <<= 32;
-        while ( (int64_t) remhi < 0 ) {
-                        result -= 0x1ll << 32;  
-                        remlo += b;
-                        remhi += bhi + ( remlo < b );
-                }
-        remhi = (remhi << 32) | (remlo >> 32);
-        if (( bhi << 32) <= remhi) {
-                result |= 0xFFFFFFFF;
-        } else {
-                result |= remhi / bhi;
-        }
-        
-        
+        float128 result;
+        int64_t aexp, bexp, cexp;
+        uint64_t afrac_hi, afrac_lo, bfrac_hi, bfrac_lo, cfrac_hi, cfrac_lo;
+        uint64_t shift_out;
+        uint64_t rem_hihi, rem_hilo, rem_lohi, rem_lolo;
+        uint64_t tmp_hihi, tmp_hilo, tmp_lohi, tmp_lolo;
+
+        result.parts.sign = a.parts.sign ^ b.parts.sign;
+
+        if (isFloat128NaN(a)) {
+                if (isFloat128SigNaN(b)) {
+                        /*FIXME: SigNaN*/
+                        return b;
+                }
+
+                if (isFloat128SigNaN(a)) {
+                        /*FIXME: SigNaN*/
+                }
+                /*NaN*/
+                return a;
+        }
+
+        if (isFloat128NaN(b)) {
+                if (isFloat128SigNaN(b)) {
+                        /*FIXME: SigNaN*/
+                }
+                /*NaN*/
+                return b;
+        }
+
+        if (isFloat128Infinity(a)) {
+                if (isFloat128Infinity(b) || isFloat128Zero(b)) {
+                        /*FIXME: inf / inf */
+                        result.binary.hi = FLOAT128_NAN_HI;
+                        result.binary.lo = FLOAT128_NAN_LO;
+                        return result;
+                }
+                /* inf / num */
+                result.parts.exp = a.parts.exp;
+                result.parts.frac_hi = a.parts.frac_hi;
+                result.parts.frac_lo = a.parts.frac_lo;
+                return result;
+        }
+
+        if (isFloat128Infinity(b)) {
+                if (isFloat128Zero(a)) {
+                        /* FIXME 0 / inf */
+                        result.parts.exp = 0;
+                        result.parts.frac_hi = 0;
+                        result.parts.frac_lo = 0;
+                        return result;
+                }
+                /* FIXME: num / inf*/
+                result.parts.exp = 0;
+                result.parts.frac_hi = 0;
+                result.parts.frac_lo = 0;
+                return result;
+        }
+
+        if (isFloat128Zero(b)) {
+                if (isFloat128Zero(a)) {
+                        /*FIXME: 0 / 0*/
+                        result.binary.hi = FLOAT128_NAN_HI;
+                        result.binary.lo = FLOAT128_NAN_LO;
+                        return result;
+                }
+                /* FIXME: division by zero */
+                result.parts.exp = 0;
+                result.parts.frac_hi = 0;
+                result.parts.frac_lo = 0;
+                return result;
+        }
+
+        afrac_hi = a.parts.frac_hi;
+        afrac_lo = a.parts.frac_lo;
+        aexp = a.parts.exp;
+        bfrac_hi = b.parts.frac_hi;
+        bfrac_lo = b.parts.frac_lo;
+        bexp = b.parts.exp;
+
+        /* denormalized numbers */
+        if (aexp == 0) {
+                if (eq128(afrac_hi, afrac_lo, 0x0ll, 0x0ll)) {
+                        result.parts.exp = 0;
+                        result.parts.frac_hi = 0;
+                        result.parts.frac_lo = 0;
+                        return result;
+                }
+
+                /* normalize it*/
+                aexp++;
+                /* afrac is nonzero => it must stop */
+                and128(afrac_hi, afrac_lo,
+                    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+                    &tmp_hihi, &tmp_lolo);
+                while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
+                        lshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
+                        aexp--;
+                }
+        }
+
+        if (bexp == 0) {
+                bexp++;
+                /* bfrac is nonzero => it must stop */
+                and128(bfrac_hi, bfrac_lo,
+                    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+                    &tmp_hihi, &tmp_lolo);
+                while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
+                        lshift128(bfrac_hi, bfrac_lo, 1, &bfrac_hi, &bfrac_lo);
+                        bexp--;
+                }
+        }
+
+        or128(afrac_hi, afrac_lo,
+            FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+            &afrac_hi, &afrac_lo);
+        lshift128(afrac_hi, afrac_lo,
+            (128 - FLOAT128_FRACTION_SIZE - 1), &afrac_hi, &afrac_lo);
+        or128(bfrac_hi, bfrac_lo,
+            FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+            &bfrac_hi, &bfrac_lo);
+        lshift128(bfrac_hi, bfrac_lo,
+            (128 - FLOAT128_FRACTION_SIZE - 1), &bfrac_hi, &bfrac_lo);
+
+        if (le128(bfrac_hi, bfrac_lo, afrac_hi, afrac_lo)) {
+                rshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
+                aexp++;
+        }
+
+        cexp = aexp - bexp + FLOAT128_BIAS - 2;
+
+        cfrac_hi = div128est(afrac_hi, afrac_lo, bfrac_hi);
+
+        mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_hi,
+            &tmp_lolo /* dummy */, &tmp_hihi, &tmp_hilo, &tmp_lohi);
+
+        /* sub192(afrac_hi, afrac_lo, 0, 
+         *     tmp_hihi, tmp_hilo, tmp_lohi
+         *     &rem_hihi, &rem_hilo, &rem_lohi); */
+        sub128(afrac_hi, afrac_lo, tmp_hihi, tmp_hilo, &rem_hihi, &rem_hilo);
+        if (tmp_lohi > 0) {
+                sub128(rem_hihi, rem_hilo, 0x0ll, 0x1ll, &rem_hihi, &rem_hilo);
+        }
+        rem_lohi = -tmp_lohi;
+
+        while ((int64_t) rem_hihi < 0) {
+                --cfrac_hi;
+                /* add192(rem_hihi, rem_hilo, rem_lohi, 
+                 *     0, bfrac_hi, bfrac_lo,
+                 *     &rem_hihi, &rem_hilo, &rem_lohi); */
+                add128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo, &rem_hilo, &rem_lohi);
+                if (lt128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo)) {
+                        ++rem_hihi;
+                }
+        }
+
+        cfrac_lo = div128est(rem_hilo, rem_lohi, bfrac_lo);
+
+        if ((cfrac_lo & 0x3FFF) <= 4) {
+                mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_lo,
+                &tmp_hihi /* dummy */, &tmp_hilo, &tmp_lohi, &tmp_lolo);
+
+                /* sub192(rem_hilo, rem_lohi, 0,
+                 *     tmp_hilo, tmp_lohi, tmp_lolo,
+                 *     &rem_hilo, &rem_lohi, &rem_lolo); */
+                sub128(rem_hilo, rem_lohi, tmp_hilo, tmp_lohi, &rem_hilo, &rem_lohi);
+                if (tmp_lolo > 0) {
+                        sub128(rem_hilo, rem_lohi, 0x0ll, 0x1ll, &rem_hilo, &rem_lohi);
+                }
+                rem_lolo = -tmp_lolo;
+
+                while ((int64_t) rem_hilo < 0) {
+                        --cfrac_lo;
+                        /* add192(rem_hilo, rem_lohi, rem_lolo,
+                         *     0, bfrac_hi, bfrac_lo,
+                         *     &rem_hilo, &rem_lohi, &rem_lolo); */
+                        add128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo, &rem_lohi, &rem_lolo);
+                        if (lt128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo)) {
+                                ++rem_hilo;
+                        }
+                }
+
+                cfrac_lo |= ((rem_hilo | rem_lohi | rem_lolo) != 0 );
+        }
+
+        shift_out = cfrac_lo << (64 - (128 - FLOAT128_FRACTION_SIZE - 1));
+        rshift128(cfrac_hi, cfrac_lo, (128 - FLOAT128_FRACTION_SIZE - 1),
+            &cfrac_hi, &cfrac_lo);
+
+        result = finishFloat128(cexp, cfrac_hi, cfrac_lo, result.parts.sign, shift_out);
         return result;
 }