Index: uspace/lib/softfloat/generic/div.c
===================================================================
--- uspace/lib/softfloat/generic/div.c	(revision 750636ac4a65360c44846b7681af6ae46854dd3a)
+++ uspace/lib/softfloat/generic/div.c	(revision 4449c6cb631ad8343054effc2668e54a01579700)
@@ -1,4 +1,5 @@
 /*
  * Copyright (c) 2005 Josef Cejka
+ * Copyright (c) 2011 Petr Koupy
  * All rights reserved.
  *
@@ -30,5 +31,5 @@
  * @{
  */
-/** @file
+/** @file Division functions.
  */
 
@@ -40,4 +41,11 @@
 #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,5 +108,4 @@
 		return result;
 	}
-
 	
 	afrac = a.parts.fraction;
@@ -110,13 +117,13 @@
 	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,6 +133,6 @@
 	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,8 +140,8 @@
 	}
 
-	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,6 +151,6 @@
 	
 	cfrac = (afrac << 32) / bfrac;
-	if ((  cfrac & 0x3F ) == 0) { 
-		cfrac |= ( bfrac * cfrac != afrac << 32 );
+	if ((cfrac & 0x3F) == 0) { 
+		cfrac |= (bfrac * cfrac != afrac << 32);
 	}
 	
@@ -151,9 +158,9 @@
 	
 	/* 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,8 +169,8 @@
 		++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,10 +188,9 @@
 		cfrac >>= 1;
 		while (cexp < 0) {
-			cexp ++;
+			cexp++;
 			cfrac >>= 1;
-		}
-		
+		}	
 	} else {
-		result.parts.exp = (uint32_t)cexp;
+		result.parts.exp = (uint32_t) cexp;
 	}
 	
@@ -194,4 +200,11 @@
 }
 
+/**
+ * Divide two double-precision floats.
+ *
+ * @param a Nominator.
+ * @param b Denominator.
+ * @return Result of division.
+ */
 float64 divFloat64(float64 a, float64 b) 
 {
@@ -200,9 +213,9 @@
 	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,5 +275,4 @@
 	}
 
-	
 	afrac = a.parts.fraction;
 	aexp = a.parts.exp;
@@ -275,9 +287,9 @@
 			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,6 +299,6 @@
 	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,8 +306,8 @@
 	}
 
-	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,18 +316,15 @@
 	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,38 +332,205 @@
 	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;
 }
