Index: uspace/lib/softfloat/generic/sub.c
===================================================================
--- uspace/lib/softfloat/generic/sub.c	(revision 750636ac4a65360c44846b7681af6ae46854dd3a)
+++ uspace/lib/softfloat/generic/sub.c	(revision 3a11f175908a16f3fd12702eb1389589bb9a14c7)
@@ -1,4 +1,5 @@
 /*
  * Copyright (c) 2005 Josef Cejka
+ * Copyright (c) 2011 Petr Koupy
  * All rights reserved.
  *
@@ -30,5 +31,5 @@
  * @{
  */
-/** @file
+/** @file Substraction functions.
  */
 
@@ -36,6 +37,12 @@
 #include <sub.h>
 #include <comparison.h>
-
-/** Subtract two float32 numbers with same signs
+#include <common.h>
+
+/**
+ * Subtract two single-precision floats with the same signs.
+ *
+ * @param a First input operand.
+ * @param b Second input operand.
+ * @return Result of substraction.
  */
 float32 subFloat32(float32 a, float32 b)
@@ -52,7 +59,7 @@
 			/* TODO: fix SigNaN */
 			if (isFloat32SigNaN(b)) {
-			};
-			return b;
-		};
+			}
+			return b;
+		}
 		
 		if (b.parts.exp == FLOAT32_MAX_EXPONENT) { 
@@ -72,7 +79,7 @@
 			/* TODO: fix SigNaN */
 			if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) {
-			};
-			return a;
-		};
+			}
+			return a;
+		}
 		
 		if (a.parts.exp == FLOAT32_MAX_EXPONENT) { 
@@ -82,5 +89,5 @@
 				result.binary = FLOAT32_NAN;
 				return result;
-			};
+			}
 			return a;
 		}
@@ -92,16 +99,16 @@
 		frac2 = b.parts.fraction;
 		exp2 = b.parts.exp;	
-	};
+	}
 	
 	if (exp1 == 0) {
 		/* both are denormalized */
-		result.parts.fraction = frac1-frac2;
+		result.parts.fraction = frac1 - frac2;
 		if (result.parts.fraction > frac1) {
 			/* TODO: underflow exception */
 			return result;
-		};
+		}
 		result.parts.exp = 0;
 		return result;
-	};
+	}
 
 	/* add hidden bit */
@@ -114,5 +121,5 @@
 		/* normalized */
 		frac2 |= FLOAT32_HIDDEN_BIT_MASK;
-	};
+	}
 	
 	/* create some space for rounding */
@@ -121,8 +128,9 @@
 	
 	if (expdiff > FLOAT32_FRACTION_SIZE + 1) {
-	     goto done;	
-	     };
+		goto done;
+	}
 	
 	frac1 = frac1 - (frac2 >> expdiff);
+
 done:
 	/* TODO: find first nonzero digit and shift result and detect possibly underflow */
@@ -130,6 +138,6 @@
 		--exp1;
 		frac1 <<= 1;
-			/* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
-	};
+		/* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
+	}
 	
 	/* rounding - if first bit after fraction is set then round up */
@@ -139,7 +147,7 @@
 		++exp1;
 		frac1 >>= 1;
-	};
-	
-	/*Clear hidden bit and shift */
+	}
+	
+	/* Clear hidden bit and shift */
 	result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); 
 	result.parts.exp = exp1;
@@ -148,5 +156,10 @@
 }
 
-/** Subtract two float64 numbers with same signs
+/**
+ * Subtract two double-precision floats with the same signs.
+ *
+ * @param a First input operand.
+ * @param b Second input operand.
+ * @return Result of substraction.
  */
 float64 subFloat64(float64 a, float64 b)
@@ -164,7 +177,7 @@
 			/* TODO: fix SigNaN */
 			if (isFloat64SigNaN(b)) {
-			};
-			return b;
-		};
+			}
+			return b;
+		}
 		
 		if (b.parts.exp == FLOAT64_MAX_EXPONENT) { 
@@ -184,7 +197,7 @@
 			/* TODO: fix SigNaN */
 			if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) {
-			};
-			return a;
-		};
+			}
+			return a;
+		}
 		
 		if (a.parts.exp == FLOAT64_MAX_EXPONENT) { 
@@ -194,5 +207,5 @@
 				result.binary = FLOAT64_NAN;
 				return result;
-			};
+			}
 			return a;
 		}
@@ -204,5 +217,5 @@
 		frac2 = b.parts.fraction;
 		exp2 = b.parts.exp;	
-	};
+	}
 	
 	if (exp1 == 0) {
@@ -212,8 +225,8 @@
 			/* TODO: underflow exception */
 			return result;
-		};
+		}
 		result.parts.exp = 0;
 		return result;
-	};
+	}
 
 	/* add hidden bit */
@@ -226,5 +239,5 @@
 		/* normalized */
 		frac2 |= FLOAT64_HIDDEN_BIT_MASK;
-	};
+	}
 	
 	/* create some space for rounding */
@@ -233,8 +246,9 @@
 	
 	if (expdiff > FLOAT64_FRACTION_SIZE + 1) {
-	     goto done;	
-	     };
+		goto done;
+	}
 	
 	frac1 = frac1 - (frac2 >> expdiff);
+
 done:
 	/* TODO: find first nonzero digit and shift result and detect possibly underflow */
@@ -242,6 +256,6 @@
 		--exp1;
 		frac1 <<= 1;
-			/* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
-	};
+		/* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
+	}
 	
 	/* rounding - if first bit after fraction is set then round up */
@@ -251,7 +265,7 @@
 		++exp1;
 		frac1 >>= 1;
-	};
-	
-	/*Clear hidden bit and shift */
+	}
+	
+	/* Clear hidden bit and shift */
 	result.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK)); 
 	result.parts.exp = exp1;
@@ -260,4 +274,152 @@
 }
 
+/**
+ * Subtract two quadruple-precision floats with the same signs.
+ *
+ * @param a First input operand.
+ * @param b Second input operand.
+ * @return Result of substraction.
+ */
+float128 subFloat128(float128 a, float128 b)
+{
+	int expdiff;
+	uint32_t exp1, exp2;
+	uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
+	float128 result;
+
+	result.binary.hi = 0;
+	result.binary.lo = 0;
+
+	expdiff = a.parts.exp - b.parts.exp;
+	if ((expdiff < 0 ) || ((expdiff == 0) &&
+	    lt128(a.parts.frac_hi, a.parts.frac_lo, b.parts.frac_hi, b.parts.frac_lo))) {
+		if (isFloat128NaN(b)) {
+			/* TODO: fix SigNaN */
+			if (isFloat128SigNaN(b)) {
+			}
+			return b;
+		}
+
+		if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
+			b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */
+			return b;
+		}
+
+		result.parts.sign = !a.parts.sign;
+
+		frac1_hi = b.parts.frac_hi;
+		frac1_lo = b.parts.frac_lo;
+		exp1 = b.parts.exp;
+		frac2_hi = a.parts.frac_hi;
+		frac2_lo = a.parts.frac_lo;
+		exp2 = a.parts.exp;
+		expdiff *= -1;
+	} else {
+		if (isFloat128NaN(a)) {
+			/* TODO: fix SigNaN */
+			if (isFloat128SigNaN(a) || isFloat128SigNaN(b)) {
+			}
+			return a;
+		}
+
+		if (a.parts.exp == FLOAT128_MAX_EXPONENT) {
+			if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
+				/* inf - inf => nan */
+				/* TODO: fix exception */
+				result.binary.hi = FLOAT128_NAN_HI;
+				result.binary.lo = FLOAT128_NAN_LO;
+				return result;
+			}
+			return a;
+		}
+
+		result.parts.sign = a.parts.sign;
+
+		frac1_hi = a.parts.frac_hi;
+		frac1_lo = a.parts.frac_lo;
+		exp1 = a.parts.exp;
+		frac2_hi = b.parts.frac_hi;
+		frac2_lo = b.parts.frac_lo;
+		exp2 = b.parts.exp;
+	}
+
+	if (exp1 == 0) {
+		/* both are denormalized */
+		sub128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &tmp_hi, &tmp_lo);
+		result.parts.frac_hi = tmp_hi;
+		result.parts.frac_lo = tmp_lo;
+		if (lt128(frac1_hi, frac1_lo, result.parts.frac_hi, result.parts.frac_lo)) {
+			/* TODO: underflow exception */
+			return result;
+		}
+		result.parts.exp = 0;
+		return result;
+	}
+
+	/* add hidden bit */
+	or128(frac1_hi, frac1_lo,
+	    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	    &frac1_hi, &frac1_lo);
+
+	if (exp2 == 0) {
+		/* denormalized */
+		--expdiff;
+	} else {
+		/* normalized */
+		or128(frac2_hi, frac2_lo,
+		    FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+		    &frac2_hi, &frac2_lo);
+	}
+
+	/* create some space for rounding */
+	lshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
+	lshift128(frac2_hi, frac2_lo, 6, &frac2_hi, &frac2_lo);
+
+	if (expdiff > FLOAT128_FRACTION_SIZE + 1) {
+		goto done;
+	}
+
+	rshift128(frac2_hi, frac2_lo, expdiff, &tmp_hi, &tmp_lo);
+	sub128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo);
+
+done:
+	/* TODO: find first nonzero digit and shift result and detect possibly underflow */
+	lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6,
+	    &tmp_hi, &tmp_lo);
+	and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	while ((exp1 > 0) && (!lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo))) {
+		--exp1;
+		lshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
+		/* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
+
+		lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6,
+		    &tmp_hi, &tmp_lo);
+		and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	}
+
+	/* rounding - if first bit after fraction is set then round up */
+	add128(frac1_hi, frac1_lo, 0x0ll, 0x20ll, &frac1_hi, &frac1_lo);
+
+	lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7,
+	   &tmp_hi, &tmp_lo);
+	and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
+		++exp1;
+		rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
+	}
+
+	/* Clear hidden bit and shift */
+	rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
+	not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	    &tmp_hi, &tmp_lo);
+	and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	result.parts.frac_hi = tmp_hi;
+	result.parts.frac_lo = tmp_lo;
+
+	result.parts.exp = exp1;
+
+	return result;
+}
+
 /** @}
  */
