Changeset c67aff2 in mainline for uspace/lib/softfloat/generic/common.c

Timestamp:

2011-08-06T07:04:50Z (13 years ago)

Author:

Petr Koupy <petr.koupy@…>

Branches:

lfn, master, serial, ticket/834-toolchain-update, topic/msim-upgrade, topic/simplify-dev-export

Children:

d3e241a, e0e922d

Parents:

9a6034a

Message:

Quadruple-precision softfloat, coding style improvements. Details below…

Highlights:

• completed double-precision support
• added quadruple-precision support
• added SPARC quadruple-precision wrappers
• added doxygen comments
• corrected and unified coding style

Current state of the softfloat library:

Support for single, double and quadruple precision is currently almost complete (apart from power, square root, complex multiplication and complex division) and provides the same set of features (i.e. the support for all three precisions is now aligned). In order to extend softfloat library consistently, addition of quadruple precision was done in the same spirit as already existing single and double precision written by Josef Cejka in 2006 - that is relaxed standard-compliance for corner cases while mission-critical code sections heavily inspired by the widely used softfloat library written by John R. Hauser (although I personally think it would be more appropriate for HelenOS to use something less optimized, shorter and more readable).

Most of the quadruple-precision code is just an adapted double-precision code to work on 128-bit variables. That means if there is TODO, FIXME or some defect in single or double-precision code, it is most likely also in the quadruple-precision code. Please note that quadruple-precision functions are currently not tested - it is challenging task for itself, especially when the ports that use them are either not finished (mips64) or badly supported by simulators (sparc64). To test whole softfloat library, one would probably have to either write very non-trivial native tester, or use some existing one (e.g. TestFloat from J. R. Hauser) and port it to HelenOS (or rip the softfloat library out of HelenOS and test it on a host system). At the time of writing this, the code dependent on quadruple-precision functions (on mips64 and sparc64) is just a libposix strtold() function (and its callers, most notably scanf backend).

File:

• uspace/lib/softfloat/generic/common.c (modified) (10 diffs)

uspace/lib/softfloat/generic/common.c

@@ -1 +1 @@
 /*
  * Copyright (c) 2005 Josef Cejka
+ * Copyright (c) 2011 Petr Koupy
  * All rights reserved.
  *
@@ -30 +31 @@
  * @{
  */
-/** @file
+/** @file Common helper operations.
  */
 
@@ -36 +37 @@
 #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 +57 @@
 };
 
-
-
-/** 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 +73 @@
 
         /* 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 +97 @@
                 
                 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 +107 @@
         ++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 +113 @@
 
         /* check overflow */
-        if (cexp >= FLOAT64_MAX_EXPONENT ) {
+        if (cexp >= FLOAT64_MAX_EXPONENT) {
                 /* FIXME: overflow, return infinity */
                 result.parts.exp = FLOAT64_MAX_EXPONENT;
@@ -116 +120 @@
         }
 
-        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 +279 @@
 }
 
-/** 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;
 }