C: Add option to compile using __int128 rather than GMP

Only requires a very small change to c_backend.ml. Most of this commit
is duplication of the builtins and runtime in lib/int128. But the
actual differences in those files is also fairly minor could be
handled by some simple ifdefs for the integer builtins.

1 files changed, 1558 insertions, 0 deletions

diff --git a/lib/int128/sail.c b/lib/int128/sail.c
new file mode 100644
index 00000000..6113c7f2
--- /dev/null
+++ b/lib/int128/sail.c
@@ -0,0 +1,1558 @@
+#define _GNU_SOURCE
+#include<assert.h>
+#include<inttypes.h>
+#include<stdbool.h>
+#include<stdio.h>
+#include<stdlib.h>
+#include<string.h>
+#include<time.h>
+
+#include <x86intrin.h>
+
+#include"sail.h"
+
+void mpz_set_si128(mpz_t rop, __int128 op)
+{
+  mpz_set_si(rop, (int64_t) (op >> (__int128) 64));
+  mpz_mul_2exp(rop, rop, 64);
+  mpz_add_ui(rop, rop, (uint64_t) op);
+}
+
+void mpz_init_set_si128(mpz_t rop, __int128 op)
+{
+  mpz_init(rop);
+  mpz_set_si128(rop, op);
+}
+
+bool EQUAL(unit)(const unit a, const unit b)
+{
+  return true;
+}
+
+unit UNDEFINED(unit)(const unit u)
+{
+  return UNIT;
+}
+
+unit skip(const unit u)
+{
+  return UNIT;
+}
+
+/* ***** Sail bit type ***** */
+
+bool eq_bit(const fbits a, const fbits b)
+{
+  return a == b;
+}
+
+/* ***** Sail booleans ***** */
+
+bool not(const bool b) {
+  return !b;
+}
+
+bool EQUAL(bool)(const bool a, const bool b) {
+  return a == b;
+}
+
+bool UNDEFINED(bool)(const unit u) {
+  return false;
+}
+
+/* ***** Sail strings ***** */
+
+void CREATE(sail_string)(sail_string *str)
+{
+  char *istr = (char *) malloc(1 * sizeof(char));
+  istr[0] = '\0';
+  *str = istr;
+}
+
+void RECREATE(sail_string)(sail_string *str)
+{
+  free(*str);
+  char *istr = (char *) malloc(1 * sizeof(char));
+  istr[0] = '\0';
+  *str = istr;
+}
+
+void COPY(sail_string)(sail_string *str1, const sail_string str2)
+{
+  size_t len = strlen(str2);
+  *str1 = realloc(*str1, len + 1);
+  *str1 = strcpy(*str1, str2);
+}
+
+void KILL(sail_string)(sail_string *str)
+{
+  free(*str);
+}
+
+void dec_str(sail_string *str, const sail_int n)
+{
+  if (INT64_MIN <= n && n <= INT64_MAX) {
+    int ret = asprintf(str, "%" PRId64, (int64_t) n);
+    if (ret == -1) {
+      printf("dec_str failed");
+      exit(1);
+    }
+  } else {
+    printf("dec_str");
+    exit(1);
+  }
+}
+
+void hex_str(sail_string *str, const sail_int n)
+{
+  //free(*str);
+  //gmp_asprintf(str, "0x%Zx", n);
+}
+
+bool eq_string(const sail_string str1, const sail_string str2)
+{
+  return strcmp(str1, str2) == 0;
+}
+
+bool EQUAL(sail_string)(const sail_string str1, const sail_string str2)
+{
+  return strcmp(str1, str2) == 0;
+}
+
+void undefined_string(sail_string *str, const unit u) {}
+
+void concat_str(sail_string *stro, const sail_string str1, const sail_string str2)
+{
+  *stro = realloc(*stro, strlen(str1) + strlen(str2) + 1);
+  (*stro)[0] = '\0';
+  strcat(*stro, str1);
+  strcat(*stro, str2);
+}
+
+bool string_startswith(sail_string s, sail_string prefix)
+{
+  return strstr(s, prefix) == s;
+}
+
+sail_int string_length(sail_string s)
+{
+  return (sail_int) strlen(s);
+}
+
+void string_drop(sail_string *dst, sail_string s, sail_int ns)
+{
+  size_t len = strlen(s);
+  mach_int n = CREATE_OF(mach_int, sail_int)(ns);
+  if (len >= n) {
+    *dst = realloc(*dst, (len - n) + 1);
+    memcpy(*dst, s + n, len - n);
+    (*dst)[len - n] = '\0';
+  } else {
+    *dst = realloc(*dst, 1);
+    **dst = '\0';
+  }
+}
+
+void string_take(sail_string *dst, sail_string s, sail_int ns)
+{
+  size_t len = strlen(s);
+  mach_int n = CREATE_OF(mach_int, sail_int)(ns);
+  mach_int to_copy;
+  if (len <= n) {
+    to_copy = len;
+  } else {
+    to_copy = n;
+  }
+  *dst = realloc(*dst, to_copy + 1);
+  memcpy(*dst, s, to_copy);
+  *dst[to_copy] = '\0';
+}
+
+/* ***** Sail integers ***** */
+
+uint64_t sail_int_get_ui(const sail_int op)
+{
+  return (uint64_t) op;
+}
+
+bool EQUAL(mach_int)(const mach_int op1, const mach_int op2)
+{
+  return op1 == op2;
+}
+
+sail_int CREATE_OF(sail_int, mach_int)(const mach_int op)
+{
+  return (sail_int) op;
+}
+
+mach_int CREATE_OF(mach_int, sail_int)(const sail_int op)
+{
+  return (mach_int) op;
+}
+
+mach_int CONVERT_OF(mach_int, sail_int)(const sail_int op)
+{
+  return (mach_int) op;
+}
+
+sail_int CONVERT_OF(sail_int, mach_int)(const mach_int op)
+{
+  return (sail_int) op;
+}
+
+sail_int CONVERT_OF(sail_int, sail_string)(const sail_string str)
+{
+  mpz_t tmp;
+  mpz_init(tmp);
+  mpz_set_str(tmp, str, 10);
+  uint64_t lo = mpz_get_ui(tmp);
+  mpz_div_2exp(tmp, tmp, 64);
+  uint64_t hi = mpz_get_ui(tmp);
+  mpz_clear(tmp);
+  
+  unsigned __int128 r = (((unsigned __int128) hi) << 64) + ((unsigned __int128) lo);
+  return (__int128) r;
+}
+
+bool eq_int(const sail_int op1, const sail_int op2)
+{
+  return op1 == op2;
+}
+
+bool EQUAL(sail_int)(const sail_int op1, const sail_int op2)
+{
+  return op1 == op2;
+}
+
+bool lt(const sail_int op1, const sail_int op2)
+{
+  return op1 < op2;
+}
+
+bool gt(const sail_int op1, const sail_int op2)
+{
+  return op1 > op2;
+}
+
+bool lteq(const sail_int op1, const sail_int op2)
+{
+  return op1 <= op2;
+}
+
+bool gteq(const sail_int op1, const sail_int op2)
+{
+  return op1 >= op2;
+}
+
+sail_int shl_int(const sail_int op1, const sail_int op2)
+{
+  return op1 << op2;
+}
+
+mach_int shl_mach_int(const mach_int op1, const mach_int op2)
+{
+  return op1 << op2;
+}
+
+sail_int shr_int(const sail_int op1, const sail_int op2)
+{
+  return op1 >> op2;
+}
+
+mach_int shr_mach_int(const mach_int op1, const mach_int op2)
+{
+  return op1 >> op2;
+}
+
+sail_int undefined_int(const int n)
+{
+  return (__int128) n;
+}
+
+sail_int undefined_range(const sail_int l, const sail_int u)
+{
+  return l;
+}
+
+sail_int add_int(const sail_int op1, const sail_int op2)
+{
+  return op1 + op2;
+}
+
+sail_int sub_int(const sail_int op1, const sail_int op2)
+{
+  return op1 - op2;
+}
+
+sail_int sub_nat(const sail_int op1, const sail_int op2)
+{
+  sail_int rop = op1 - op2;
+  if (rop < 0) return (sail_int) 0;
+  return rop;
+}
+
+sail_int mult_int(const sail_int op1, const sail_int op2)
+{
+  return op1 * op2;
+}
+
+// FIXME: Make sure all division operators do the right thing
+sail_int ediv_int(const sail_int op1, const sail_int op2)
+{
+  return op1 / op2;
+}
+
+sail_int emod_int(const sail_int op1, const sail_int op2)
+{
+  return op1 % op2;
+}
+
+sail_int tdiv_int(const sail_int op1, const sail_int op2)
+{
+  return op1 / op2;
+}
+
+sail_int tmod_int(const sail_int op1, const sail_int op2)
+{
+  return op1 % op2;
+}
+
+sail_int max_int(const sail_int op1, const sail_int op2)
+{
+  if (op1 < op2) {
+    return op2;
+  } else {
+    return op1;
+  }
+}
+
+sail_int min_int(const sail_int op1, const sail_int op2)
+{
+  if (op1 > op2) {
+    return op2;
+  } else {
+    return op1;
+  }
+}
+
+sail_int neg_int(const sail_int op)
+{
+  return -op;
+}
+
+sail_int abs_int(const sail_int op)
+{
+  if (op < 0) {
+    return -op;
+  } else {
+    return op;
+  }
+}
+
+sail_int pow_int(sail_int base, sail_int exp)
+{
+  sail_int result = 1;
+  while (true)
+  {
+    if (exp & 1) {
+      result *= base;
+    }
+    exp >>= 1;
+    if (!exp) {
+      break;
+    }
+    base *= base;
+  }
+  return result;
+}
+
+sail_int pow2(const sail_int exp)
+{
+  return pow_int(2, exp);
+}
+
+/* ***** Sail bitvectors ***** */
+
+bool EQUAL(fbits)(const fbits op1, const fbits op2)
+{
+  return op1 == op2;
+}
+
+void CREATE(lbits)(lbits *rop)
+{
+  rop->bits = malloc(sizeof(mpz_t));
+  rop->len = 0;
+  mpz_init(*rop->bits);
+}
+
+void RECREATE(lbits)(lbits *rop)
+{
+  rop->len = 0;
+  mpz_set_ui(*rop->bits, 0);
+}
+
+void COPY(lbits)(lbits *rop, const lbits op)
+{
+  rop->len = op.len;
+  mpz_set(*rop->bits, *op.bits);
+}
+
+void KILL(lbits)(lbits *rop)
+{
+  mpz_clear(*rop->bits);
+  free(rop->bits);
+}
+
+void CREATE_OF(lbits, fbits)(lbits *rop, const uint64_t op, const uint64_t len, const bool direction)
+{
+  rop->bits = malloc(sizeof(mpz_t));
+  rop->len = len;
+  mpz_init_set_ui(*rop->bits, op);
+}
+
+fbits CREATE_OF(fbits, lbits)(const lbits op, const bool direction)
+{
+  return mpz_get_ui(*op.bits);
+}
+
+sbits CREATE_OF(sbits, lbits)(const lbits op, const bool direction)
+{
+  sbits rop;
+  rop.bits = mpz_get_ui(*op.bits);
+  rop.len = op.len;
+  return rop;
+}
+
+sbits CREATE_OF(sbits, fbits)(const fbits op, const uint64_t len, const bool direction)
+{
+  sbits rop;
+  rop.bits = op;
+  rop.len = len;
+  return rop;
+}
+
+void RECREATE_OF(lbits, fbits)(lbits *rop, const uint64_t op, const uint64_t len, const bool direction)
+{
+  rop->len = len;
+  mpz_set_ui(*rop->bits, op);
+}
+
+void CREATE_OF(lbits, sbits)(lbits *rop, const sbits op, const bool direction)
+{
+  rop->bits = malloc(sizeof(mpz_t));
+  rop->len = op.len;
+  mpz_init_set_ui(*rop->bits, op.bits);
+}
+
+void RECREATE_OF(lbits, sbits)(lbits *rop, const sbits op, const bool direction)
+{
+  rop->len = op.len;
+  mpz_set_ui(*rop->bits, op.bits);
+}
+
+// Bitvector conversions
+
+fbits CONVERT_OF(fbits, lbits)(const lbits op, const bool direction)
+{
+  return mpz_get_ui(*op.bits);
+}
+
+fbits CONVERT_OF(fbits, sbits)(const sbits op, const bool direction)
+{
+  return op.bits;
+}
+
+void CONVERT_OF(lbits, fbits)(lbits *rop, const fbits op, const uint64_t len, const bool direction)
+{
+  rop->len = len;
+  // use safe_rshift to correctly handle the case when we have a 0-length vector.
+  mpz_set_ui(*rop->bits, op & safe_rshift(UINT64_MAX, 64 - len));
+}
+
+void CONVERT_OF(lbits, sbits)(lbits *rop, const sbits op, const bool direction)
+{
+  rop->len = op.len;
+  mpz_set_ui(*rop->bits, op.bits & safe_rshift(UINT64_MAX, 64 - op.len));
+}
+
+inline
+sbits CONVERT_OF(sbits, fbits)(const fbits op, const uint64_t len, const bool direction)
+{
+  sbits rop;
+  rop.len = len;
+  rop.bits = op;
+  return rop;
+}
+
+inline
+sbits CONVERT_OF(sbits, lbits)(const lbits op, const bool direction)
+{
+  sbits rop;
+  rop.len = op.len;
+  rop.bits = mpz_get_ui(*op.bits);
+  return rop;
+}
+
+void UNDEFINED(lbits)(lbits *rop, const sail_int len, const fbits bit)
+{
+  zeros(rop, len);
+}
+
+fbits UNDEFINED(fbits)(const unit u) { return 0; }
+
+sbits undefined_sbits(void)
+{
+  sbits rop;
+  rop.bits = UINT64_C(0);
+  rop.len = UINT64_C(0);
+  return rop;
+}
+
+fbits safe_rshift(const fbits x, const fbits n)
+{
+  if (n >= 64) {
+    return 0ul;
+  } else {
+    return x >> n;
+  }
+}
+
+void normalize_lbits(lbits *rop)
+{
+  mpz_t tmp;
+  mpz_init(tmp);
+  
+  mpz_set_ui(tmp, 1);
+  mpz_mul_2exp(tmp, tmp, rop->len);
+  mpz_sub_ui(tmp, tmp, 1);
+  mpz_and(*rop->bits, *rop->bits, tmp);
+
+  mpz_clear(tmp);
+}
+
+void append_64(lbits *rop, const lbits op, const fbits chunk)
+{
+  rop->len = rop->len + 64ul;
+  mpz_mul_2exp(*rop->bits, *op.bits, 64ul);
+  mpz_add_ui(*rop->bits, *rop->bits, chunk);
+}
+
+void add_bits(lbits *rop, const lbits op1, const lbits op2)
+{
+  rop->len = op1.len;
+  mpz_add(*rop->bits, *op1.bits, *op2.bits);
+  normalize_lbits(rop);
+}
+
+void sub_bits(lbits *rop, const lbits op1, const lbits op2)
+{
+  assert(op1.len == op2.len);
+  rop->len = op1.len;
+  mpz_sub(*rop->bits, *op1.bits, *op2.bits);
+  normalize_lbits(rop);
+}
+
+void add_bits_int(lbits *rop, const lbits op1, const sail_int op2)
+{
+  assert(op2 >= 0);
+  rop->len = op1.len;
+  mpz_add_ui(*rop->bits, *op1.bits, (uint64_t) op2);
+  normalize_lbits(rop);
+}
+
+void sub_bits_int(lbits *rop, const lbits op1, const sail_int op2)
+{
+  assert(op2 >= 0);
+  rop->len = op1.len;
+  mpz_sub_ui(*rop->bits, *op1.bits, (uint64_t) op2);
+  normalize_lbits(rop);
+}
+
+void and_bits(lbits *rop, const lbits op1, const lbits op2)
+{
+  assert(op1.len == op2.len);
+  rop->len = op1.len;
+  mpz_and(*rop->bits, *op1.bits, *op2.bits);
+}
+
+void or_bits(lbits *rop, const lbits op1, const lbits op2)
+{
+  assert(op1.len == op2.len);
+  rop->len = op1.len;
+  mpz_ior(*rop->bits, *op1.bits, *op2.bits);
+}
+
+void xor_bits(lbits *rop, const lbits op1, const lbits op2)
+{
+  assert(op1.len == op2.len);
+  rop->len = op1.len;
+  mpz_xor(*rop->bits, *op1.bits, *op2.bits);
+}
+
+void not_bits(lbits *rop, const lbits op)
+{
+  rop->len = op.len;
+  mpz_set(*rop->bits, *op.bits);
+  for (mp_bitcnt_t i = 0; i < op.len; i++) {
+    mpz_combit(*rop->bits, i);
+  }
+}
+
+void mults_vec(lbits *rop, const lbits op1, const lbits op2)
+{
+  return;
+}
+
+void mult_vec(lbits *rop, const lbits op1, const lbits op2)
+{
+  rop->len = op1.len * 2;
+  mpz_mul(*rop->bits, *op1.bits, *op2.bits);
+  normalize_lbits(rop); /* necessary? */
+}
+
+
+void zeros(lbits *rop, const sail_int op)
+{
+  rop->len = (mp_bitcnt_t) op;
+  mpz_set_ui(*rop->bits, 0);
+}
+
+void zero_extend(lbits *rop, const lbits op, const sail_int len)
+{
+  assert(op.len <= (uint64_t) len);
+  rop->len = (uint64_t) len;
+  mpz_set(*rop->bits, *op.bits);
+}
+
+fbits fast_zero_extend(const sbits op, const uint64_t n)
+{
+  return op.bits;
+}
+
+void sign_extend(lbits *rop, const lbits op, const sail_int len)
+{
+  assert(op.len <= (uint64_t) len);
+  rop->len = (uint64_t) len;
+  if(mpz_tstbit(*op.bits, op.len - 1)) {
+    mpz_set(*rop->bits, *op.bits);
+    for(mp_bitcnt_t i = rop->len - 1; i >= op.len; i--) {
+      mpz_setbit(*rop->bits, i);
+    }
+  } else {
+    mpz_set(*rop->bits, *op.bits);
+  }
+}
+
+fbits fast_sign_extend(const fbits op, const uint64_t n, const uint64_t m)
+{
+  uint64_t rop = op;
+  if (op & (UINT64_C(1) << (n - 1))) {
+    for (uint64_t i = m - 1; i >= n; i--) {
+      rop = rop | (UINT64_C(1) << i);
+    }
+    return rop;
+  } else {
+    return rop;
+  }
+}
+
+fbits fast_sign_extend2(const sbits op, const uint64_t m)
+{
+  uint64_t rop = op.bits;
+  if (op.bits & (UINT64_C(1) << (op.len - 1))) {
+    for (uint64_t i = m - 1; i >= op.len; i--) {
+      rop = rop | (UINT64_C(1) << i);
+    }
+    return rop;
+  } else {
+    return rop;
+  }
+}
+
+sail_int length_lbits(const lbits op)
+{
+  return (sail_int) op.len;
+}
+
+bool eq_bits(const lbits op1, const lbits op2)
+{
+  assert(op1.len == op2.len);
+  for (mp_bitcnt_t i = 0; i < op1.len; i++) {
+    if (mpz_tstbit(*op1.bits, i) != mpz_tstbit(*op2.bits, i)) return false;
+  }
+  return true;
+}
+
+bool EQUAL(lbits)(const lbits op1, const lbits op2)
+{
+  return eq_bits(op1, op2);
+}
+
+bool neq_bits(const lbits op1, const lbits op2)
+{
+  assert(op1.len == op2.len);
+  for (mp_bitcnt_t i = 0; i < op1.len; i++) {
+    if (mpz_tstbit(*op1.bits, i) != mpz_tstbit(*op2.bits, i)) return true;
+  }
+  return false;
+}
+
+void vector_subrange_lbits(lbits *rop,
+			       const lbits op,
+			       const sail_int n_mpz,
+			       const sail_int m_mpz)
+{
+  uint64_t n = (uint64_t) n_mpz;
+  uint64_t m = (uint64_t) m_mpz;
+
+  rop->len = n - (m - 1ul);
+  mpz_fdiv_q_2exp(*rop->bits, *op.bits, m);
+  normalize_lbits(rop);
+}
+
+void sail_truncate(lbits *rop, const lbits op, const sail_int len)
+{
+  rop->len = (mp_bitcnt_t) len;
+  mpz_set(*rop->bits, *op.bits);
+  normalize_lbits(rop);
+}
+
+void sail_truncateLSB(lbits *rop, const lbits op, const sail_int len)
+{
+  uint64_t rlen = (uint64_t) len;
+  assert(op.len >= rlen);
+  rop->len = rlen;
+  // similar to vector_subrange_lbits above -- right shift LSBs away
+  mpz_fdiv_q_2exp(*rop->bits, *op.bits, op.len - rlen);
+  normalize_lbits(rop);
+}
+
+fbits bitvector_access(const lbits op, const sail_int n)
+{
+  return (fbits) mpz_tstbit(*op.bits, (uint64_t) n);
+}
+
+sail_int sail_unsigned(const lbits op)
+{
+  return (sail_int) mpz_get_ui(*op.bits);
+}
+
+sail_int sail_signed(const lbits op)
+{
+  if (op.len <= 64) {
+    uint64_t b = mpz_get_ui(*op.bits);
+    uint64_t sign_bit = UINT64_C(1) << (op.len - UINT64_C(1));
+    if ((b & sign_bit) > 0) {
+      return ((sail_int) (b & ~sign_bit)) - ((sail_int) sign_bit);
+    } else {
+      return (sail_int) b;
+    }
+  } else if (op.len <= 128) {
+    uint64_t b_lo = mpz_get_ui(*op.bits);
+    mpz_t tmp;
+    mpz_init(tmp);
+    mpz_tdiv_q_2exp(tmp, *op.bits, 64);
+    uint64_t b_hi = mpz_get_ui(tmp);
+    mpz_clear(tmp);
+    uint64_t sign_bit = UINT64_C(1) << (op.len - UINT64_C(65));
+    if (b_hi & sign_bit) {
+      unsigned __int128 b = b_hi & ~sign_bit;
+      b <<= 64;
+      b |= (unsigned __int128) b_lo;
+      unsigned __int128 sb = (unsigned __int128) sign_bit << (unsigned __int128) 64;
+      return (sail_int) b + (sail_int) (~sb + 1);
+    } else {
+      unsigned __int128 b = b_hi;
+      b <<= 64;
+      b |= (unsigned __int128) b_lo;
+      return (__int128) b;
+    }
+  } else {
+    printf("sail_signed >128\n");
+    exit(1);
+  }
+}
+
+mach_int fast_unsigned(const fbits op)
+{
+  return (mach_int) op;
+}
+
+mach_int fast_signed(const fbits op, const uint64_t n)
+{
+  if (op & (UINT64_C(1) << (n - 1))) {
+    uint64_t rop = op & ~(UINT64_C(1) << (n - 1));
+    return (mach_int) (rop - (UINT64_C(1) << (n - 1)));
+  } else {
+    return (mach_int) op;
+  }
+}
+
+void append(lbits *rop, const lbits op1, const lbits op2)
+{
+  rop->len = op1.len + op2.len;
+  mpz_mul_2exp(*rop->bits, *op1.bits, op2.len);
+  mpz_ior(*rop->bits, *rop->bits, *op2.bits);
+}
+
+sbits append_sf(const sbits op1, const fbits op2, const uint64_t len)
+{
+  sbits rop;
+  rop.bits = (op1.bits << len) | op2;
+  rop.len = op1.len + len;
+  return rop;
+}
+
+sbits append_fs(const fbits op1, const uint64_t len, const sbits op2)
+{
+  sbits rop;
+  rop.bits = (op1 << op2.len) | op2.bits;
+  rop.len = len + op2.len;
+  return rop;
+}
+
+sbits append_ss(const sbits op1, const sbits op2)
+{
+  sbits rop;
+  rop.bits = (op1.bits << op2.len) | op2.bits;
+  rop.len = op1.len + op2.len;
+  return rop;
+}
+
+void replicate_bits(lbits *rop, const lbits op1, const sail_int op2)
+{
+  uint64_t op2_ui = (uint64_t) op2;
+  rop->len = op1.len * op2_ui;
+  mpz_set_ui(*rop->bits, 0);
+  for (int i = 0; i < op2_ui; i++) {
+    mpz_mul_2exp(*rop->bits, *rop->bits, op1.len);
+    mpz_ior(*rop->bits, *rop->bits, *op1.bits);
+  }
+}
+
+uint64_t fast_replicate_bits(const uint64_t shift, const uint64_t v, const int64_t times)
+{
+  uint64_t r = v;
+  for (int i = 1; i < times; ++i) {
+    r |= (r << shift);
+  }
+  return r;
+}
+
+// Takes a slice of the (two's complement) binary representation of
+// integer n, starting at bit start, and of length len. With the
+// argument in the following order:
+//
+// get_slice_int(len, n, start)
+//
+// For example:
+//
+// get_slice_int(8, 1680, 4) =
+//
+//                    11           0
+//                    V            V
+// get_slice_int(8, 0b0110_1001_0000, 4) = 0b0110_1001
+//                    <-------^
+//                    (8 bit) 4
+//
+__attribute__((target ("bmi2")))
+void get_slice_int(lbits *rop, const sail_int len, const sail_int n, const sail_int start)
+{
+  assert(len <= 128);
+  
+  unsigned __int128 nbits = (unsigned __int128) (n >> start);
+
+  if (len <= 64) {
+    mpz_set_ui(*rop->bits, _bzhi_u64((uint64_t) nbits, (uint64_t) len));
+    rop->len = (uint64_t) len;
+  } else {
+    print("get_slice_int");
+    exit(1);
+  }
+}
+
+sail_int set_slice_int(const sail_int len, const sail_int n, const sail_int start, const lbits slice)
+{
+  printf("set_slice_int");
+  exit(1);
+  return 0;
+}
+
+void update_lbits(lbits *rop, const lbits op, const sail_int n_mpz, const uint64_t bit)
+{
+  uint64_t n = (uint64_t) n_mpz;
+
+  mpz_set(*rop->bits, *op.bits);
+  rop->len = op.len;
+
+  if (bit == UINT64_C(0)) {
+    mpz_clrbit(*rop->bits, n);
+  } else {
+    mpz_setbit(*rop->bits, n);
+  }
+}
+
+void vector_update_subrange_lbits(lbits *rop,
+				 const lbits op,
+				 const sail_int n_mpz,
+				 const sail_int m_mpz,
+				 const lbits slice)
+{
+  uint64_t n = (uint64_t) n_mpz;
+  uint64_t m = (uint64_t) m_mpz;
+
+  mpz_set(*rop->bits, *op.bits);
+  rop->len = op.len;
+
+  for (uint64_t i = 0; i < n - (m - 1ul); i++) {
+    if (mpz_tstbit(*slice.bits, i)) {
+      mpz_setbit(*rop->bits, i + m);
+    } else {
+      mpz_clrbit(*rop->bits, i + m);
+    }
+  }
+}
+
+fbits fast_update_subrange(const fbits op,
+			   const mach_int n,
+			   const mach_int m,
+			   const fbits slice)
+{
+  fbits rop = op;
+  for (mach_int i = 0; i < n - (m - UINT64_C(1)); i++) {
+    uint64_t bit = UINT64_C(1) << ((uint64_t) i);
+    if (slice & bit) {
+      rop |= (bit << m);
+    } else {
+      rop &= ~(bit << m);
+    }
+  }
+  return rop;
+}
+
+__attribute__((target ("bmi2")))
+void slice(lbits *rop, const lbits op, const sail_int start_big, const sail_int len_big)
+{
+  uint64_t start = (uint64_t) start_big;
+  uint64_t len = (uint64_t) len_big;
+
+  if (len + start <= 64) {
+    mpz_set_ui(*rop->bits, _bzhi_u64(mpz_get_ui(*op.bits) >> start, len));
+    rop->len = len;
+  } else {
+    mpz_set_ui(*rop->bits, 0);
+    rop->len = len;
+
+    for (uint64_t i = 0; i < len; i++) {
+      if (mpz_tstbit(*op.bits, i + start)) mpz_setbit(*rop->bits, i);
+    }
+  }
+}
+
+__attribute__((target ("bmi2")))
+sbits sslice(const fbits op, const mach_int start, const mach_int len)
+{
+  sbits rop;
+  rop.bits = _bzhi_u64(op >> start, (uint64_t) len);
+  rop.len = (uint64_t) len;
+  return rop;
+}
+
+void set_slice(lbits *rop,
+	       const sail_int len_mpz,
+	       const sail_int slen_mpz,
+	       const lbits op,
+	       const sail_int start_mpz,
+	       const lbits slice)
+{
+  uint64_t start = (uint64_t) start_mpz;
+
+  mpz_set(*rop->bits, *op.bits);
+  rop->len = op.len;
+
+  for (uint64_t i = 0; i < slice.len; i++) {
+    if (mpz_tstbit(*slice.bits, i)) {
+      mpz_setbit(*rop->bits, i + start);
+    } else {
+      mpz_clrbit(*rop->bits, i + start);
+    }
+  }
+}
+
+void shift_bits_left(lbits *rop, const lbits op1, const lbits op2)
+{
+  rop->len = op1.len;
+  mpz_mul_2exp(*rop->bits, *op1.bits, mpz_get_ui(*op2.bits));
+  normalize_lbits(rop);
+}
+
+void shift_bits_right(lbits *rop, const lbits op1, const lbits op2)
+{
+  rop->len = op1.len;
+  mpz_tdiv_q_2exp(*rop->bits, *op1.bits, mpz_get_ui(*op2.bits));
+}
+
+/* FIXME */
+void shift_bits_right_arith(lbits *rop, const lbits op1, const lbits op2)
+{
+  rop->len = op1.len;
+  mp_bitcnt_t shift_amt = mpz_get_ui(*op2.bits);
+  mp_bitcnt_t sign_bit = op1.len - 1;
+  mpz_fdiv_q_2exp(*rop->bits, *op1.bits, shift_amt);
+  if(mpz_tstbit(*op1.bits, sign_bit) != 0) {
+    /* */
+    for(; shift_amt > 0; shift_amt--) {
+      mpz_setbit(*rop->bits, sign_bit - shift_amt + 1);
+    }
+  }
+}
+
+void shiftl(lbits *rop, const lbits op1, const sail_int op2)
+{
+  rop->len = op1.len;
+  mpz_mul_2exp(*rop->bits, *op1.bits, (uint64_t) op2);
+  normalize_lbits(rop);
+}
+
+void shiftr(lbits *rop, const lbits op1, const sail_int op2)
+{
+  rop->len = op1.len;
+  mpz_tdiv_q_2exp(*rop->bits, *op1.bits, (uint64_t) op2);
+}
+
+void reverse_endianness(lbits *rop, const lbits op)
+{
+  rop->len = op.len;
+  if (rop->len == 64ul) {
+    uint64_t x = mpz_get_ui(*op.bits);
+    x = (x & 0xFFFFFFFF00000000) >> 32 | (x & 0x00000000FFFFFFFF) << 32;
+    x = (x & 0xFFFF0000FFFF0000) >> 16 | (x & 0x0000FFFF0000FFFF) << 16;
+    x = (x & 0xFF00FF00FF00FF00) >> 8  | (x & 0x00FF00FF00FF00FF) << 8;
+    mpz_set_ui(*rop->bits, x);
+  } else if (rop->len == 32ul) {
+    uint64_t x = mpz_get_ui(*op.bits);
+    x = (x & 0xFFFF0000FFFF0000) >> 16 | (x & 0x0000FFFF0000FFFF) << 16;
+    x = (x & 0xFF00FF00FF00FF00) >> 8  | (x & 0x00FF00FF00FF00FF) << 8;
+    mpz_set_ui(*rop->bits, x);
+  } else if (rop->len == 16ul) {
+    uint64_t x = mpz_get_ui(*op.bits);
+    x = (x & 0xFF00FF00FF00FF00) >> 8  | (x & 0x00FF00FF00FF00FF) << 8;
+    mpz_set_ui(*rop->bits, x);
+  } else if (rop->len == 8ul) {
+    mpz_set(*rop->bits, *op.bits);
+  } else {
+    mpz_t tmp1;
+    mpz_t tmp2;
+    mpz_init(tmp1);
+    mpz_init(tmp2);
+    
+    /* For other numbers of bytes we reverse the bytes.
+     * XXX could use mpz_import/export for this. */
+    mpz_set_ui(tmp1, 0xff); // byte mask
+    mpz_set_ui(*rop->bits, 0); // reset accumulator for result
+    for(mp_bitcnt_t byte = 0; byte < op.len; byte+=8) {
+      mpz_tdiv_q_2exp(tmp2, *op.bits, byte); // shift byte to bottom
+      mpz_and(tmp2, tmp2, tmp1); // and with mask
+      mpz_mul_2exp(*rop->bits, *rop->bits, 8); // shift result left 8
+      mpz_ior(*rop->bits, *rop->bits, tmp2); // or byte into result
+    }
+  }
+}
+
+bool eq_sbits(const sbits op1, const sbits op2)
+{
+  return op1.bits == op2.bits;
+}
+
+bool neq_sbits(const sbits op1, const sbits op2)
+{
+  return op1.bits != op2.bits;
+}
+
+__attribute__((target ("bmi2")))
+sbits not_sbits(const sbits op)
+{
+  sbits rop;
+  rop.bits = (~op.bits) & _bzhi_u64(UINT64_MAX, op.len);
+  rop.len = op.len;
+  return rop;
+}
+
+sbits xor_sbits(const sbits op1, const sbits op2)
+{
+  sbits rop;
+  rop.bits = op1.bits ^ op2.bits;
+  rop.len = op1.len;
+  return rop;
+}
+
+sbits or_sbits(const sbits op1, const sbits op2)
+{
+  sbits rop;
+  rop.bits = op1.bits | op2.bits;
+  rop.len = op1.len;
+  return rop;
+}
+
+sbits and_sbits(const sbits op1, const sbits op2)
+{
+  sbits rop;
+  rop.bits = op1.bits & op2.bits;
+  rop.len = op1.len;
+  return rop;
+}
+
+__attribute__((target ("bmi2")))
+sbits add_sbits(const sbits op1, const sbits op2)
+{
+  sbits rop;
+  rop.bits = (op1.bits + op2.bits) & _bzhi_u64(UINT64_MAX, op1.len);
+  rop.len = op1.len;
+  return rop;
+}
+
+__attribute__((target ("bmi2")))
+sbits sub_sbits(const sbits op1, const sbits op2)
+{
+  sbits rop;
+  rop.bits = (op1.bits - op2.bits) & _bzhi_u64(UINT64_MAX, op1.len);
+  rop.len = op1.len;
+  return rop;
+}
+
+/* ***** Sail Reals ***** */
+
+void CREATE(real)(real *rop)
+{
+  mpq_init(*rop);
+}
+
+void RECREATE(real)(real *rop)
+{
+  mpq_set_ui(*rop, 0, 1);
+}
+
+void KILL(real)(real *rop)
+{
+  mpq_clear(*rop);
+}
+
+void COPY(real)(real *rop, const real op)
+{
+  mpq_set(*rop, op);
+}
+
+void UNDEFINED(real)(real *rop, unit u)
+{
+  mpq_set_ui(*rop, 0, 1);
+}
+
+void neg_real(real *rop, const real op)
+{
+  mpq_neg(*rop, op);
+}
+
+void mult_real(real *rop, const real op1, const real op2) {
+  mpq_mul(*rop, op1, op2);
+}
+
+void sub_real(real *rop, const real op1, const real op2)
+{
+  mpq_sub(*rop, op1, op2);
+}
+
+void add_real(real *rop, const real op1, const real op2)
+{
+  mpq_add(*rop, op1, op2);
+}
+
+void div_real(real *rop, const real op1, const real op2)
+{
+  mpq_div(*rop, op1, op2);
+}
+
+#define SQRT_PRECISION 30
+
+/*
+ * sqrt_real first checks if op has the form n/1 - in which case, if n
+ * is a perfect square (i.e. it's square root is an integer), then it
+ * will return the exact square root using mpz_sqrt. If that's not the
+ * case we use the Babylonian method to calculate the square root to
+ * SQRT_PRECISION decimal places.
+ */
+void sqrt_real(real *rop, const real op)
+{
+  /* First check if op is a perfect square and use mpz_sqrt if so */
+  if (mpz_cmp_ui(mpq_denref(op), 1) == 0 && mpz_perfect_square_p(mpq_numref(op))) {
+    mpz_sqrt(mpq_numref(*rop), mpq_numref(op));
+    mpz_set_ui(mpq_denref(*rop), 1);
+    return;
+  }
+
+  mpq_t tmp;
+  mpz_t tmp_z;
+  mpq_t p; /* previous estimate, p */
+  mpq_t n; /* next estimate, n */
+  /* convergence is the precision (in decimal places) we want to reach as a fraction 1/(10^precision) */
+  mpq_t convergence;
+
+  mpq_init(tmp);
+  mpz_init(tmp_z);
+  mpq_init(p);
+  mpq_init(n);
+  mpq_init(convergence);
+
+  /* calculate an initial guess using mpz_sqrt */
+  mpz_cdiv_q(tmp_z, mpq_numref(op), mpq_denref(op));
+  mpz_sqrt(tmp_z, tmp_z);
+  mpq_set_z(p, tmp_z);
+
+  /* initialise convergence based on SQRT_PRECISION */
+  mpz_set_ui(tmp_z, 10);
+  mpz_pow_ui(tmp_z, tmp_z, SQRT_PRECISION);
+  mpz_set_ui(mpq_numref(convergence), 1);
+  mpq_set_den(convergence, tmp_z);
+
+  while (true) {
+    // n = (p + op / p) / 2
+    mpq_div(tmp, op, p);
+    mpq_add(tmp, tmp, p);
+    mpq_div_2exp(n, tmp, 1);
+
+    /* calculate the difference between n and p */
+    mpq_sub(tmp, p, n);
+    mpq_abs(tmp, tmp);
+
+    /* if the difference is small enough, return */
+    if (mpq_cmp(tmp, convergence) < 0) {
+      mpq_set(*rop, n);
+      break;
+    }
+
+    mpq_swap(n, p);
+  }
+
+  mpq_clear(tmp);
+  mpz_clear(tmp_z);
+  mpq_clear(p);
+  mpq_clear(n);
+  mpq_clear(convergence);
+}
+
+void abs_real(real *rop, const real op)
+{
+  mpq_abs(*rop, op);
+}
+
+sail_int round_up(const real op)
+{
+  mpz_t rop;
+  mpz_init(rop);
+  mpz_cdiv_q(rop, mpq_numref(op), mpq_denref(op));
+  sail_int r = mpz_get_si(rop);
+  mpz_clear(rop);
+  return r;
+}
+
+sail_int round_down(const real op)
+{
+  mpz_t rop;
+  mpz_init(rop);
+  mpz_fdiv_q(rop, mpq_numref(op), mpq_denref(op));
+  sail_int r = mpz_get_si(rop);
+  mpz_clear(rop);
+  return r;
+}
+
+void to_real(real *rop, const sail_int op)
+{
+  mpz_t op_mpz;
+  mpz_init_set_si128(op_mpz, op);
+  
+  mpq_set_z(*rop, op_mpz);
+  mpq_canonicalize(*rop);
+
+  mpz_clear(op_mpz);
+}
+
+bool EQUAL(real)(const real op1, const real op2)
+{
+  return mpq_cmp(op1, op2) == 0;
+}
+
+bool lt_real(const real op1, const real op2)
+{
+  return mpq_cmp(op1, op2) < 0;
+}
+
+bool gt_real(const real op1, const real op2)
+{
+  return mpq_cmp(op1, op2) > 0;
+}
+
+bool lteq_real(const real op1, const real op2)
+{
+  return mpq_cmp(op1, op2) <= 0;
+}
+
+bool gteq_real(const real op1, const real op2)
+{
+  return mpq_cmp(op1, op2) >= 0;
+}
+
+void real_power(real *rop, const real base, const sail_int exp)
+{
+  int64_t exp_si = (int64_t) exp;
+
+  mpz_set_ui(mpq_numref(*rop), 1);
+  mpz_set_ui(mpq_denref(*rop), 1);
+
+  real b;
+  mpq_init(b);
+  mpq_set(b, base);
+  int64_t pexp = llabs(exp_si);
+  while (pexp != 0) {
+    // invariant: rop * b^pexp == base^abs(exp)
+    if (pexp & 1) { // b^(e+1) = b * b^e
+      mpq_mul(*rop, *rop, b);
+      pexp -= 1;
+    } else { // b^(2e) = (b*b)^e
+      mpq_mul(b, b, b);
+      pexp >>= 1;
+    }
+  }
+  if (exp_si < 0) {
+    mpq_inv(*rop, *rop);
+  }
+  mpq_clear(b);
+}
+
+void CREATE_OF(real, sail_string)(real *rop, const sail_string op)
+{
+  mpq_init(*rop);
+  CONVERT_OF(real, sail_string)(rop, op);
+}
+
+void CONVERT_OF(real, sail_string)(real *rop, const sail_string op)
+{
+  int decimal;
+  int total;
+  mpz_t tmp1;
+  mpz_t tmp2;
+  mpz_t tmp3;
+  mpq_t tmp_real;
+  mpz_init(tmp1);
+  mpz_init(tmp2);
+  mpz_init(tmp3);
+  mpq_init(tmp_real);
+  
+  gmp_sscanf(op, "%Zd.%n%Zd%n", tmp1, &decimal, tmp2, &total);
+
+  int len = total - decimal;
+  mpz_ui_pow_ui(tmp3, 10, len);
+  mpz_set(mpq_numref(*rop), tmp2);
+  mpz_set(mpq_denref(*rop), tmp3);
+  mpq_canonicalize(*rop);
+  mpz_set(mpq_numref(tmp_real), tmp1);
+  mpz_set_ui(mpq_denref(tmp_real), 1);
+  mpq_add(*rop, *rop, tmp_real);
+
+  mpz_clear(tmp1);
+  mpz_clear(tmp2);
+  mpz_clear(tmp3);
+  mpq_clear(tmp_real);
+}
+
+unit print_real(const sail_string str, const real op)
+{
+  gmp_printf("%s%Qd\n", str, op);
+  return UNIT;
+}
+
+unit prerr_real(const sail_string str, const real op)
+{
+  gmp_fprintf(stderr, "%s%Qd\n", str, op);
+  return UNIT;
+}
+
+void random_real(real *rop, const unit u)
+{
+  if (rand() & 1) {
+    mpz_set_si(mpq_numref(*rop), rand());
+  } else {
+    mpz_set_si(mpq_numref(*rop), -rand());
+  }
+  mpz_set_si(mpq_denref(*rop), rand());
+  mpq_canonicalize(*rop);
+}
+
+/* ***** Printing functions ***** */
+
+void string_of_int(sail_string *str, const sail_int i)
+{
+  free(*str);
+  //gmp_asprintf(str, "%Zd", i);
+}
+
+/* asprintf is a GNU extension, but it should exist on BSD */
+void string_of_fbits(sail_string *str, const fbits op)
+{
+  free(*str);
+  int bytes = asprintf(str, "0x%" PRIx64, op);
+  if (bytes == -1) {
+    fprintf(stderr, "Could not print bits 0x%" PRIx64 "\n", op);
+  }
+}
+
+void string_of_lbits(sail_string *str, const lbits op)
+{
+  free(*str);
+  if ((op.len % 4) == 0) {
+    gmp_asprintf(str, "0x%*0ZX", op.len / 4, *op.bits);
+  } else {
+    *str = (char *) malloc((op.len + 3) * sizeof(char));
+    (*str)[0] = '0';
+    (*str)[1] = 'b';
+    for (int i = 1; i <= op.len; ++i) {
+      (*str)[i + 1] = mpz_tstbit(*op.bits, op.len - i) + 0x30;
+    }
+    (*str)[op.len + 2] = '\0';
+  }
+}
+
+void decimal_string_of_fbits(sail_string *str, const fbits op)
+{
+  free(*str);
+  int bytes = asprintf(str, "%" PRId64, op);
+  if (bytes == -1) {
+    fprintf(stderr, "Could not print bits %" PRId64 "\n", op);
+  }
+}
+
+void decimal_string_of_lbits(sail_string *str, const lbits op)
+{
+  free(*str);
+  gmp_asprintf(str, "%Z", *op.bits);
+}
+
+void fprint_bits(const sail_string pre,
+		 const lbits op,
+		 const sail_string post,
+		 FILE *stream)
+{
+  fputs(pre, stream);
+
+  if (op.len % 4 == 0) {
+    fputs("0x", stream);
+    mpz_t buf;
+    mpz_init_set(buf, *op.bits);
+
+    char *hex = malloc((op.len / 4) * sizeof(char));
+
+    for (int i = 0; i < op.len / 4; ++i) {
+      char c = (char) ((0xF & mpz_get_ui(buf)) + 0x30);
+      hex[i] = (c < 0x3A) ? c : c + 0x7;
+      mpz_fdiv_q_2exp(buf, buf, 4);
+    }
+
+    for (int i = op.len / 4; i > 0; --i) {
+      fputc(hex[i - 1], stream);
+    }
+
+    free(hex);
+    mpz_clear(buf);
+  } else {
+    fputs("0b", stream);
+    for (int i = op.len; i > 0; --i) {
+      fputc(mpz_tstbit(*op.bits, i - 1) + 0x30, stream);
+    }
+  }
+
+  fputs(post, stream);
+}
+
+unit print_bits(const sail_string str, const lbits op)
+{
+  fprint_bits(str, op, "\n", stdout);
+  return UNIT;
+}
+
+unit prerr_bits(const sail_string str, const lbits op)
+{
+  fprint_bits(str, op, "\n", stderr);
+  return UNIT;
+}
+
+unit print(const sail_string str)
+{
+  printf("%s", str);
+  return UNIT;
+}
+
+unit print_endline(const sail_string str)
+{
+  printf("%s\n", str);
+  return UNIT;
+}
+
+unit prerr(const sail_string str)
+{
+  fprintf(stderr, "%s", str);
+  return UNIT;
+}
+
+unit prerr_endline(const sail_string str)
+{
+  fprintf(stderr, "%s\n", str);
+  return UNIT;
+}
+
+unit print_int(const sail_string str, const sail_int op)
+{
+  mpz_t op_mpz;
+  mpz_init_set_si128(op_mpz, op);
+  
+  fputs(str, stdout);
+  mpz_out_str(stdout, 10, op_mpz);
+  putchar('\n');
+
+  mpz_clear(op_mpz);
+  return UNIT;
+}
+
+unit prerr_int(const sail_string str, const sail_int op)
+{
+  fputs(str, stderr);
+  //mpz_out_str(stderr, 10, op);
+  fputs("\n", stderr);
+  return UNIT;
+}
+
+unit sail_putchar(const sail_int op)
+{
+  char c = (char) op;
+  putchar(c);
+  fflush(stdout);
+  return UNIT;
+}
+
+sail_int get_time_ns(const unit u)
+{
+  struct timespec t;
+  clock_gettime(CLOCK_REALTIME, &t);
+  __int128 rop = (__int128) t.tv_sec;
+  rop *= 1000000000;
+  rop += (__int128) t.tv_nsec;
+  return rop;
+}
+
+// Monomorphisation
+sail_int make_the_value(const sail_int op)
+{
+  return op;
+}
+
+sail_int size_itself_int(const sail_int op)
+{
+  return op;
+}