Index: uspace/lib/softfloat/generic/common.c
===================================================================
--- uspace/lib/softfloat/generic/common.c	(revision 750636ac4a65360c44846b7681af6ae46854dd3a)
+++ uspace/lib/softfloat/generic/common.c	(revision 128c78b842343a9ece0922104cf36feb7da16be9)
@@ -1,4 +1,5 @@
 /*
  * Copyright (c) 2005 Josef Cejka
+ * Copyright (c) 2011 Petr Koupy
  * All rights reserved.
  *
@@ -30,5 +31,5 @@
  * @{
  */
-/** @file
+/** @file Common helper operations.
  */
 
@@ -36,5 +37,5 @@
 #include <common.h>
 
-/* Table for fast leading zeroes counting */
+/* Table for fast leading zeroes counting. */
 char zeroTable[256] = {
 	8, 7, 7, 6, 6, 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, \
@@ -56,11 +57,12 @@
 };
 
-
-
-/** Take fraction shifted by 10 bits to left, round it, normalize it and detect exceptions
- * @param cexp exponent with bias
- * @param cfrac fraction shifted 10 places left with added hidden bit
- * @param sign
- * @return valied float64
+/** 
+ * Take fraction shifted by 10 bits to the left, round it, normalize it
+ * and detect exceptions
+ * 
+ * @param cexp Exponent with bias.
+ * @param cfrac Fraction shifted 10 bits to the left with added hidden bit.
+ * @param sign Resulting sign.
+ * @return Finished double-precision float.
  */
 float64 finishFloat64(int32_t cexp, uint64_t cfrac, char sign)
@@ -71,11 +73,13 @@
 
 	/* 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 ) )))) {
+	while ((cexp > 0) && (cfrac) &&
+	    (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1))))) {
 		cexp--; 
 		cfrac <<= 1;
-			/* TODO: fix underflow */
-	};
-	
-	if ((cexp < 0) || ( cexp == 0 && (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1)))))) {
+		/* TODO: fix underflow */
+	}
+	
+	if ((cexp < 0) || (cexp == 0 &&
+	    (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1)))))) {
 		/* FIXME: underflow */
 		result.parts.exp = 0;
@@ -93,6 +97,6 @@
 		
 		if (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1)))) {
-			
-			result.parts.fraction = ((cfrac >>(64 - FLOAT64_FRACTION_SIZE - 2) ) & (~FLOAT64_HIDDEN_BIT_MASK)); 
+			result.parts.fraction =
+			    ((cfrac >> (64 - FLOAT64_FRACTION_SIZE - 2)) & (~FLOAT64_HIDDEN_BIT_MASK));
 			return result;
 		}	
@@ -103,5 +107,5 @@
 	++cexp;
 
-	if (cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ))) {
+	if (cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1))) {
 		++cexp;
 		cfrac >>= 1;
@@ -109,5 +113,5 @@
 
 	/* check overflow */
-	if (cexp >= FLOAT64_MAX_EXPONENT ) {
+	if (cexp >= FLOAT64_MAX_EXPONENT) {
 		/* FIXME: overflow, return infinity */
 		result.parts.exp = FLOAT64_MAX_EXPONENT;
@@ -116,19 +120,156 @@
 	}
 
-	result.parts.exp = (uint32_t)cexp;
-	
-	result.parts.fraction = ((cfrac >>(64 - FLOAT64_FRACTION_SIZE - 2 ) ) & (~FLOAT64_HIDDEN_BIT_MASK)); 
+	result.parts.exp = (uint32_t) cexp;
+	
+	result.parts.fraction = 
+	    ((cfrac >> (64 - FLOAT64_FRACTION_SIZE - 2)) & (~FLOAT64_HIDDEN_BIT_MASK));
 	
 	return result;	
 }
 
-/** Counts leading zeroes in 64bit unsigned integer
- * @param i 
+/**
+ * Take fraction, round it, normalize it and detect exceptions
+ * 
+ * @param cexp Exponent with bias.
+ * @param cfrac_hi High part of the fraction shifted 14 bits to the left
+ *     with added hidden bit.
+ * @param cfrac_lo Low part of the fraction shifted 14 bits to the left
+ *     with added hidden bit.
+ * @param sign Resulting sign.
+ * @param shift_out Bits right-shifted out from fraction by the caller.
+ * @return Finished quadruple-precision float.
+ */
+float128 finishFloat128(int32_t cexp, uint64_t cfrac_hi, uint64_t cfrac_lo, 
+    char sign, uint64_t shift_out)
+{
+	float128 result;
+	uint64_t tmp_hi, tmp_lo;
+
+	result.parts.sign = sign;
+
+	/* find first nonzero digit and shift result and detect possibly underflow */
+	lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	    1, &tmp_hi, &tmp_lo);
+	and128(cfrac_hi, cfrac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	while ((cexp > 0) && (lt128(0x0ll, 0x0ll, cfrac_hi, cfrac_lo)) &&
+	    (!lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo))) {
+		cexp--;
+		lshift128(cfrac_hi, cfrac_lo, 1, &cfrac_hi, &cfrac_lo);
+		/* TODO: fix underflow */
+
+		lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+		    1, &tmp_hi, &tmp_lo);
+		and128(cfrac_hi, cfrac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	}
+
+	lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	    1, &tmp_hi, &tmp_lo);
+	and128(cfrac_hi, cfrac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	if ((cexp < 0) || (cexp == 0 &&
+	    (!lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)))) {
+		/* FIXME: underflow */
+		result.parts.exp = 0;
+		if ((cexp + FLOAT128_FRACTION_SIZE + 1) < 0) { /* +1 is place for rounding */
+			result.parts.frac_hi = 0x0ll;
+			result.parts.frac_lo = 0x0ll;
+			return result;
+		}
+
+		while (cexp < 0) {
+			cexp++;
+			rshift128(cfrac_hi, cfrac_lo, 1, &cfrac_hi, &cfrac_lo);
+		}
+
+		if (shift_out & (0x1ull < 64)) {
+			add128(cfrac_hi, cfrac_lo, 0x0ll, 0x1ll, &cfrac_hi, &cfrac_lo);
+		}
+
+		lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+		    1, &tmp_hi, &tmp_lo);
+		and128(cfrac_hi, cfrac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+		if (!lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
+			not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+			    &tmp_hi, &tmp_lo);
+			and128(cfrac_hi, cfrac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+			result.parts.frac_hi = tmp_hi;
+			result.parts.frac_lo = tmp_lo;
+			return result;
+		}
+	} else {
+		if (shift_out & (0x1ull < 64)) {
+			add128(cfrac_hi, cfrac_lo, 0x0ll, 0x1ll, &cfrac_hi, &cfrac_lo);
+		}
+	}
+
+	++cexp;
+
+	lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	    1, &tmp_hi, &tmp_lo);
+	and128(cfrac_hi, cfrac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
+		++cexp;
+		rshift128(cfrac_hi, cfrac_lo, 1, &cfrac_hi, &cfrac_lo);
+	}
+
+	/* check overflow */
+	if (cexp >= FLOAT128_MAX_EXPONENT) {
+		/* FIXME: overflow, return infinity */
+		result.parts.exp = FLOAT128_MAX_EXPONENT;
+		result.parts.frac_hi = 0x0ll;
+		result.parts.frac_lo = 0x0ll;
+		return result;
+	}
+
+	result.parts.exp = (uint32_t) cexp;
+
+	not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	    &tmp_hi, &tmp_lo);
+	and128(cfrac_hi, cfrac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	result.parts.frac_hi = tmp_hi;
+	result.parts.frac_lo = tmp_lo;
+
+	return result;	
+}
+
+/**
+ * Counts leading zeroes in byte.
+ *
+ * @param i Byte for which to count leading zeroes.
+ * @return Number of detected leading zeroes.
+ */
+int countZeroes8(uint8_t i)
+{
+	return zeroTable[i];
+}
+
+/** 
+ * Counts leading zeroes in 32bit unsigned integer.
+ *
+ * @param i Integer for which to count leading zeroes.
+ * @return Number of detected leading zeroes.
+ */
+int countZeroes32(uint32_t i)
+{
+	int j;
+	for (j = 0; j < 32; j += 8) {
+		if (i & (0xFF << (24 - j))) {
+			return (j + countZeroes8(i >> (24 - j)));
+		}
+	}
+
+	return 32;
+}
+
+/**
+ * Counts leading zeroes in 64bit unsigned integer.
+ *
+ * @param i Integer for which to count leading zeroes.
+ * @return Number of detected leading zeroes.
  */
 int countZeroes64(uint64_t i)
 {
 	int j;
-	for (j =0; j < 64; j += 8) {
-		if ( i & (0xFFll << (56 - j))) {
+	for (j = 0; j < 64; j += 8) {
+		if (i & (0xFFll << (56 - j))) {
 			return (j + countZeroes8(i >> (56 - j)));
 		}
@@ -138,75 +279,410 @@
 }
 
-/** Counts leading zeroes in 32bit unsigned integer
- * @param i 
- */
-int countZeroes32(uint32_t i)
-{
-	int j;
-	for (j =0; j < 32; j += 8) {
-		if ( i & (0xFF << (24 - j))) {
-			return (j + countZeroes8(i >> (24 - j)));
-		}
-	}
-
-	return 32;
-}
-
-/** Counts leading zeroes in byte
- * @param i 
- */
-int countZeroes8(uint8_t i)
-{
-	return zeroTable[i];
-}
-
-/** Round and normalize number expressed by exponent and fraction with first bit (equal to hidden bit) at 30. bit
- * @param exp exponent 
- * @param fraction part with hidden bit shifted to 30. bit
+/**
+ * Round and normalize number expressed by exponent and fraction with
+ * first bit (equal to hidden bit) at 30th bit.
+ *
+ * @param exp Exponent part.
+ * @param fraction Fraction with hidden bit shifted to 30th bit.
  */
 void roundFloat32(int32_t *exp, uint32_t *fraction)
 {
 	/* rounding - if first bit after fraction is set then round up */
-	(*fraction) += (0x1 << 6);
-	
-	if ((*fraction) & (FLOAT32_HIDDEN_BIT_MASK << 8)) { 
+	(*fraction) += (0x1 << (32 - FLOAT32_FRACTION_SIZE - 3));
+	
+	if ((*fraction) & 
+	    (FLOAT32_HIDDEN_BIT_MASK << (32 - FLOAT32_FRACTION_SIZE - 1))) {
 		/* rounding overflow */
 		++(*exp);
 		(*fraction) >>= 1;
-	};
-	
-	if (((*exp) >= FLOAT32_MAX_EXPONENT ) || ((*exp) < 0)) {
+	}
+	
+	if (((*exp) >= FLOAT32_MAX_EXPONENT) || ((*exp) < 0)) {
 		/* overflow - set infinity as result */
 		(*exp) = FLOAT32_MAX_EXPONENT;
 		(*fraction) = 0;
-		return;
-	}
-
-	return;
-}
-
-/** Round and normalize number expressed by exponent and fraction with first bit (equal to hidden bit) at 62. bit
- * @param exp exponent 
- * @param fraction part with hidden bit shifted to 62. bit
+	}
+}
+
+/**
+ * Round and normalize number expressed by exponent and fraction with
+ * first bit (equal to hidden bit) at 62nd bit.
+ *
+ * @param exp Exponent part.
+ * @param fraction Fraction with hidden bit shifted to 62nd bit.
  */
 void roundFloat64(int32_t *exp, uint64_t *fraction)
 {
 	/* rounding - if first bit after fraction is set then round up */
-	(*fraction) += (0x1 << 9);
-	
-	if ((*fraction) & (FLOAT64_HIDDEN_BIT_MASK << 11)) { 
+	(*fraction) += (0x1 << (64 - FLOAT64_FRACTION_SIZE - 3));
+	
+	if ((*fraction) & 
+	    (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 3))) {
 		/* rounding overflow */
 		++(*exp);
 		(*fraction) >>= 1;
-	};
-	
-	if (((*exp) >= FLOAT64_MAX_EXPONENT ) || ((*exp) < 0)) {
+	}
+	
+	if (((*exp) >= FLOAT64_MAX_EXPONENT) || ((*exp) < 0)) {
 		/* overflow - set infinity as result */
 		(*exp) = FLOAT64_MAX_EXPONENT;
 		(*fraction) = 0;
-		return;
-	}
-
-	return;
+	}
+}
+
+/**
+ * Round and normalize number expressed by exponent and fraction with
+ * first bit (equal to hidden bit) at 126th bit.
+ *
+ * @param exp Exponent part.
+ * @param frac_hi High part of fraction part with hidden bit shifted to 126th bit.
+ * @param frac_lo Low part of fraction part with hidden bit shifted to 126th bit.
+ */
+void roundFloat128(int32_t *exp, uint64_t *frac_hi, uint64_t *frac_lo)
+{
+	uint64_t tmp_hi, tmp_lo;
+
+	/* rounding - if first bit after fraction is set then round up */
+	lshift128(0x0ll, 0x1ll, (128 - FLOAT128_FRACTION_SIZE - 3), &tmp_hi, &tmp_lo);
+	add128(*frac_hi, *frac_lo, tmp_hi, tmp_lo, frac_hi, frac_lo);
+
+	lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
+	    (128 - FLOAT128_FRACTION_SIZE - 3), &tmp_hi, &tmp_lo);
+	and128(*frac_hi, *frac_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
+	if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
+		/* rounding overflow */
+		++(*exp);
+		rshift128(*frac_hi, *frac_lo, 1, frac_hi, frac_lo);
+	}
+
+	if (((*exp) >= FLOAT128_MAX_EXPONENT) || ((*exp) < 0)) {
+		/* overflow - set infinity as result */
+		(*exp) = FLOAT128_MAX_EXPONENT;
+		(*frac_hi) = 0;
+		(*frac_lo) = 0;
+	}
+}
+
+/**
+ * Logical shift left on the 128-bit operand.
+ *
+ * @param a_hi High part of the input operand.
+ * @param a_lo Low part of the input operand.
+ * @param shift Number of bits by witch to shift.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void lshift128(
+    uint64_t a_hi, uint64_t a_lo, int shift,
+    uint64_t *r_hi, uint64_t *r_lo)
+{
+	if (shift <= 0) {
+		/* do nothing */
+	} else if (shift >= 128) {
+		a_hi = 0;
+		a_lo = 0;
+	} else if (shift >= 64) {
+		a_hi = a_lo << (shift - 64);
+		a_lo = 0;
+	} else {
+		a_hi <<= shift;
+		a_hi |= a_lo >> (64 - shift);
+		a_lo <<= shift;
+	}
+
+	*r_hi = a_hi;
+	*r_lo = a_lo;
+}
+
+/**
+ * Logical shift right on the 128-bit operand.
+ *
+ * @param a_hi High part of the input operand.
+ * @param a_lo Low part of the input operand.
+ * @param shift Number of bits by witch to shift.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void rshift128(
+    uint64_t a_hi, uint64_t a_lo, int shift,
+    uint64_t *r_hi, uint64_t *r_lo)
+{
+	if (shift <= 0) {
+		/* do nothing */
+	} else 	if (shift >= 128) {
+		a_hi = 0;
+		a_lo = 0;
+	} else if (shift >= 64) {
+		a_lo = a_hi >> (shift - 64);
+		a_hi = 0;
+	} else {
+		a_lo >>= shift;
+		a_lo |= a_hi << (64 - shift);
+		a_hi >>= shift;
+	}
+
+	*r_hi = a_hi;
+	*r_lo = a_lo;
+}
+
+/**
+ * Bitwise AND on 128-bit operands.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void and128(
+    uint64_t a_hi, uint64_t a_lo,
+    uint64_t b_hi, uint64_t b_lo,
+    uint64_t *r_hi, uint64_t *r_lo)
+{
+	*r_hi = a_hi & b_hi;
+	*r_lo = a_lo & b_lo;
+}
+
+/**
+ * Bitwise inclusive OR on 128-bit operands.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void or128(
+    uint64_t a_hi, uint64_t a_lo,
+    uint64_t b_hi, uint64_t b_lo,
+    uint64_t *r_hi, uint64_t *r_lo)
+{
+	*r_hi = a_hi | b_hi;
+	*r_lo = a_lo | b_lo;
+}
+
+/**
+ * Bitwise exclusive OR on 128-bit operands.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void xor128(
+    uint64_t a_hi, uint64_t a_lo,
+    uint64_t b_hi, uint64_t b_lo,
+    uint64_t *r_hi, uint64_t *r_lo)
+{
+	*r_hi = a_hi ^ b_hi;
+	*r_lo = a_lo ^ b_lo;
+}
+
+/**
+ * Bitwise NOT on the 128-bit operand.
+ *
+ * @param a_hi High part of the input operand.
+ * @param a_lo Low part of the input operand.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void not128(
+    uint64_t a_hi, uint64_t a_lo,
+	uint64_t *r_hi, uint64_t *r_lo)
+{
+	*r_hi = ~a_hi;
+	*r_lo = ~a_lo;
+}
+
+/**
+ * Equality comparison of 128-bit operands.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @return 1 if operands are equal, 0 otherwise.
+ */
+int eq128(uint64_t a_hi, uint64_t a_lo, uint64_t b_hi, uint64_t b_lo)
+{
+	return (a_hi == b_hi) && (a_lo == b_lo);
+}
+
+/**
+ * Lower-or-equal comparison of 128-bit operands.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @return 1 if a is lower or equal to b, 0 otherwise.
+ */
+int le128(uint64_t a_hi, uint64_t a_lo, uint64_t b_hi, uint64_t b_lo)
+{
+	return (a_hi < b_hi) || ((a_hi == b_hi) && (a_lo <= b_lo));
+}
+
+/**
+ * Lower-than comparison of 128-bit operands.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @return 1 if a is lower than b, 0 otherwise.
+ */
+int lt128(uint64_t a_hi, uint64_t a_lo, uint64_t b_hi, uint64_t b_lo)
+{
+	return (a_hi < b_hi) || ((a_hi == b_hi) && (a_lo < b_lo));
+}
+
+/**
+ * Addition of two 128-bit unsigned integers.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void add128(uint64_t a_hi, uint64_t a_lo,
+    uint64_t b_hi, uint64_t b_lo,
+    uint64_t *r_hi, uint64_t *r_lo)
+{
+	uint64_t low = a_lo + b_lo;
+	*r_lo = low;
+	/* detect overflow to add a carry */
+	*r_hi = a_hi + b_hi + (low < a_lo);
+}
+
+/**
+ * Substraction of two 128-bit unsigned integers.
+ *
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void sub128(uint64_t a_hi, uint64_t a_lo,
+    uint64_t b_hi, uint64_t b_lo,
+    uint64_t *r_hi, uint64_t *r_lo)
+{
+	*r_lo = a_lo - b_lo;
+	/* detect underflow to substract a carry */
+	*r_hi = a_hi - b_hi - (a_lo < b_lo);
+}
+
+/**
+ * Multiplication of two 64-bit unsigned integers.
+ * 
+ * @param a First input operand.
+ * @param b Second input operand.
+ * @param r_hi Address to store high part of the result.
+ * @param r_lo Address to store low part of the result.
+ */
+void mul64(uint64_t a, uint64_t b, uint64_t *r_hi, uint64_t *r_lo)
+{
+	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);
+	*r_lo = low;
+	*r_hi = high;
+}
+
+/**
+ * Multiplication of two 128-bit unsigned integers.
+ * 
+ * @param a_hi High part of the first input operand.
+ * @param a_lo Low part of the first input operand.
+ * @param b_hi High part of the second input operand.
+ * @param b_lo Low part of the second input operand.
+ * @param r_hihi Address to store first (highest) quarter of the result.
+ * @param r_hilo Address to store second quarter of the result.
+ * @param r_lohi Address to store third quarter of the result.
+ * @param r_lolo Address to store fourth (lowest) quarter of the result.
+ */
+void mul128(uint64_t a_hi, uint64_t a_lo, uint64_t b_hi, uint64_t b_lo,
+    uint64_t *r_hihi, uint64_t *r_hilo, uint64_t *r_lohi, uint64_t *r_lolo)
+{
+	uint64_t hihi, hilo, lohi, lolo;
+	uint64_t tmp1, tmp2;
+
+	mul64(a_lo, b_lo, &lohi, &lolo);
+	mul64(a_lo, b_hi, &hilo, &tmp2);
+	add128(hilo, tmp2, 0x0ll, lohi, &hilo, &lohi);
+	mul64(a_hi, b_hi, &hihi, &tmp1);
+	add128(hihi, tmp1, 0x0ll, hilo, &hihi, &hilo);
+	mul64(a_hi, b_lo, &tmp1, &tmp2);
+	add128(tmp1, tmp2, 0x0ll, lohi, &tmp1, &lohi);
+	add128(hihi, hilo, 0x0ll, tmp1, &hihi, &hilo);
+
+	*r_hihi = hihi;
+	*r_hilo = hilo;
+	*r_lohi = lohi;
+	*r_lolo = lolo;
+}
+
+/**
+ * Estimate the quotient of 128-bit unsigned divident and 64-bit unsigned
+ * divisor.
+ * 
+ * @param a_hi High part of the divident.
+ * @param a_lo Low part of the divident.
+ * @param b Divisor.
+ * @return Quotient approximation.
+ */
+uint64_t div128est(uint64_t a_hi, uint64_t a_lo, uint64_t b)
+{
+	uint64_t b_hi, b_lo;
+	uint64_t rem_hi, rem_lo;
+	uint64_t tmp_hi, tmp_lo;
+	uint64_t result;
+
+	if (b <= a_hi) {
+		return 0xFFFFFFFFFFFFFFFFull;
+	}
+
+	b_hi = b >> 32;
+	result = ((b_hi << 32) <= a_hi) ? (0xFFFFFFFFull << 32) : (a_hi / b_hi) << 32;
+	mul64(b, result, &tmp_hi, &tmp_lo);
+	sub128(a_hi, a_lo, tmp_hi, tmp_lo, &rem_hi, &rem_lo);
+	
+	while ((int64_t) rem_hi < 0) {
+		result -= 0x1ll << 32;
+		b_lo = b << 32;
+		add128(rem_hi, rem_lo, b_hi, b_lo, &rem_hi, &rem_lo);
+	}
+
+	rem_hi = (rem_hi << 32) | (rem_lo >> 32);
+	if ((b_hi << 32) <= rem_hi) {
+		result |= 0xFFFFFFFF;
+	} else {
+		result |= rem_hi / b_hi;
+	}
+
+	return result;
 }
 
