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

 /*
  * Copyright (C) 2005 Josef Cejka
  * All rights reserved.
+ *
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions
  * are met:
+ *
  * - Redistributions of source code must retain the above copyright
  *   notice, this list of conditions and the following disclaimer.
  * - Redistributions in binary form must reproduce the above copyright
  *   notice, this list of conditions and the following disclaimer in the
  *   documentation and/or other materials provided with the distribution.
  * - The name of the author may not be used to endorse or promote products
  *   derived from this software without specific prior written permission.
+ *
  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  * 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;
     __u64 afrac, bfrac, cfrac;
     result.parts.sign = a.parts.sign ^ b.parts.sign;
     if (isFloat32NaN(a)) {
         if (isFloat32SigNaN(a)) {
             /*FIXME: SigNaN*/
+        }
         /*NaN*/
         return a;
+    }
     if (isFloat32NaN(b)) {
         if (isFloat32SigNaN(b)) {
             /*FIXME: SigNaN*/
+        }
         /*NaN*/
         return b;
+    }
     if (isFloat32Infinity(a)) {
         if (isFloat32Infinity(b)) {
             /*FIXME: inf / inf */
             result.binary = FLOAT32_NAN;
             return result;
+        }
         /* inf / num */
         result.parts.exp = a.parts.exp;
         result.parts.fraction = a.parts.fraction;
         return result;
+    }
     if (isFloat32Infinity(b)) {
         if (isFloat32Zero(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 (isFloat32Zero(b)) {
         if (isFloat32Zero(a)) {
             /*FIXME: 0 / 0*/
             result.binary = FLOAT32_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 & FLOAT32_HIDDEN_BIT_MASK) ) {
             afrac <<= 1;
             aexp--;
+        }
+    }
     if (bexp == 0) {
         bfrac <<= 1;
             /* bfrac is nonzero => it must stop */
         while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
             bfrac <<= 1;
             bexp--;
+        }
+    }
     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++;
+    }
     cexp = aexp - bexp + FLOAT32_BIAS - 2;
     cfrac = (afrac << 32) / bfrac;
     if ((  cfrac & 0x3F ) == 0) {
         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;
+        }
     /* check overflow */
     if (cexp >= FLOAT32_MAX_EXPONENT ) {
         /* FIXME: overflow, return infinity */
         result.parts.exp = FLOAT32_MAX_EXPONENT;
         result.parts.fraction = 0;
         return result;
+    }
     if (cexp < 0) {
         /* FIXME: underflow */
         result.parts.exp = 0;
         if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
             result.parts.fraction = 0;
             return result;
+        }
         cfrac >>= 1;
         while (cexp < 0) {
             cexp ++;
             cfrac >>= 1;
+        }
     } else {
         result.parts.exp = (__u32)cexp;
+    }
     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