WIP

src/libvfcinstrumentonline/Makefile.am

@@ -10,6 +10,6 @@ libvfcinstrumentonline_la_CXXFLAGS = @LLVM_CPPFLAGS@
 libvfcinstrumentonline_la_LDFLAGS = @LLVM_LDFLAGS@
 libvfcinstrumentonline_la_SOURCES = libVFCInstrumentOnline.cpp
 
-librand_la_CXXFLAGS = @LLVM_CPPFLAGS@
-librand_la_LDFLAGS = @LLVM_LDFLAGS@
+librand_la_CFLAGS = -I@INTERFLOP_INCLUDEDIR@
+librand_la_LDFLAGS = -lm
 librand_la_SOURCES = rand.c

src/libvfcinstrumentonline/libVFCInstrumentOnline.cpp

@@ -76,6 +76,11 @@ static cl::opt<bool> VfclibInstVerbose("vfclibinst-verbose",
                                        cl::value_desc("Verbose"),
                                        cl::init(false));
 
+static cl::opt<std::string>
+    VfclibInstMode("vfclibinst-mode",
+                   cl::desc("Instrumentation mode: up-down or sr"),
+                   cl::value_desc("Mode"), cl::init("up-down"));
+
 /* pointer that hold the vfcwrapper Module */
 // static Module *vfcwrapperM = nullptr;
 
@@ -349,9 +354,40 @@ struct VfclibInst : public ModulePass {
     }
   }
 
-  /* Replace arithmetic instructions with MCA */
-  Value *replaceArithmeticWithMCACall(IRBuilder<> &Builder, Instruction *I,
-                                      std::set<User *> &users) {
+  std::string sr_hook_name(Instruction &I) {
+    switch (I.getOpcode()) {
+    case Instruction::FAdd:
+      return "add2";
+    case Instruction::FSub:
+      // In LLVM IR the FSub instruction is used to represent FNeg
+      return "sub2";
+    case Instruction::FMul:
+      return "mul2";
+    case Instruction::FDiv:
+      return "div2";
+    }
+  }
+
+  Value *insertSROpCall(IRBuilder<> &Builder, Instruction *I) {
+    Module *M = I->getModule();
+    std::string typestring = validTypesMap[I->getType()->getTypeID()];
+    std::string srName = sr_hook_name(*I) + "_" + typestring;
+    FunctionType *funcType =
+        FunctionType::get(I->getType(), {I->getType(), I->getType()}, false);
+    FunctionCallee SR = M->getOrInsertFunction(srName, funcType);
+    return Builder.CreateCall(SR, {I->getOperand(0), I->getOperand(1)});
+  }
+
+  Value *replaceArithmeticWithMCACall_SR(IRBuilder<> &Builder, Instruction *I,
+                                         std::set<User *> &users) {
+    Type *fpAsIntTy = getFPAsIntType(I->getType());
+    Value *value = insertSROpCall(Builder, I);
+    return value;
+  }
+
+  Value *replaceArithmeticWithMCACall_UpOrDown(IRBuilder<> &Builder,
+                                               Instruction *I,
+                                               std::set<User *> &users) {
     Type *fpAsIntTy = getFPAsIntType(I->getType());
     Value *randomBit = insertRNGCall(Builder, I);
     users.insert(static_cast<User *>(randomBit));
@@ -363,6 +399,19 @@ struct VfclibInst : public ModulePass {
     return newFP;
   }
 
+  /* Replace arithmetic instructions with MCA */
+  Value *replaceArithmeticWithMCACall(IRBuilder<> &Builder, Instruction *I,
+                                      std::set<User *> &users) {
+    if (VfclibInstMode == "up-down") {
+      return replaceArithmeticWithMCACall_UpOrDown(Builder, I, users);
+    } else if (VfclibInstMode == "sr") {
+      return replaceArithmeticWithMCACall_SR(Builder, I, users);
+    } else {
+      errs() << "Unsupported mode: " << VfclibInstMode << "\n";
+      report_fatal_error("libVFCInstrument fatal error");
+    }
+  }
+
   Value *replaceWithMCACall(Module &M, Instruction *I, Fops opCode,
                             std::set<User *> &users) {
     if (not isValidInstruction(I)) {

src/libvfcinstrumentonline/rand.c

@@ -1,10 +1,13 @@
+#include <math.h>
 #include <pthread.h>
 #include <stdbool.h>
 #include <stdint.h>
 #include <sys/syscall.h>
 #include <sys/time.h>
 #include <unistd.h>
 
+#include "interflop/common/float_utils.h"
+
 #define __INTERNAL_RNG_STATE xoroshiro_state
 
 typedef uint64_t xoroshiro_state[2];
@@ -48,6 +51,15 @@ int64_t get_rand_double(double a) {
   return (get_rand_uint64() < half_max_int) ? 1 : -1;
 }
 
+static inline double get_rand_double01() {
+  uint64_t x = get_rand_uint64();
+  const union {
+    uint64_t i;
+    double d;
+  } u = {.i = UINT64_C(0x3FF) << 52 | x >> 12};
+  return u.d - 1.0;
+}
+
 __attribute__((constructor)) static void init(void) {
   if (already_initialized == 0) {
     already_initialized = 1;
@@ -60,4 +72,197 @@ __attribute__((constructor)) static void init(void) {
   seed = t1.tv_sec ^ t1.tv_usec ^ syscall(__NR_gettid);
   rng_state[0] = next_seed(seed);
   rng_state[1] = next_seed(seed);
+}
+
+uint32_t get_exponent_b64(double x) { return _get_exponent_binary64(x); }
+uint32_t get_exponent_b32(float x) { return _get_exponent_binary32(x); }
+
+double predecessor_b64(double x) {
+  uint64_t x_bits = *(uint64_t *)&x;
+  x_bits--;
+  return *(double *)&x_bits;
+}
+
+float predecessor_b32(float x) {
+  uint32_t x_bits = *(uint32_t *)&x;
+  x_bits--;
+  return *(float *)&x_bits;
+}
+
+double abs_b64(double x) {
+  uint64_t x_bits = *(uint64_t *)&x;
+  x_bits &= 0x7FFFFFFFFFFFFFFF;
+  return *(double *)&x_bits;
+}
+
+float abs_b32(float x) {
+  uint32_t x_bits = *(uint32_t *)&x;
+  x_bits &= 0x7FFFFFFF;
+  return *(float *)&x_bits;
+}
+
+double pow2_b64(uint32_t n) {
+  uint64_t n_bits = (uint64_t)(n + 1023) << 52;
+  return *(double *)&n_bits;
+}
+
+float pow2_b32(uint32_t n) {
+  uint32_t n_bits = (n + 127) << 23;
+  return *(float *)&n_bits;
+}
+
+double sr_round_b64(double sigma, double tau, double z) {
+  // Compute round
+  double eps = 0x1.0p-52;
+  bool sign_tau, sign_sigma;
+  sign_tau = tau < 0;
+  sign_sigma = sigma < 0;
+  uint32_t eta;
+  if (sign_tau != sign_sigma) {
+    eta = get_exponent_b64(predecessor_b64(abs_b64(sigma)));
+  } else {
+    eta = get_exponent_b64(sigma);
+  }
+  double ulp = sign_tau * pow2_b64(eta) * eps;
+  double pi = ulp * z;
+  double round;
+  if (abs_b64(tau + pi) >= abs_b64(ulp)) {
+    round = ulp;
+  } else {
+    round = 0;
+  }
+  return round;
+}
+
+float sr_round_b32(float sigma, float tau, float z) {
+  // Compute round
+  float eps = 0x1.0p-23;
+  bool sign_tau, sign_sigma;
+  sign_tau = tau < 0;
+  sign_sigma = sigma < 0;
+  uint32_t eta;
+  if (sign_tau != sign_sigma) {
+    eta = get_exponent_b32(predecessor_b32(abs_b32(sigma)));
+  } else {
+    eta = get_exponent_b32(sigma);
+  }
+  float ulp = sign_tau * pow2_b32(eta) * eps;
+  float pi = ulp * z;
+  float round;
+  if (abs_b32(tau + pi) >= abs_b32(ulp)) {
+    round = ulp;
+  } else {
+    round = 0;
+  }
+  return round;
+}
+
+double twosum_b64(double a, double b, double *tau, double *sigma) {
+  // Compute tau and sigma
+  *sigma = a + b;
+  double z = *sigma - a;
+  *tau = (a - (*sigma - z)) + (b - z);
+}
+
+float twosum_b32(float a, float b, float *tau, float *sigma) {
+  // Compute tau and sigma
+  *sigma = a + b;
+  float z = *sigma - a;
+  *tau = (a - (*sigma - z)) + (b - z);
+}
+
+double twoprodfma_b64(double a, double b, double *tau, double *sigma) {
+  // Compute tau and sigma
+  *sigma = a * b;
+  *tau = __builtin_fma(a, b, -(*sigma));
+}
+
+float twoprodfma_b32(float a, float b, float *tau, float *sigma) {
+  // Compute tau and sigma
+  *sigma = a * b;
+  *tau = __builtin_fmaf(a, b, -(*sigma));
+}
+
+double add2_double(double a, double b) {
+  // compute SR(a+b)
+  double z, tau, sigma, round;
+  z = get_rand_double01(); // return float in [0,1)
+  twosum_b64(a, b, &tau, &sigma);
+  round = sr_round_b64(sigma, tau, z);
+  return sigma + round;
+}
+
+double sub2_double(float a, float b) { return add2_double(a, -b); }
+
+double mul2_double(double a, double b) {
+  // if a and b satisfy the condition (5.1), compute SR(a*b)
+  double z, tau, sigma, round;
+  z = get_rand_double01(); // return float in [0,1)
+  twoprodfma_b64(a, b, &tau, &sigma);
+  round = sr_round_b64(sigma, tau, z);
+  return sigma + round;
+}
+
+double div2_double(double a, double b) {
+  // compute SR(a/b)
+  double z, sigma, tau, round;
+  z = get_rand_double01(); // return float in [0,1)
+  sigma = a / b;
+  tau = -sigma * b + a;
+  tau = tau / b;
+  round = sr_round_b64(sigma, tau, z);
+  return sigma + round;
+}
+
+double sqrt2_double(double a) {
+  // compute SR(sqrt(a))
+  double z, sigma, tau, round;
+  z = get_rand_double01(); // return float in [0,1)
+  sigma = sqrt(a);
+  tau = -sigma * sigma + a;
+  tau = tau / (2 * sigma);
+  round = sr_round_b64(sigma, tau, z);
+  return sigma + round;
+}
+
+float add2_float(float a, float b) {
+  // compute SR(a+b)
+  float z, tau, sigma, round;
+  z = get_rand_double01(); // return float in [0,1)
+  twosum_b32(a, b, &tau, &sigma);
+  round = sr_round_b32(sigma, tau, z);
+  return sigma + round;
+}
+
+float sub2_float(float a, float b) { return add2_float(a, -b); }
+
+float mul2_float(float a, float b) {
+  // if a and b satisfy the condition (5.1), compute SR(a*b)
+  float z, tau, sigma, round;
+  z = get_rand_double01(); // return float in [0,1)
+  twoprodfma_b32(a, b, &tau, &sigma);
+  round = sr_round_b32(sigma, tau, z);
+  return sigma + round;
+}
+
+float div2_float(float a, float b) {
+  // compute SR(a/b)
+  float z, sigma, tau, round;
+  z = get_rand_double01(); // return float in [0,1)
+  sigma = a / b;
+  tau = -sigma * b + a;
+  tau = tau / b;
+  round = sr_round_b32(sigma, tau, z);
+  return sigma + round;
+}
+
+float sqrt2_float(float a) {
+  // compute SR(sqrt(a))
+  float z, sigma, tau, round;
+  z = get_rand_double01(); // return float in [0,1)
+  sigma = sqrtf(a);
+  tau = -sigma * sigma + a;
+  tau = tau / (2 * sigma);
+  round = sr_round_b32(sigma, tau, z);
+  return sigma + round;
 }

tests/test_online_instrumentation/Makefile

@@ -8,8 +8,8 @@ endif
 
 
 CC=verificarlo-c
-CFLAGS=-g -Wall -DREAL=$(type) $(optimization) --function=operator --online-instrumentation # --show-cmd --verbose
-LDFLAGS=--online-instrumentation
+CFLAGS=-g -Wall -DREAL=$(type) $(optimization) --function=operator --online-instrumentation=up-down # --show-cmd --verbose
+LDFLAGS=--online-instrumentation=up-down
 SOURCE=test.c
 BINARY=test_$(type)_$(optimization)
 

tests/test_online_instrumentation_sr/check_status.py

@@ -0,0 +1,5 @@
+import numpy as np
+import sys
+with open(sys.argv[1]) as fi:
+    x = [float.fromhex(l.strip()) for l in fi]
+    print(int(len(set(x)) == 1))

tests/test_online_instrumentation_sr/clean.sh

@@ -0,0 +1,3 @@
+#!/bin/bash
+
+rm -Rf *~ test_* *.ll *.o .*.o *.txt

tests/test_online_instrumentation_sr/test.c

@@ -0,0 +1,55 @@
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "interflop/common/float_const.h"
+
+#ifndef REAL
+#warning "REAL is not defined, defaulting to double"
+#define REAL double
+#endif
+
+typedef union {
+  double f64;
+  int64_t s64;
+} binary64;
+
+typedef union {
+  float f32;
+  int32_t s32;
+} binary32;
+
+#define OPERATOR(a, b)                                                         \
+  switch (op) {                                                                \
+  case '+':                                                                    \
+    return a + b;                                                              \
+  case '-':                                                                    \
+    return a - b;                                                              \
+  case 'x':                                                                    \
+    return a * b;                                                              \
+  case '/':                                                                    \
+    return a / b;                                                              \
+  default:                                                                     \
+    fprintf(stderr, "Error: unknown operation\n");                             \
+    abort();                                                                   \
+  }
+
+__attribute__((noinline)) REAL operator(char op, REAL a, REAL b) {
+  OPERATOR(a, b)
+}
+
+int main(int argc, const char *argv[]) {
+
+  if (argc != 4) {
+    fprintf(stderr, "usage: ./test <op> a b\n");
+    exit(1);
+  }
+
+  const char op = argv[1][0];
+  REAL a = atof(argv[2]);
+  REAL b = atof(argv[3]);
+
+  printf("%.13a\n", operator(op, a, b));
+
+  return EXIT_SUCCESS;
+}