Hybrid GCD core routine that combines Stein and Euclid

hybrid_gcd.c

// result of internal core
typedef struct {
  uint64_t r;   // common odd factors
  uint32_t s;   // common powers-of-two (as shift amount)
} gcd_t;

// core GCD routines: for 'm' significant bit & uniform inputs.
// Probablity comments in each routine also assume the above.
// 
// Expected number of iterations per method:
// ∙ math modeled (asymptotic: m → ∞)
//     Euclid: 0.5841608166566490218792*m
//     Stein:  0.7059712461019163915293*m
//     Hybrid: N/A
// ∙ empirical data (m on [8,64] & uniform)
//     Euclid: 0.5786530030011198*m
//     Stein:  0.7076493159462484*m
//     Hybrid: 0.3697942106047032*m
// ∙ empirical data (better linear model for m ≥ 8)
//     Euclid: 0.5841647345086854*m - 0.23986238968109566
//     Stein:  0.7059484193673851*m + 0.07402049926536608
//     Hybrid: 0.3667655628078821*m + 0.13180226523203306
// On entry the probablity of the inputs being coprime:
// ∙ 1/Zeta(2) = 6/π² ≈ 0.607927 (for unbounded integers)
// ∙ empirical data agrees for m ≥ 8
// I think it'd be better that the inputs were a tuple
// but I'm sticking with expected style for a C/C++ 

// Hybrid (Euclid/Stein) for hardware with fast integer divides
// ∙ like Stein starts by removing all common powers of two
// ∙ each iteration perform a Euclid then a Stein step
gcd_t gcd_core(uint64_t a, uint64_t b)
{
  uint64_t t;

  // perform the Stein remove all common powers-of-two
  uint32_t s = ctz_64(a|b);
  a >>= s;
  b >>= s;
  
  hint_cswap(b>a,a,b);  // swap on condition (macro with compiler hint to be brachless. Doesn't matter much out of loop)

  // the probablity the remaining terms are coprime
  // is 8/π²  ≈  0.810569
  while (b > 0) {
    // Euclid step
    a %= b;
    t = a; a = b; b = t;
    // Stein step
    b >>= ctz_64(b);
  }
  return (gcd_t){.r=a, .s=s};
}