WebSVN – HelenOS-historic – Diff – /uspace/trunk/softfloat/generic/div.c

  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 #include<sftypes.h>
 #include<add.h>
+#include<div.h>
 #include<comparison.h>
+#include<mul.h>
 float32 divFloat32(float32 a, float32 b)
+{
     float32 result;
     __s32 aexp, bexp, cexp;
         cfrac |= ( bfrac * cfrac != afrac << 32 );
+    }
     /* pack and round */
-    /* TODO: find first nonzero digit and shift result and detect possibly underflow */
+    /* find first nonzero digit and shift result and detect possibly underflow */
     while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
         cexp--;
         cfrac <<= 1;
             /* TODO: fix underflow */
     };
     cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
     if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
         ++cexp;
         cfrac >>= 1;
     result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
     return result;
+}
+float64 divFloat64(float64 a, float64 b)
+{
+    float64 result;
+    __s32 aexp, bexp, cexp;
+    __u64 afrac, bfrac, cfrac;
+    __u64 remlo, remhi;
+    result.parts.sign = a.parts.sign ^ b.parts.sign;
+    if (isFloat64NaN(a)) {
+        if (isFloat64SigNaN(a)) {
+            /*FIXME: SigNaN*/
+        }
+        /*NaN*/
+        return a;
+    }
+    if (isFloat64NaN(b)) {
+        if (isFloat64SigNaN(b)) {
+            /*FIXME: SigNaN*/
+        }
+        /*NaN*/
+        return b;
+    }
+    if (isFloat64Infinity(a)) {
+        if (isFloat64Infinity(b)) {
+            /*FIXME: inf / inf */
+            result.binary = FLOAT64_NAN;
+            return result;
+        }
+        /* inf / num */
+        result.parts.exp = a.parts.exp;
+        result.parts.fraction = a.parts.fraction;
+        return result;
+    }
+    if (isFloat64Infinity(b)) {
+        if (isFloat64Zero(a)) {
+            /* FIXME 0 / inf */
+            result.parts.exp = 0;
+            result.parts.fraction = 0;
+            return result;
+        }
+        /* FIXME: num / inf*/
+        result.parts.exp = 0;
+        result.parts.fraction = 0;
+        return result;
+    }
+    if (isFloat64Zero(b)) {
+        if (isFloat64Zero(a)) {
+            /*FIXME: 0 / 0*/
+            result.binary = FLOAT64_NAN;
+            return result;
+        }
+        /* FIXME: division by zero */
+        result.parts.exp = 0;
+        result.parts.fraction = 0;
+        return result;
+    }
+    afrac = a.parts.fraction;
+    aexp = a.parts.exp;
+    bfrac = b.parts.fraction;
+    bexp = b.parts.exp;
+    /* denormalized numbers */
+    if (aexp == 0) {
+        if (afrac == 0) {
+        result.parts.exp = 0;
+        result.parts.fraction = 0;
+        return result;
+        }
+        /* normalize it*/
+        afrac <<= 1;
+            /* afrac is nonzero => it must stop */
+        while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
+            afrac <<= 1;
+            aexp--;
+        }
+    }
+    if (bexp == 0) {
+        bfrac <<= 1;
+            /* bfrac is nonzero => it must stop */
+        while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
+            bfrac <<= 1;
+            bexp--;
+        }
+    }
+    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++;
+    }
+    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;
+        while ((__s64) remhi < 0) {
+            cfrac--;
+            remlo += bfrac;
+            remhi += ( remlo < bfrac );
+        }
+        cfrac |= ( remlo != 0 );
+    }
+    /* pack and round */
+    /* find first nonzero digit and shift result and detect possibly underflow */
+    while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ) )))) {
+        cexp--;
+        cfrac <<= 1;
+            /* TODO: fix underflow */
+    };
+    cfrac >>= 1;
+    ++cexp;
+    cfrac += (0x1 << (64 - FLOAT64_FRACTION_SIZE - 3));
+    if (cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ))) {
+        ++cexp;
+        cfrac >>= 1;
+        }
+    /* check overflow */
+    if (cexp >= FLOAT64_MAX_EXPONENT ) {
+        /* FIXME: overflow, return infinity */
+        result.parts.exp = FLOAT64_MAX_EXPONENT;
+        result.parts.fraction = 0;
+        return result;
+    }
+    if (cexp < 0) {
+        /* FIXME: underflow */
+        result.parts.exp = 0;
+        if ((cexp + FLOAT64_FRACTION_SIZE) < 0) {
+            result.parts.fraction = 0;
+            return result;
+        }
+        cfrac >>= 1;
+        while (cexp < 0) {
+            cexp ++;
+            cfrac >>= 1;
+        }
+        return result;
+    } else {
+        cexp ++; /*normalized*/
+        result.parts.exp = (__u32)cexp;
+    }
+    result.parts.fraction = ((cfrac >>(64 - FLOAT64_FRACTION_SIZE - 2 ) ) & (~FLOAT64_HIDDEN_BIT_MASK));
+    return result;
+}
+__u64 divFloat64estim(__u64 a, __u64 b)
+{
+    __u64 bhi;
+    __u64 remhi, remlo;
+    __u64 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 ( (__s64) 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;
+    }
+    return result;
+}

Subversion Repositories HelenOS-historic

(root)/uspace/trunk/softfloat/generic/div.c – Rev 804 → 828