test harness for comparing __uint128_t "gcd()" against gmplib "mpz_mod()"

gcd.c

#include <stdio.h>
#include <immintrin.h>
#include <time.h>
#include <sys/time.h>

#include <stdlib.h>
#include <string.h>
#include <inttypes.h>

#include "gmp.h"

#define N 1000000

void P(__int128_t i)
{
  printf("%016lx %016lx\n", (long) (i>>64), (long) (i & UINT64_MAX));
}

inline int ctz(__uint128_t u)
{
  unsigned long long h = u;
  return (h!=0) ?      __builtin_ctzll( h )
                : 64 + __builtin_ctzll( u>>64 );
}

unsigned long long _gcd(unsigned long long u, unsigned long long v)
{
  for(;;) {
    if (u > v) { unsigned long long a=u; u=v; v=a; }

    v -= u;

    if (v == 0)  return u;

    v >>= __builtin_ctzll(v);
  }
}

__uint128_t gcd(__uint128_t u, __uint128_t v)
{
       if (u == 0) { return v; }
  else if (v == 0) { return u; }

  int i = ctz(u);  u >>= i;
  int j = ctz(v);  v >>= j;
  int k = (i < j) ? i : j;

  for(;;) {
    if (u > v) { __uint128_t a=u; u=v; v=a; }

    if (v <= UINT64_MAX)  return _gcd(u, v) << k;

    v -= u;

    if (v == 0)  return u << k;

    v >>= ctz(v);
  }
}

int main() {
  struct timespec t0, t1;
  __int128_t s;
  __uint128_t a = 0x9abfd87547c0e48cL;
  a <<= 64;
  a |= 0x30173357e778cd8dL;
  __uint128_t b = 0xfa63c8d9fa216a8fL;
  b <<= 64;
  b |= 0xc8a7213b333270f8L;

  s = 0;
  clock_gettime(CLOCK_REALTIME, &t0);
  for(int i=0; i<N; ++i)
    s += gcd(a, b);
  clock_gettime(CLOCK_REALTIME, &t1);
  printf("%ldns(1)\n", (t1.tv_sec-t0.tv_sec) * 1000000000 + t1.tv_nsec - t0.tv_nsec);
  P(s);

  {
    mpz_t a,b,s,t;

    mpz_init(t);

    mpz_init_set_ui(s, 0);

    mpz_init_set_ui(a, 0x9abfd87547c0e48cL);
    mpz_mul_2exp(t, a, 64);
    mpz_add_ui(a, t, 0x30173357e778cd8dL);

    mpz_init_set_ui(b, 0xfa63c8d9fa216a8fL);
    mpz_mul_2exp(t, b, 64);
    mpz_add_ui(b, t, 0xc8a7213b333270f8L);

    clock_gettime(CLOCK_REALTIME, &t0);
    for(int i=0; i<N; ++i)
    {
      mpz_gcd(t, a, b);
      mpz_add(s, s, t);
    }
    clock_gettime(CLOCK_REALTIME, &t1);
    printf("%ldns(2)\n", (t1.tv_sec-t0.tv_sec) * 1000000000 + t1.tv_nsec - t0.tv_nsec);
    gmp_printf ("%Zx\n", s);
  } 

  exit (0);
}