How can I use mach_absolute_time without overflowing?

Question

On Darwin, the POSIX standard clock_gettime(CLOCK_MONOTONIC) timer is not available. Instead, the highest resolution monotonic timer is obtained through the mach_absolute_time function from mach/mach_time.h.

The result returned may be an unadjusted tick count from the processor, in which case the time units could be a strange multiple. For example, on a CPU with a 33MHz tick count, Darwin returns 1000000000/33333335 as the exact units of the returned result (ie, multiply the mach_absolute_time by that fraction to obtain a nanosecond value).

We usually wish to convert from exact ticks to "standard" (decimal) units, but unfortunately, naively multiplying the absolute time by the fraction will overflow even in 64-bit arithmetic. This is an error that Apple's sole piece of documentation on mach_absolute_time falls into (Technical Q&A QA1398).¹

How should I write a function that correctly uses mach_absolute_time?

---

1. Note that this is not a theoretical problem: the sample code in QA1398 completely fails to work on PowerPC-based Macs. On Intel Macs, mach_timebase_info always returns 1/1 as the scaling factor because the CPU's raw tick count is unreliable (dynamic speed-stepping), so the API does the scaling for you. On PowerPC Macs, mach_timebase_info returns either 1000000000/33333335 or 1000000000/25000000, so Apple's provided code definitely overflows every few minutes. Oops.

Answer 1

Most-precise (best) answer

Perform the arithmetic at 128-bit precision to avoid the overflow!

// Returns monotonic time in nanos, measured from the first time the function
// is called in the process.
uint64_t monotonicTimeNanos() {
  uint64_t now = mach_absolute_time();
  static struct Data {
    Data(uint64_t bias_) : bias(bias_) {
      kern_return_t mtiStatus = mach_timebase_info(&tb);
      assert(mtiStatus == KERN_SUCCESS);
    }
    uint64_t scale(uint64_t i) {
      return scaleHighPrecision(i - bias, tb.numer, tb.denom);
    }
    static uint64_t scaleHighPrecision(uint64_t i, uint32_t numer,
                                       uint32_t denom) {
      U64 high = (i >> 32) * numer;
      U64 low = (i & 0xffffffffull) * numer / denom;
      U64 highRem = ((high % denom) << 32) / denom;
      high /= denom;
      return (high << 32) + highRem + low;
    }
    mach_timebase_info_data_t tb;
    uint64_t bias;
  } data(now);
  return data.scale(now);
}

A simple low-resolution answer

// Returns monotonic time in nanos, measured from the first time the function
// is called in the process.  The clock may run up to 0.1% faster or slower
// than the "exact" tick count.
uint64_t monotonicTimeNanos() {
  uint64_t now = mach_absolute_time();
  static struct Data {
    Data(uint64_t bias_) : bias(bias_) {
      kern_return_t mtiStatus = mach_timebase_info(&tb);
      assert(mtiStatus == KERN_SUCCESS);
      if (tb.denom > 1024) {
        double frac = (double)tb.numer/tb.denom;
        tb.denom = 1024;
        tb.numer = tb.denom * frac + 0.5;
        assert(tb.numer > 0);
      }
    }
    mach_timebase_info_data_t tb;
    uint64_t bias;
  } data(now);
  return (now - data.bias) * data.tb.numer / data.tb.denom;
}

A fiddly solution using low-precision arithmetic but using continued fractions to avoid loss of accuracy

// This function returns the rational number inside the given interval with
// the smallest denominator (and smallest numerator breaks ties; correctness
// proof neglects floating-point errors).
static mach_timebase_info_data_t bestFrac(double a, double b) {
  if (floor(a) < floor(b))
  { mach_timebase_info_data_t rv = {(int)ceil(a), 1}; return rv; }
  double m = floor(a);
  mach_timebase_info_data_t next = bestFrac(1/(b-m), 1/(a-m));
  mach_timebase_info_data_t rv = {(int)m*next.numer + next.denum, next.numer};
  return rv;
}

// Returns monotonic time in nanos, measured from the first time the function
// is called in the process.  The clock may run up to 0.1% faster or slower
// than the "exact" tick count. However, although the bound on the error is
// the same as for the pragmatic answer, the error is actually minimized over
// the given accuracy bound.
uint64_t monotonicTimeNanos() {
  uint64_t now = mach_absolute_time();
  static struct Data {
    Data(uint64_t bias_) : bias(bias_) {
      kern_return_t mtiStatus = mach_timebase_info(&tb);
      assert(mtiStatus == KERN_SUCCESS);
      double frac = (double)tb.numer/tb.denom;
      uint64_t spanTarget = 315360000000000000llu; // 10 years
      if (getExpressibleSpan(tb.numer, tb.denom) >= spanTarget)
        return;
      for (double errorTarget = 1/1024.0; errorTarget > 0.000001;) {
        mach_timebase_info_data_t newFrac =
            bestFrac((1-errorTarget)*frac, (1+errorTarget)*frac);
        if (getExpressibleSpan(newFrac.numer, newFrac.denom) < spanTarget)
          break;
        tb = newFrac;
        errorTarget = fabs((double)tb.numer/tb.denom - frac) / frac / 8;
      }
      assert(getExpressibleSpan(tb.numer, tb.denom) >= spanTarget);
    }
    mach_timebase_info_data_t tb;
    uint64_t bias;
  } data(now);
  return (now - data.bias) * data.tb.numer / data.tb.denom;
}

The derivation

We aim to reduce the fraction returned by mach_timebase_info to one that is essentially the same, but with a small denominator. The size of the timespan that we can handle is limited only by the size of the denominator, not the numerator of the fraction we shall multiply by:

uint64_t getExpressibleSpan(uint32_t numer, uint32_t denom) {
  // This is just less than the smallest thing we can multiply numer by without
  // overflowing. ceilLog2(numer) = 64 - number of leading zeros of numer
  uint64_t maxDiffWithoutOverflow = ((uint64_t)1 << (64 - ceilLog2(numer))) - 1;
  return maxDiffWithoutOverflow * numer / denom;
}

If denom=33333335 as returned by mach_timebase_info, we can handle differences of up to 18 seconds only before the multiplication by numer overflows. As getExpressibleSpan shows, by calculating a rough lower bound for this, the size of numer doesn't matter: halving numer doubles maxDiffWithoutOverflow. The only goal therefore is to produce a fraction close to numer/denom that has a smaller denominator. The simplest method to do this is using continued fractions.

The continued fractions method is rather handy. bestFrac clearly works correctly if the provided interval contains an integer: it returns the least integer in the interval over 1. Otherwise, it calls itself recursively with a strictly larger interval and returns m+1/next. The final result is a continued fraction that can be shown by induction to have the correct property: it's optimal, the fraction inside the given interval with the least denominator.

Finally, we reduce the fraction Darwin passes us to a smaller one to use when rescaling the mach_absolute_time to nanoseconds. We may introduce an error here because we can't reduce the fraction in general without losing accuracy. We set ourselves the target of 0.1% error, and check that we've reduced the fraction enough for common timespans (up to ten years) to be handled correctly.

Arguably the method is over-complicated for what it does, but it handles correctly anything the API can throw at it, and the resulting code is still short and extremely fast (bestFrac typically recurses only three or four iterations deep before returning a denominator less than 1000 for random intervals [a,a*1.002]).