gstutils: Refactor gst_util_uint64_scale()
authorKipp Cannon <kcannon@ligo.caltech.edu>
Tue, 11 Aug 2009 07:10:47 +0000 (09:10 +0200)
committerSebastian Dröge <sebastian.droege@collabora.co.uk>
Thu, 13 Aug 2009 14:32:27 +0000 (16:32 +0200)
This will later make it possible to provide rounding versions
of it without much code duplication.

Partially fixes bug #590919.

gst/gstutils.c

index f9155e1..302e3f8 100644 (file)
@@ -204,83 +204,19 @@ typedef union
   } l;
 } GstUInt64;
 
-/* based on Hacker's Delight p152 */
-static guint64
-gst_util_div128_64 (GstUInt64 c1, GstUInt64 c0, guint64 denom)
-{
-  GstUInt64 q1, q0, rhat;
-  GstUInt64 v, cmp1, cmp2;
-  guint s;
-
-  v.ll = denom;
-
-  /* count number of leading zeroes, we know they must be in the high
-   * part of denom since denom > G_MAXUINT32. */
-  s = v.l.high | (v.l.high >> 1);
-  s |= (s >> 2);
-  s |= (s >> 4);
-  s |= (s >> 8);
-  s = ~(s | (s >> 16));
-  s = s - ((s >> 1) & 0x55555555);
-  s = (s & 0x33333333) + ((s >> 2) & 0x33333333);
-  s = (s + (s >> 4)) & 0x0f0f0f0f;
-  s += (s >> 8);
-  s = (s + (s >> 16)) & 0x3f;
-
-  if (s > 0) {
-    /* normalize divisor and dividend */
-    v.ll <<= s;
-    c1.ll = (c1.ll << s) | (c0.l.high >> (32 - s));
-    c0.ll <<= s;
-  }
-
-  q1.ll = c1.ll / v.l.high;
-  rhat.ll = c1.ll - q1.ll * v.l.high;
-
-  cmp1.l.high = rhat.l.low;
-  cmp1.l.low = c0.l.high;
-  cmp2.ll = q1.ll * v.l.low;
-
-  while (q1.l.high || cmp2.ll > cmp1.ll) {
-    q1.ll--;
-    rhat.ll += v.l.high;
-    if (rhat.l.high)
-      break;
-    cmp1.l.high = rhat.l.low;
-    cmp2.ll -= v.l.low;
-  }
-  c1.l.high = c1.l.low;
-  c1.l.low = c0.l.high;
-  c1.ll -= q1.ll * v.ll;
-  q0.ll = c1.ll / v.l.high;
-  rhat.ll = c1.ll - q0.ll * v.l.high;
-
-  cmp1.l.high = rhat.l.low;
-  cmp1.l.low = c0.l.low;
-  cmp2.ll = q0.ll * v.l.low;
-
-  while (q0.l.high || cmp2.ll > cmp1.ll) {
-    q0.ll--;
-    rhat.ll += v.l.high;
-    if (rhat.l.high)
-      break;
-    cmp1.l.high = rhat.l.low;
-    cmp2.ll -= v.l.low;
-  }
-  q0.l.high += q1.l.low;
-
-  return q0.ll;
-}
-
-static guint64
-gst_util_uint64_scale_int64_unchecked (guint64 val, guint64 num, guint64 denom)
+/* multiply two 64-bit unsigned ints into a 128-bit unsigned int.  the high
+ * and low 64 bits of the product are placed in c1 and c0 respectively.
+ * this operation cannot overflow. */
+static void
+gst_util_uint64_mul_uint64 (GstUInt64 * c1, GstUInt64 * c0, guint64 arg1,
+    guint64 arg2)
 {
-  GstUInt64 a0, a1, b0, b1, c0, ct, c1, result;
+  GstUInt64 a1, b0;
   GstUInt64 v, n;
 
   /* prepare input */
-  v.ll = val;
-  n.ll = num;
+  v.ll = arg1;
+  n.ll = arg2;
 
   /* do 128 bits multiply
    *                   nh   nl
@@ -291,79 +227,169 @@ gst_util_uint64_scale_int64_unchecked (guint64 val, guint64 num, guint64 denom)
    * b0 =         vh * nl
    * b1 =  + vh * nh
    *       -------------------
-   * c1,c0
+   *        c1h  c1l  c0h  c0l
+   *
+   * "a0" is optimized away, result is stored directly in c0.  "b1" is
+   * optimized away, result is stored directly in c1.
    */
-  a0.ll = (guint64) v.l.low * n.l.low;
+  c0->ll = (guint64) v.l.low * n.l.low;
   a1.ll = (guint64) v.l.low * n.l.high;
   b0.ll = (guint64) v.l.high * n.l.low;
-  b1.ll = (guint64) v.l.high * n.l.high;
-
-  /* and sum together with carry into 128 bits c1, c0 */
-  c0.l.low = a0.l.low;
-  ct.ll = (guint64) a0.l.high + a1.l.low + b0.l.low;
-  c0.l.high = ct.l.low;
-  c1.ll = (guint64) a1.l.high + b0.l.high + ct.l.high + b1.ll;
-
-  /* if high bits bigger than denom, we overflow */
-  if (G_UNLIKELY (c1.ll >= denom))
-    goto overflow;
-
-  /* shortcut for division by 1, c1.ll should be 0 because of the
-   * overflow check above. */
-  if (denom == 1)
-    return c0.ll;
-
-  /* and 128/64 bits division, result fits 64 bits */
-  if (denom <= G_MAXUINT32) {
-    guint32 den = (guint32) denom;
-
-    /* easy case, (c1,c0)128/(den)32 division */
-    c1.l.high %= den;
-    c1.l.high = c1.ll % den;
-    c1.l.low = c0.l.high;
-    c0.l.high = c1.ll % den;
-    result.l.high = c1.ll / den;
-    result.l.low = c0.ll / den;
-  } else {
-    result.ll = gst_util_div128_64 (c1, c0, denom);
+
+  /* add the high word of a0 to the low words of a1 and b0 using c1 as
+   * scrach space to capture the carry.  the low word of the result becomes
+   * the final high word of c0 */
+  c1->ll = (guint64) c0->l.high + a1.l.low + b0.l.low;
+  c0->l.high = c1->l.low;
+
+  /* add the carry from the result above (found in the high word of c1) and
+   * the high words of a1 and b0 to b1, the result is c1. */
+  c1->ll = (guint64) v.l.high * n.l.high + c1->l.high + a1.l.high + b0.l.high;
+}
+
+/* compute the quotient and remainder of 2^64 / d.  returns 0 if the
+ * quotient overflows (meaning d = 1). */
+static guint64
+gst_util_two_to_the_64_over_d (guint64 d, guint64 * remainder)
+{
+  guint64 quotient = G_MAXUINT64 / d;
+  *remainder = G_MAXUINT64 % d + 1;
+  if (*remainder == d) {
+    quotient++;
+    *remainder = 0;
   }
-  return result.ll;
+  return quotient;
+}
 
-overflow:
-  {
-    return G_MAXUINT64;
+/* divide a 128-bit unsigned int by a 64-bit unsigned int when we know the
+ * quotient fits into 64 bits. */
+static guint64
+gst_util_div128_64 (guint64 c1, guint64 c0, guint64 denom, guint64 * remainder)
+{
+  /* we are trying to compute
+   *
+   * c1 * 2^64 + c0
+   * --------------
+   *       d
+   *
+   * this can be re-written as:
+   *
+   * c1 * 2^64 + c0        2^64   c0
+   * -------------- = c1 * ---- + --
+   *       d                d      d
+   *
+   *                       (            2^64 % d )   c0
+   *                = c1 * (2^64 // d + ---------) + --
+   *                       (                d    )    d
+   *
+   *                                     c1 * (2^64 % d) + c0
+   *                = c1 * (2^64 // d) + --------------------
+   *                                              d
+   *
+   * where "//" indicates the integer quotient and "%" indicates remainder.
+   * note that 2^64 // d != 0 because d fits in 64 bits, and therefore if
+   * c1 != 0 the first term on the right-hand-side is != 0 and therefore
+   * the numerator in the fraction on the right-hand-side must be less than
+   * the numerator in the fraction on the left-hand-side.
+   *
+   * this provides us with an algorithm to compute both the quotient and
+   * remainder iteratively --- essentially a base-2^64 version of long
+   * division.  initializing the quotient to 0, the first term on the
+   * right-hand side is computed and added to the quotient (this can't
+   * overflow because we know the final answer fits in 64 bits).  the
+   * numerator of the second term is then computed and the high and low
+   * words stored in c1 and c0 respectively.  this is repeated until c1 is
+   * 0, at which point the problem has been reduced to computing the
+   * quotient and remainder of a 64-bit unsigned integer (c0) divided by a
+   * 64-bit unsigned integer (denom) which can be completed using regular
+   * integer arithmetic operations.
+   *
+   * note that gst_util_two_to_the_64_over_d() returns 0 if that quotient
+   * overflows.  this can only happen if d = 1, but because we know that
+   * our quotient must fit into 64 bits c1 must be 0 when d = 1, so the
+   * algorithm produces the correct result.
+   */
+
+  guint64 quotient = 0;
+
+  while (c1) {
+    guint64 a;
+    /* add c1 * (2^64 // d) to quotient, store 2^64 % d in a */
+    quotient += c1 * gst_util_two_to_the_64_over_d (denom, &a);
+    /* store the high and low words of c1 * (2^64 % d) in c1 and a
+     * respectively */
+    gst_util_uint64_mul_uint64 ((GstUInt64 *) & c1, (GstUInt64 *) & a, c1, a);
+    /* add a to c0, with a carry into c1 if the result rolls over */
+    if (G_MAXUINT64 - c0 < a)
+      c1++;
+    c0 += a;
   }
+
+  /* c1 is 0.  use regular integer arithmetic with c0 to complete result */
+  *remainder = c0 % denom;
+  return quotient + c0 / denom;
 }
 
-static inline guint64
-gst_util_uint64_scale_int_unchecked (guint64 val, gint num, gint denom)
+/* multiply a 64-bit unsigned int by a 32-bit unsigned int into a 96-bit
+ * unsigned int.  the high 64 bits and low 32 bits of the product are
+ * placed in c1 and c0 respectively.  this operation cannot overflow. */
+static void
+gst_util_uint64_mul_uint32 (GstUInt64 * c1, GstUInt64 * c0, guint64 arg1,
+    guint32 arg2)
 {
-  GstUInt64 result;
-  GstUInt64 low, high;
+  GstUInt64 a;
 
-  /* do 96 bits mult/div */
-  low.ll = val;
-  result.ll = ((guint64) low.l.low) * num;
-  high.ll = ((guint64) low.l.high) * num + (result.l.high);
+  a.ll = arg1;
 
-  low.ll = high.ll / denom;
-  result.l.high = high.ll % denom;
-  result.ll /= denom;
+  c0->ll = (guint64) a.l.low * arg2;
+  c1->ll = (guint64) a.l.high * arg2 + c0->l.high;
+  c0->l.high = 0;
+}
 
-  /* avoid overflow */
-  if (G_UNLIKELY (low.ll + result.l.high > G_MAXUINT32))
-    goto overflow;
+/* divide a 96-bit unsigned int by a 32-bit unsigned int when we know the
+ * quotient fits into 64 bits.  the high 64 bits and low 32 bits of the
+ * numerator are expected in c1 and c0 respectively. */
+static guint64
+gst_util_div96_32 (guint64 c1, guint64 c0, guint32 denom, guint32 * remainder)
+{
+  c0 += (c1 % denom) << 32;
+  *remainder = c0 % denom;
+  return ((c1 / denom) << 32) + (c0 / denom);
+}
 
-  result.l.high += low.l.low;
+static guint64
+gst_util_uint64_scale_uint64_unchecked (guint64 val, guint64 num,
+    guint64 denom, guint64 * remainder)
+{
+  guint64 c1, c0;
 
-  return result.ll;
+  /* compute 128-bit numerator product */
+  gst_util_uint64_mul_uint64 ((GstUInt64 *) & c1, (GstUInt64 *) & c0, val, num);
 
-overflow:
-  {
+  /* high word as big as or bigger than denom --> overflow */
+  if (G_UNLIKELY (c1 >= denom))
     return G_MAXUINT64;
-  }
+
+  /* compute quotient, fits in 64 bits */
+  return gst_util_div128_64 (c1, c0, denom, remainder);
 }
 
+static inline guint64
+gst_util_uint64_scale_uint32_unchecked (guint64 val, guint32 num,
+    guint32 denom, guint32 * remainder)
+{
+  GstUInt64 c1, c0;
+
+  /* compute 96-bit numerator product */
+  gst_util_uint64_mul_uint32 (&c1, &c0, val, num);
+
+  /* high 32 bits as big as or bigger than denom --> overflow */
+  if (G_UNLIKELY (c1.l.high >= denom))
+    return G_MAXUINT64;
+
+  /* compute quotient, fits in 64 bits */
+  return gst_util_div96_32 (c1.ll, c0.ll, denom, remainder);
+}
 
 /**
  * gst_util_uint64_scale:
@@ -371,16 +397,19 @@ overflow:
  * @num: the numerator of the scale ratio
  * @denom: the denominator of the scale ratio
  *
- * Scale @val by @num / @denom, trying to avoid overflows.
+ * Scale @val by the rational number @num / @denom, avoiding overflows and
+ * underflows and without loss of precision.
  *
- * This function can potentially be very slow if denom > G_MAXUINT32.
+ * This function can potentially be very slow if val and num are both
+ * greater than G_MAXUINT32.
  *
- * Returns: @val * @num / @denom, trying to avoid overflows.
- * In the case of an overflow, this function returns G_MAXUINT64.
+ * Returns: @val * @num / @denom.  In the case of an overflow, this
+ * function returns G_MAXUINT64.
  */
 guint64
 gst_util_uint64_scale (guint64 val, guint64 num, guint64 denom)
 {
+  guint64 remainder;
   g_return_val_if_fail (denom != 0, G_MAXUINT64);
 
   if (G_UNLIKELY (num == 0))
@@ -389,27 +418,22 @@ gst_util_uint64_scale (guint64 val, guint64 num, guint64 denom)
   if (G_UNLIKELY (num == denom))
     return val;
 
-  /* if the denom is high, we need to do a 64 muldiv */
-  if (G_UNLIKELY (denom > G_MAXINT32))
-    goto do_int64;
+  /* deneom is low --> try to use 96 bit muldiv */
+  if (G_LIKELY (denom <= G_MAXUINT32)) {
+    guint32 remainder;
+    /* num is low --> use 96 bit muldiv */
+    if (G_LIKELY (num <= G_MAXUINT32))
+      return gst_util_uint64_scale_uint32_unchecked (val, (guint32) num,
+          (guint32) denom, &remainder);
 
-  /* if num and denom are low we can do a 32 bit muldiv */
-  if (G_LIKELY (num <= G_MAXINT32))
-    goto do_int32;
-
-  /* val and num are high, we need 64 muldiv */
-  if (G_UNLIKELY (val > G_MAXINT32))
-    goto do_int64;
-
-  /* val is low and num is high, we can swap them and do 32 muldiv */
-  return gst_util_uint64_scale_int_unchecked (num, (gint) val, (gint) denom);
-
-do_int32:
-  return gst_util_uint64_scale_int_unchecked (val, (gint) num, (gint) denom);
+    /* num is high but val is low --> swap and use 96-bit muldiv */
+    if (G_LIKELY (val <= G_MAXUINT32))
+      return gst_util_uint64_scale_uint32_unchecked (num, (guint32) val,
+          (guint32) denom, &remainder);
+  }
 
-do_int64:
-  /* to the more heavy implementations... */
-  return gst_util_uint64_scale_int64_unchecked (val, num, denom);
+  /* val is high and num is high --> use 128-bit muldiv */
+  return gst_util_uint64_scale_uint64_unchecked (val, num, denom, &remainder);
 }
 
 /**
@@ -418,17 +442,17 @@ do_int64:
  * @num: numerator of the scale factor.
  * @denom: denominator of the scale factor.
  *
- * Scale a guint64 by a factor expressed as a fraction (num/denom), avoiding
- * overflows and loss of precision.
- *
- * @num and @denom must be positive integers. @denom cannot be 0.
+ * Scale @val by the rational number @num / @denom, avoiding overflows and
+ * underflows and without loss of precision.  @num must be non-negative and
+ * @denom must be positive.
  *
- * Returns: @val * @num / @denom, avoiding overflow and loss of precision.
- * In the case of an overflow, this function returns G_MAXUINT64.
+ * Returns: @val * @num / @denom.  In the case of an overflow, this
+ * function returns G_MAXUINT64.
  */
 guint64
 gst_util_uint64_scale_int (guint64 val, gint num, gint denom)
 {
+  guint32 remainder;
   g_return_val_if_fail (denom > 0, G_MAXUINT64);
   g_return_val_if_fail (num >= 0, G_MAXUINT64);
 
@@ -438,11 +462,16 @@ gst_util_uint64_scale_int (guint64 val, gint num, gint denom)
   if (G_UNLIKELY (num == denom))
     return val;
 
-  if (val <= G_MAXUINT32)
-    /* simple case */
-    return val * num / denom;
+  if (val <= G_MAXUINT32) {
+    /* simple case, use two statements to prevent optimizer from screwing
+     * up result.  num and denom are not negative so casts are OK */
+    val *= (guint64) num;
+    return val / (guint64) denom;
+  }
 
-  return gst_util_uint64_scale_int_unchecked (val, num, denom);
+  /* num and denom are not negative so casts are OK */
+  return gst_util_uint64_scale_uint32_unchecked (val, (guint32) num,
+      (guint32) denom, &remainder);
 }
 
 /**