gcd.h

// result for blog post on my machine
// 10.7 Ma/s  93.60 ns/d   4.17 GHz  390.21 c/d  778.90 i/d  24.39 c/b  48.68 i/b   2.00 i/c 
//                                                                                  ^^^^^^^^
//
// binary_gcd_paolo_1 wins for #instructions, but is really bad in terms of ipc
// 5.3 Ma/s  189.59 ns/d   4.16 GHz  788.92 c/d  324.72 i/d  49.31 c/b  20.30 i/b   0.41 i/c 
//                                               ^^^^^^^^^^
//
// binary_gcd_paolo_2 wins overall on my machine, saving more in #instructions than it loses ipc
// 12.0 Ma/s  82.99 ns/d   4.15 GHz  344.06 c/d  592.35 i/d  21.50 c/b  37.02 i/b   1.72 i/c 
// ^^^^^^^^^

// avoid swap altogether, at the cost of more mispredictions
template <typename int_type>
int_type binary_gcd_paolo_1(int_type u, int_type v) {
  if (u == 0) { return v; }
  if (v == 0) { return u; }
  auto shift = std::countr_zero(u | v);
  u >>= std::countr_zero(u);
  v >>= std::countr_zero(v);
  if (u <= v) goto le;
  for(;;) {
    do {
      u -= v;
      u >>= std::countr_zero(u);
    } while (u > v);
le:
    if (__builtin_expect(v == u, 0)) break;
    do {
      v -= u;
      v >>= std::countr_zero(v);
    } while (v > u);
    if (__builtin_expect(v == u, 0)) break;
  }
  return u << shift;
}

// just inline std::swap
template <typename int_type>
int_type binary_gcd_paolo_2(int_type u, int_type v) {
  if (u == 0) {
    return v;
  }
  if (v == 0) {
    return u;
  }
  auto shift = std::countr_zero(u | v);
  u >>= std::countr_zero(u);
  do {
    int_type t = v >> std::countr_zero(v);
    if (u > t)
      v = u - t, u = t;
    else
      v = t - u;
  } while (v != 0);
  return u << shift;
}