Index: uspace/lib/softfloat/generic/mul.c
===================================================================
--- uspace/lib/softfloat/generic/mul.c	(revision 750636ac4a65360c44846b7681af6ae46854dd3a)
+++ uspace/lib/softfloat/generic/mul.c	(revision 4863bb52d2de7a13a21dab7be2475975a848d49f)
@@ -1,4 +1,5 @@
 /*
  * Copyright (c) 2005 Josef Cejka
+ * Copyright (c) 2011 Petr Koupy
  * All rights reserved.
  *
@@ -30,5 +31,5 @@
  * @{
  */
-/** @file
+/** @file Multiplication functions.
  */
 
@@ -38,6 +39,10 @@
 #include <common.h>
 
-/** Multiply two 32 bit float numbers
- *
+/**
+ * Multiply two single-precision floats.
+ *
+ * @param a First input operand.
+ * @param b Second input operand.
+ * @return Result of multiplication.
  */
 float32 mulFloat32(float32 a, float32 b)
@@ -49,5 +54,5 @@
 	result.parts.sign = a.parts.sign ^ b.parts.sign;
 	
-	if (isFloat32NaN(a) || isFloat32NaN(b) ) {
+	if (isFloat32NaN(a) || isFloat32NaN(b)) {
 		/* TODO: fix SigNaNs */
 		if (isFloat32SigNaN(a)) {
@@ -55,14 +60,14 @@
 			result.parts.exp = a.parts.exp;
 			return result;
-		};
+		}
 		if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */
 			result.parts.fraction = b.parts.fraction;
 			result.parts.exp = b.parts.exp;
 			return result;
-		};
+		}
 		/* set NaN as result */
 		result.binary = FLOAT32_NAN;
 		return result;
-	};
+	}
 		
 	if (isFloat32Infinity(a)) { 
@@ -98,5 +103,5 @@
 		result.parts.sign = a.parts.sign ^ b.parts.sign;
 		return result;
-	};
+	}
 	
 	if (exp < 0) { 
@@ -106,5 +111,5 @@
 		result.parts.exp = 0x0;
 		return result;
-	};
+	}
 	
 	frac1 = a.parts.fraction;
@@ -113,5 +118,5 @@
 	} else {
 		++exp;
-	};
+	}
 	
 	frac2 = b.parts.fraction;
@@ -121,16 +126,16 @@
 	} else {
 		++exp;
-	};
+	}
 
 	frac1 <<= 1; /* one bit space for rounding */
 
 	frac1 = frac1 * frac2;
-/* round and return */
-	
-	while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) { 
-		/* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/
+
+	/* round and return */
+	while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 2)))) { 
+		/* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left) */
 		++exp;
 		frac1 >>= 1;
-	};
+	}
 
 	/* rounding */
@@ -141,7 +146,7 @@
 		++exp;
 		frac1 >>= 1;
-	};
-
-	if (exp >= FLOAT32_MAX_EXPONENT ) {	
+	}
+
+	if (exp >= FLOAT32_MAX_EXPONENT) {	
 		/* TODO: fix overflow */
 		/* return infinity*/
@@ -159,21 +164,24 @@
 			frac1 >>= 1;
 			++exp;
-		};
+		}
 		if (frac1 == 0) {
 			/* FIXME : underflow */
-		result.parts.exp = 0;
-		result.parts.fraction = 0;
-		return result;
-		};
-	};
+			result.parts.exp = 0;
+			result.parts.fraction = 0;
+			return result;
+		}
+	}
 	result.parts.exp = exp; 
-	result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);
+	result.parts.fraction = frac1 & ((1 << FLOAT32_FRACTION_SIZE) - 1);
 	
 	return result;	
-	
 }
 
-/** Multiply two 64 bit float numbers
- *
+/**
+ * Multiply two double-precision floats.
+ *
+ * @param a First input operand.
+ * @param b Second input operand.
+ * @return Result of multiplication.
  */
 float64 mulFloat64(float64 a, float64 b)
@@ -185,5 +193,5 @@
 	result.parts.sign = a.parts.sign ^ b.parts.sign;
 	
-	if (isFloat64NaN(a) || isFloat64NaN(b) ) {
+	if (isFloat64NaN(a) || isFloat64NaN(b)) {
 		/* TODO: fix SigNaNs */
 		if (isFloat64SigNaN(a)) {
@@ -191,14 +199,14 @@
 			result.parts.exp = a.parts.exp;
 			return result;
-		};
+		}
 		if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */
 			result.parts.fraction = b.parts.fraction;
 			result.parts.exp = b.parts.exp;
 			return result;
-		};
+		}
 		/* set NaN as result */
 		result.binary = FLOAT64_NAN;
 		return result;
-	};
+	}
 		
 	if (isFloat64Infinity(a)) { 
@@ -233,5 +241,5 @@
 	} else {
 		++exp;
-	};
+	}
 	
 	frac2 = b.parts.fraction;
@@ -241,52 +249,128 @@
 	} else {
 		++exp;
-	};
+	}
 
 	frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1);
 	frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2);
 
-	mul64integers(frac1, frac2, &frac1, &frac2);
-
-	frac2 |= (frac1 != 0);
-	if (frac2 & (0x1ll << 62)) {
-		frac2 <<= 1;
+	mul64(frac1, frac2, &frac1, &frac2);
+
+	frac1 |= (frac2 != 0);
+	if (frac1 & (0x1ll << 62)) {
+		frac1 <<= 1;
 		exp--;
 	}
 
-	result = finishFloat64(exp, frac2, result.parts.sign);
+	result = finishFloat64(exp, frac1, result.parts.sign);
 	return result;
 }
 
-/** Multiply two 64 bit numbers and return result in two parts
- * @param a first operand
- * @param b second operand
- * @param lo lower part from result
- * @param hi higher part of result
- */
-void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi)
+/**
+ * Multiply two quadruple-precision floats.
+ *
+ * @param a First input operand.
+ * @param b Second input operand.
+ * @return Result of multiplication.
+ */
+float128 mulFloat128(float128 a, float128 b)
 {
-	uint64_t low, high, middle1, middle2;
-	uint32_t alow, blow;
-
-	alow = a & 0xFFFFFFFF;
-	blow = b & 0xFFFFFFFF;
-	
-	a >>= 32;
-	b >>= 32;
-	
-	low = ((uint64_t)alow) * blow;
-	middle1 = a * blow;
-	middle2 = alow * b;
-	high = a * b;
-
-	middle1 += middle2;
-	high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32);
-	middle1 <<= 32;
-	low += middle1;
-	high += (low < middle1);
-	*lo = low;
-	*hi = high;
-	
-	return;
+	float128 result;
+	uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
+	int32_t exp;
+
+	result.parts.sign = a.parts.sign ^ b.parts.sign;
+
+	if (isFloat128NaN(a) || isFloat128NaN(b)) {
+		/* TODO: fix SigNaNs */
+		if (isFloat128SigNaN(a)) {
+			result.parts.frac_hi = a.parts.frac_hi;
+			result.parts.frac_lo = a.parts.frac_lo;
+			result.parts.exp = a.parts.exp;
+			return result;
+		}
+		if (isFloat128SigNaN(b)) { /* TODO: fix SigNaN */
+			result.parts.frac_hi = b.parts.frac_hi;
+			result.parts.frac_lo = b.parts.frac_lo;
+			result.parts.exp = b.parts.exp;
+			return result;
+		}
+		/* set NaN as result */
+		result.binary.hi = FLOAT128_NAN_HI;
+		result.binary.lo = FLOAT128_NAN_LO;
+		return result;
+	}
+
+	if (isFloat128Infinity(a)) {
+		if (isFloat128Zero(b)) {
+			/* FIXME: zero * infinity */
+			result.binary.hi = FLOAT128_NAN_HI;
+			result.binary.lo = FLOAT128_NAN_LO;
+			return result;
+		}
+		result.parts.frac_hi = a.parts.frac_hi;
+		result.parts.frac_lo = a.parts.frac_lo;
+		result.parts.exp = a.parts.exp;
+		return result;
+	}
+
+	if (isFloat128Infinity(b)) {
+		if (isFloat128Zero(a)) {
+			/* FIXME: zero * infinity */
+			result.binary.hi = FLOAT128_NAN_HI;
+			result.binary.lo = FLOAT128_NAN_LO;
+			return result;
+		}
+		result.parts.frac_hi = b.parts.frac_hi;
+		result.parts.frac_lo = b.parts.frac_lo;
+		result.parts.exp = b.parts.exp;
+		return result;
+	}
+
+	/* exp is signed so we can easy detect underflow */
+	exp = a.parts.exp + b.parts.exp - FLOAT128_BIAS - 1;
+
+	frac1_hi = a.parts.frac_hi;
+	frac1_lo = a.parts.frac_lo;
+
+	if (a.parts.exp > 0) {
+		or128(frac1_hi, frac1_lo,
+	        FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	        &frac1_hi, &frac1_lo);
+	} else {
+		++exp;
+	}
+
+	frac2_hi = b.parts.frac_hi;
+	frac2_lo = b.parts.frac_lo;
+
+	if (b.parts.exp > 0) {
+		or128(frac2_hi, frac2_lo,
+		    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+		    &frac2_hi, &frac2_lo);
+	} else {
+		++exp;
+	}
+
+	lshift128(frac2_hi, frac2_lo,
+	    128 - FLOAT128_FRACTION_SIZE, &frac2_hi, &frac2_lo);
+
+	tmp_hi = frac1_hi;
+	tmp_lo = frac1_lo;
+	mul128(frac1_hi, frac1_lo, frac2_hi, frac2_lo,
+	    &frac1_hi, &frac1_lo, &frac2_hi, &frac2_lo);
+	add128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo);
+	frac2_hi |= (frac2_lo != 0x0ll);
+
+	if ((FLOAT128_HIDDEN_BIT_MASK_HI << 1) <= frac1_hi) {
+		frac2_hi >>= 1;
+		if (frac1_lo & 0x1ll) {
+			frac2_hi |= (0x1ull < 64);
+		}
+		rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
+		++exp;
+	}
+
+	result = finishFloat128(exp, frac1_hi, frac1_lo, result.parts.sign, frac2_hi);
+	return result;
 }