diff --git a/src/r_interface.cpp b/src/r_interface.cpp
index 1cbbff4..363c127 100644
--- a/src/r_interface.cpp
+++ b/src/r_interface.cpp
@@ -1,8 +1,7 @@
-// needs to go first to avoid omp.h conflicts with Rinternals.h macros
-#include "cpu/mdfs.h"
-
 #include "r_interface.h"
 
+#include "cpu/mdfs.h"
+
 #ifdef WITH_CUDA
 #include "gpu/cucubes.h"
 #endif
@@ -24,52 +23,52 @@ SEXP r_compute_max_ig(
         SEXP Rin_use_cuda)
 {
     #ifndef WITH_CUDA
-    if (asLogical(Rin_use_cuda)) {
-        error("CUDA acceleration not compiled");
+    if (Rf_asLogical(Rin_use_cuda)) {
+        Rf_error("CUDA acceleration not compiled");
     }
     #endif
 
-    const int* dataDims = INTEGER(getAttrib(Rin_data, R_DimSymbol));
+    const int* dataDims = INTEGER(Rf_getAttrib(Rin_data, R_DimSymbol));
     int* contrastDataDims = nullptr;
-    if (!isNull(Rin_contrast_data)) {
-        contrastDataDims = INTEGER(getAttrib(Rin_contrast_data, R_DimSymbol));
+    if (!Rf_isNull(Rin_contrast_data)) {
+        contrastDataDims = INTEGER(Rf_getAttrib(Rin_contrast_data, R_DimSymbol));
     }
 
     const int obj_count = dataDims[0];
     const int variable_count = dataDims[1];
     int contrast_variable_count = 0;
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         contrast_variable_count = contrastDataDims[1];
     }
 
     #ifdef WITH_CUDA
-    if (asLogical(Rin_use_cuda)) {
-        SEXP Rout_max_igs = PROTECT(allocVector(REALSXP, variable_count));
+    if (Rf_asLogical(Rin_use_cuda)) {
+        SEXP Rout_max_igs = PROTECT(Rf_allocVector(REALSXP, variable_count));
 
         try {
             run_cucubes(
                     obj_count,
                     variable_count,
-                    asInteger(Rin_dimensions),
-                    asInteger(Rin_divisions),
-                    asInteger(Rin_discretizations),
-                    asInteger(Rin_seed),
-                    asReal(Rin_range),
-                    asReal(Rin_pseudocount),
+                    Rf_asInteger(Rin_dimensions),
+                    Rf_asInteger(Rin_divisions),
+                    Rf_asInteger(Rin_discretizations),
+                    Rf_asInteger(Rin_seed),
+                    Rf_asReal(Rin_range),
+                    Rf_asReal(Rin_pseudocount),
                     REAL(Rin_data),
                     INTEGER(Rin_decision),
                     REAL(Rout_max_igs));
         } catch (const cudaException& e) {
             // TODO: ensure cleanup inside library
-            error("CUDA exception: %s (in %s:%d)", cudaGetErrorString(e.code), e.file, e.line);
+            Rf_error("CUDA exception: %s (in %s:%d)", cudaGetErrorString(e.code), e.file, e.line);
         } catch (const NotImplementedException& e) {
             // TODO: is it possible to get this?
-            error("Not-implemented exception: %s", e.msg.c_str());
+            Rf_error("Not-implemented exception: %s", e.msg.c_str());
         }
 
         const int result_members_count = 1;
 
-        SEXP Rout_result = PROTECT(allocVector(VECSXP, result_members_count));
+        SEXP Rout_result = PROTECT(Rf_allocVector(VECSXP, result_members_count));
         SET_VECTOR_ELT(Rout_result, 0, Rout_max_igs);
 
         UNPROTECT(1 + result_members_count);
@@ -78,55 +77,55 @@ SEXP r_compute_max_ig(
     }
     #endif
 
-    const int discretizations = asInteger(Rin_discretizations);
-    const int divisions = asInteger(Rin_divisions);
+    const int discretizations = Rf_asInteger(Rin_discretizations);
+    const int divisions = Rf_asInteger(Rin_divisions);
 
     const int* decision = INTEGER(Rin_decision);
 
     RawData rawdata(RawDataInfo(obj_count, variable_count), REAL(Rin_data), decision);
     RawData* contrast_rawdata = nullptr;
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         contrast_rawdata = new RawData(RawDataInfo(obj_count, contrast_variable_count), REAL(Rin_contrast_data), nullptr);
     }
 
     std::unique_ptr<const DiscretizationInfo> dfi(new DiscretizationInfo(
-        asInteger(Rin_seed),
+        Rf_asInteger(Rin_seed),
         discretizations,
         divisions,
-        asReal(Rin_range)
+        Rf_asReal(Rin_range)
     ));
 
     MDFSInfo mdfs_info(
-        asInteger(Rin_dimensions),
+        Rf_asInteger(Rin_dimensions),
         divisions,
         discretizations,
-        asReal(Rin_pseudocount),
+        Rf_asReal(Rin_pseudocount),
         0.0f,
         INTEGER(Rin_interesting_vars),
-        length(Rin_interesting_vars),
-        asLogical(Rin_require_all_vars),
+        Rf_length(Rin_interesting_vars),
+        Rf_asLogical(Rin_require_all_vars),
         nullptr,
         false
     );
 
-    SEXP Rout_max_igs = PROTECT(allocVector(REALSXP, variable_count));
+    SEXP Rout_max_igs = PROTECT(Rf_allocVector(REALSXP, variable_count));
     SEXP Rout_contrast_max_igs = nullptr;
     SEXP Rout_tuples = nullptr;
     SEXP Rout_dids = nullptr;
 
-    if (!isNull(Rin_contrast_data)) {
-        Rout_contrast_max_igs = PROTECT(allocVector(REALSXP, contrast_variable_count));
+    if (!Rf_isNull(Rin_contrast_data)) {
+        Rout_contrast_max_igs = PROTECT(Rf_allocVector(REALSXP, contrast_variable_count));
     }
 
-    const bool return_tuples = asLogical(Rin_return_tuples);
+    const bool return_tuples = Rf_asLogical(Rin_return_tuples);
     MDFSOutput mdfs_output(MDFSOutputType::MaxIGs, mdfs_info.dimensions, variable_count, contrast_variable_count);
     if (return_tuples) {
-        Rout_tuples = PROTECT(allocMatrix(INTSXP, mdfs_info.dimensions, variable_count));
-        Rout_dids = PROTECT(allocVector(INTSXP, variable_count));
+        Rout_tuples = PROTECT(Rf_allocMatrix(INTSXP, mdfs_info.dimensions, variable_count));
+        Rout_dids = PROTECT(Rf_allocVector(INTSXP, variable_count));
         mdfs_output.setMaxIGsTuples(INTEGER(Rout_tuples), INTEGER(Rout_dids)); // tuples are set row-first during computation, we transpose the result in R to speed up C code
     }
 
-    mdfs[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, contrast_rawdata, std::move(dfi), mdfs_output);
+    mdfs[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, contrast_rawdata, std::move(dfi), mdfs_output);
 
     mdfs_output.copyMaxIGsAsDouble(REAL(Rout_max_igs));
 
@@ -137,12 +136,12 @@ SEXP r_compute_max_ig(
         result_members_count += 1; // for disc nr
     }
 
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         mdfs_output.copyContrastMaxIGsAsDouble(REAL(Rout_contrast_max_igs));
         result_members_count += 1;
     }
 
-    SEXP Rout_result = PROTECT(allocVector(VECSXP, result_members_count));
+    SEXP Rout_result = PROTECT(Rf_allocVector(VECSXP, result_members_count));
     SET_VECTOR_ELT(Rout_result, 0, Rout_max_igs);
 
     if (return_tuples) {
@@ -150,7 +149,7 @@ SEXP r_compute_max_ig(
         SET_VECTOR_ELT(Rout_result, 2, Rout_dids);
     }
 
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         if (return_tuples) {
             SET_VECTOR_ELT(Rout_result, 3, Rout_contrast_max_igs);
         } else {
@@ -160,7 +159,7 @@ SEXP r_compute_max_ig(
 
     UNPROTECT(1 + result_members_count);
 
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         delete contrast_rawdata;
     }
 
@@ -181,71 +180,71 @@ SEXP r_compute_max_ig_discrete(
         SEXP Rin_use_cuda)
 {
     #ifndef WITH_CUDA
-    if (asLogical(Rin_use_cuda)) {
-        error("CUDA acceleration not compiled");
+    if (Rf_asLogical(Rin_use_cuda)) {
+        Rf_error("CUDA acceleration not compiled");
     }
     #endif
 
-    const int* dataDims = INTEGER(getAttrib(Rin_data, R_DimSymbol));
+    const int* dataDims = INTEGER(Rf_getAttrib(Rin_data, R_DimSymbol));
     int* contrastDataDims = nullptr;
-    if (!isNull(Rin_contrast_data)) {
-        contrastDataDims = INTEGER(getAttrib(Rin_contrast_data, R_DimSymbol));
+    if (!Rf_isNull(Rin_contrast_data)) {
+        contrastDataDims = INTEGER(Rf_getAttrib(Rin_contrast_data, R_DimSymbol));
     }
 
     const int obj_count = dataDims[0];
     const int variable_count = dataDims[1];
     int contrast_variable_count = 0;
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         contrast_variable_count = contrastDataDims[1];
     }
 
     #ifdef WITH_CUDA
-    if (asLogical(Rin_use_cuda)) {
-        error("CUDA not supported yet for the discrete variant");
+    if (Rf_asLogical(Rin_use_cuda)) {
+        Rf_error("CUDA not supported yet for the discrete variant");
     }
     #endif
 
-    const int divisions = asInteger(Rin_divisions);
+    const int divisions = Rf_asInteger(Rin_divisions);
 
     const int* decision = INTEGER(Rin_decision);
 
     RawData rawdata(RawDataInfo(obj_count, variable_count), INTEGER(Rin_data), decision);
     RawData* contrast_rawdata = nullptr;
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         contrast_rawdata = new RawData(RawDataInfo(obj_count, contrast_variable_count), INTEGER(Rin_contrast_data), nullptr);
     }
 
     MDFSInfo mdfs_info(
-        asInteger(Rin_dimensions),
+        Rf_asInteger(Rin_dimensions),
         divisions,
         1, // only one discretization
-        asReal(Rin_pseudocount),
+        Rf_asReal(Rin_pseudocount),
         0.0f,
         INTEGER(Rin_interesting_vars),
-        length(Rin_interesting_vars),
-        asLogical(Rin_require_all_vars),
+        Rf_length(Rin_interesting_vars),
+        Rf_asLogical(Rin_require_all_vars),
         nullptr,
         false
     );
 
-    SEXP Rout_max_igs = PROTECT(allocVector(REALSXP, variable_count));
+    SEXP Rout_max_igs = PROTECT(Rf_allocVector(REALSXP, variable_count));
     SEXP Rout_contrast_max_igs = nullptr;
     SEXP Rout_tuples = nullptr;
     SEXP Rout_dids = nullptr;
 
-    if (!isNull(Rin_contrast_data)) {
-        Rout_contrast_max_igs = PROTECT(allocVector(REALSXP, contrast_variable_count));
+    if (!Rf_isNull(Rin_contrast_data)) {
+        Rout_contrast_max_igs = PROTECT(Rf_allocVector(REALSXP, contrast_variable_count));
     }
 
-    const bool return_tuples = asLogical(Rin_return_tuples);
+    const bool return_tuples = Rf_asLogical(Rin_return_tuples);
     MDFSOutput mdfs_output(MDFSOutputType::MaxIGs, mdfs_info.dimensions, variable_count, contrast_variable_count);
     if (return_tuples) {
-        Rout_tuples = PROTECT(allocMatrix(INTSXP, mdfs_info.dimensions, variable_count));
-        Rout_dids = PROTECT(allocVector(INTSXP, variable_count));
+        Rout_tuples = PROTECT(Rf_allocMatrix(INTSXP, mdfs_info.dimensions, variable_count));
+        Rout_dids = PROTECT(Rf_allocVector(INTSXP, variable_count));
         mdfs_output.setMaxIGsTuples(INTEGER(Rout_tuples), INTEGER(Rout_dids)); // tuples are set row-first during computation, we transpose the result in R to speed up C code
     }
 
-    mdfs[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, contrast_rawdata, nullptr, mdfs_output);
+    mdfs[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, contrast_rawdata, nullptr, mdfs_output);
 
     mdfs_output.copyMaxIGsAsDouble(REAL(Rout_max_igs));
 
@@ -256,12 +255,12 @@ SEXP r_compute_max_ig_discrete(
         result_members_count += 1; // for disc nr
     }
 
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         mdfs_output.copyContrastMaxIGsAsDouble(REAL(Rout_contrast_max_igs));
         result_members_count += 1;
     }
 
-    SEXP Rout_result = PROTECT(allocVector(VECSXP, result_members_count));
+    SEXP Rout_result = PROTECT(Rf_allocVector(VECSXP, result_members_count));
     SET_VECTOR_ELT(Rout_result, 0, Rout_max_igs);
 
     if (return_tuples) {
@@ -269,7 +268,7 @@ SEXP r_compute_max_ig_discrete(
         SET_VECTOR_ELT(Rout_result, 2, Rout_dids);
     }
 
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         if (return_tuples) {
             SET_VECTOR_ELT(Rout_result, 3, Rout_contrast_max_igs);
         } else {
@@ -279,7 +278,7 @@ SEXP r_compute_max_ig_discrete(
 
     UNPROTECT(1 + result_members_count);
 
-    if (!isNull(Rin_contrast_data)) {
+    if (!Rf_isNull(Rin_contrast_data)) {
         delete contrast_rawdata;
     }
 
@@ -304,73 +303,73 @@ SEXP r_compute_all_matching_tuples(
         SEXP Rin_stat_mode,
         SEXP Rin_average)
 {
-    const int* dataDims = INTEGER(getAttrib(Rin_data, R_DimSymbol));
+    const int* dataDims = INTEGER(Rf_getAttrib(Rin_data, R_DimSymbol));
 
     const int obj_count = dataDims[0];
     const int variable_count = dataDims[1];
 
-    const int discretizations = asInteger(Rin_discretizations);
-    const int divisions = asInteger(Rin_divisions);
+    const int discretizations = Rf_asInteger(Rin_discretizations);
+    const int divisions = Rf_asInteger(Rin_divisions);
 
     const int* decision = nullptr;
-    if (!isNull(Rin_decision)) {
+    if (!Rf_isNull(Rin_decision)) {
         decision = INTEGER(Rin_decision);
     }
 
     RawData rawdata(RawDataInfo(obj_count, variable_count), REAL(Rin_data), decision);
 
     std::unique_ptr<const DiscretizationInfo> dfi(new DiscretizationInfo(
-        asInteger(Rin_seed),
+        Rf_asInteger(Rin_seed),
         discretizations,
         divisions,
-        asReal(Rin_range)
+        Rf_asReal(Rin_range)
     ));
 
     const double* I_lower = nullptr;
-    if (!isNull(Rin_I_lower)) {
+    if (!Rf_isNull(Rin_I_lower)) {
         I_lower = REAL(Rin_I_lower);
     }
 
     MDFSInfo mdfs_info(
-        asInteger(Rin_dimensions),
+        Rf_asInteger(Rin_dimensions),
         divisions,
         discretizations,
-        asReal(Rin_pseudocount),
-        asReal(Rin_ig_thr),
+        Rf_asReal(Rin_pseudocount),
+        Rf_asReal(Rin_ig_thr),
         INTEGER(Rin_interesting_vars),
-        length(Rin_interesting_vars),
-        asLogical(Rin_require_all_vars),
+        Rf_length(Rin_interesting_vars),
+        Rf_asLogical(Rin_require_all_vars),
         I_lower,
-        asLogical(Rin_average)
+        Rf_asLogical(Rin_average)
     );
 
-    MDFSOutputType out_type = mdfs_info.dimensions == 2 && asReal(Rin_ig_thr) <= 0.0 && length(Rin_interesting_vars) == 0 ? MDFSOutputType::AllTuples : MDFSOutputType::MatchingTuples;
+    MDFSOutputType out_type = mdfs_info.dimensions == 2 && Rf_asReal(Rin_ig_thr) <= 0.0 && Rf_length(Rin_interesting_vars) == 0 ? MDFSOutputType::AllTuples : MDFSOutputType::MatchingTuples;
     MDFSOutput mdfs_output(out_type, mdfs_info.dimensions, variable_count, 0);
 
-    if (isNull(Rin_decision)) {
-        switch (asInteger(Rin_stat_mode)) {
-            case 1: mdfsEntropy[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
+    if (Rf_isNull(Rin_decision)) {
+        switch (Rf_asInteger(Rin_stat_mode)) {
+            case 1: mdfsEntropy[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
             break;
-            case 2: mdfsMutualInformation[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
+            case 2: mdfsMutualInformation[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
             break;
-            case 3: mdfsVariationOfInformation[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
+            case 3: mdfsVariationOfInformation[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
             break;
-            default: error("Unknown statistic");
+            default: Rf_error("Unknown statistic");
         }
     } else {
-        switch (asInteger(Rin_stat_mode)) {
-            case 1: mdfsDecisionConditionalEntropy[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
+        switch (Rf_asInteger(Rin_stat_mode)) {
+            case 1: mdfsDecisionConditionalEntropy[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
             break;
-            case 2: mdfs[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
+            case 2: mdfs[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
             break;
-            case 3: mdfsDecisionConditionalVariationOfInformation[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
+            case 3: mdfsDecisionConditionalVariationOfInformation[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, std::move(dfi), mdfs_output);
             break;
-            default: error("Unknown statistic");
+            default: Rf_error("Unknown statistic");
         }
     }
 
-    if (out_type == MDFSOutputType::AllTuples && asLogical(Rin_return_matrix)) {
-        SEXP Rout_result = PROTECT(allocMatrix(REALSXP, variable_count, variable_count));
+    if (out_type == MDFSOutputType::AllTuples && Rf_asLogical(Rin_return_matrix)) {
+        SEXP Rout_result = PROTECT(Rf_allocMatrix(REALSXP, variable_count, variable_count));
 
         // TODO: perhaps we could avoid copying here at all and fill in this matrix already from the mdfs?
         mdfs_output.copyAllTuplesMatrix(REAL(Rout_result));
@@ -383,9 +382,9 @@ SEXP r_compute_all_matching_tuples(
         // 2D only now
         const int tuples_count = out_type == MDFSOutputType::AllTuples ? variable_count * (variable_count - 1) : mdfs_output.getMatchingTuplesCount();
 
-        SEXP Rout_igs = PROTECT(allocVector(REALSXP, tuples_count));
-        SEXP Rout_tuples = PROTECT(allocMatrix(INTSXP, tuples_count, mdfs_info.dimensions));
-        SEXP Rout_vars = PROTECT(allocVector(INTSXP, tuples_count));
+        SEXP Rout_igs = PROTECT(Rf_allocVector(REALSXP, tuples_count));
+        SEXP Rout_tuples = PROTECT(Rf_allocMatrix(INTSXP, tuples_count, mdfs_info.dimensions));
+        SEXP Rout_vars = PROTECT(Rf_allocVector(INTSXP, tuples_count));
 
         if (out_type == MDFSOutputType::AllTuples) {
             mdfs_output.copyAllTuples(INTEGER(Rout_vars), REAL(Rout_igs), INTEGER(Rout_tuples));
@@ -393,7 +392,7 @@ SEXP r_compute_all_matching_tuples(
             mdfs_output.copyMatchingTuples(INTEGER(Rout_vars), REAL(Rout_igs), INTEGER(Rout_tuples));
         }
 
-        SEXP Rout_result = PROTECT(allocVector(VECSXP, result_members_count));
+        SEXP Rout_result = PROTECT(Rf_allocVector(VECSXP, result_members_count));
         SET_VECTOR_ELT(Rout_result, 0, Rout_vars);
         SET_VECTOR_ELT(Rout_result, 1, Rout_tuples);
         SET_VECTOR_ELT(Rout_result, 2, Rout_igs);
@@ -418,65 +417,65 @@ SEXP r_compute_all_matching_tuples_discrete(
         SEXP Rin_return_matrix,
         SEXP Rin_stat_mode)
 {
-    const int* dataDims = INTEGER(getAttrib(Rin_data, R_DimSymbol));
+    const int* dataDims = INTEGER(Rf_getAttrib(Rin_data, R_DimSymbol));
 
     const int obj_count = dataDims[0];
     const int variable_count = dataDims[1];
 
-    const int divisions = asInteger(Rin_divisions);
+    const int divisions = Rf_asInteger(Rin_divisions);
 
     const int* decision = nullptr;
-    if (!isNull(Rin_decision)) {
+    if (!Rf_isNull(Rin_decision)) {
         decision = INTEGER(Rin_decision);
     }
 
     RawData rawdata(RawDataInfo(obj_count, variable_count), INTEGER(Rin_data), decision);
 
     const double* I_lower = nullptr;
-    if (!isNull(Rin_I_lower)) {
+    if (!Rf_isNull(Rin_I_lower)) {
         I_lower = REAL(Rin_I_lower);
     }
 
     MDFSInfo mdfs_info(
-        asInteger(Rin_dimensions),
+        Rf_asInteger(Rin_dimensions),
         divisions,
         1, // only one discretization
-        asReal(Rin_pseudocount),
-        asReal(Rin_ig_thr),
+        Rf_asReal(Rin_pseudocount),
+        Rf_asReal(Rin_ig_thr),
         INTEGER(Rin_interesting_vars),
-        length(Rin_interesting_vars),
-        asLogical(Rin_require_all_vars),
+        Rf_length(Rin_interesting_vars),
+        Rf_asLogical(Rin_require_all_vars),
         I_lower,
         false
     );
 
-    MDFSOutputType out_type = mdfs_info.dimensions == 2 && asReal(Rin_ig_thr) <= 0.0 && length(Rin_interesting_vars) == 0 ? MDFSOutputType::AllTuples : MDFSOutputType::MatchingTuples;
+    MDFSOutputType out_type = mdfs_info.dimensions == 2 && Rf_asReal(Rin_ig_thr) <= 0.0 && Rf_length(Rin_interesting_vars) == 0 ? MDFSOutputType::AllTuples : MDFSOutputType::MatchingTuples;
     MDFSOutput mdfs_output(out_type, mdfs_info.dimensions, variable_count, 0);
 
-    if (isNull(Rin_decision)) {
-        switch (asInteger(Rin_stat_mode)) {
-            case 1: mdfsEntropy[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
+    if (Rf_isNull(Rin_decision)) {
+        switch (Rf_asInteger(Rin_stat_mode)) {
+            case 1: mdfsEntropy[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
             break;
-            case 2: mdfsMutualInformation[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
+            case 2: mdfsMutualInformation[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
             break;
-            case 3: mdfsVariationOfInformation[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
+            case 3: mdfsVariationOfInformation[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
             break;
-            default: error("Unknown statistic");
+            default: Rf_error("Unknown statistic");
         }
     } else {
-        switch (asInteger(Rin_stat_mode)) {
-            case 1: mdfsDecisionConditionalEntropy[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
+        switch (Rf_asInteger(Rin_stat_mode)) {
+            case 1: mdfsDecisionConditionalEntropy[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
             break;
-            case 2: mdfs[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
+            case 2: mdfs[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
             break;
-            case 3: mdfsDecisionConditionalVariationOfInformation[asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
+            case 3: mdfsDecisionConditionalVariationOfInformation[Rf_asInteger(Rin_dimensions)-1](mdfs_info, &rawdata, nullptr, nullptr, mdfs_output);
             break;
-            default: error("Unknown statistic");
+            default: Rf_error("Unknown statistic");
         }
     }
 
-    if (out_type == MDFSOutputType::AllTuples && asLogical(Rin_return_matrix)) {
-        SEXP Rout_result = PROTECT(allocMatrix(REALSXP, variable_count, variable_count));
+    if (out_type == MDFSOutputType::AllTuples && Rf_asLogical(Rin_return_matrix)) {
+        SEXP Rout_result = PROTECT(Rf_allocMatrix(REALSXP, variable_count, variable_count));
 
         // TODO: perhaps we could avoid copying here at all and fill in this matrix already from the mdfs?
         mdfs_output.copyAllTuplesMatrix(REAL(Rout_result));
@@ -489,9 +488,9 @@ SEXP r_compute_all_matching_tuples_discrete(
         // 2D only now
         const int tuples_count = out_type == MDFSOutputType::AllTuples ? variable_count * (variable_count - 1) : mdfs_output.getMatchingTuplesCount();
 
-        SEXP Rout_igs = PROTECT(allocVector(REALSXP, tuples_count));
-        SEXP Rout_tuples = PROTECT(allocMatrix(INTSXP, tuples_count, mdfs_info.dimensions));
-        SEXP Rout_vars = PROTECT(allocVector(INTSXP, tuples_count));
+        SEXP Rout_igs = PROTECT(Rf_allocVector(REALSXP, tuples_count));
+        SEXP Rout_tuples = PROTECT(Rf_allocMatrix(INTSXP, tuples_count, mdfs_info.dimensions));
+        SEXP Rout_vars = PROTECT(Rf_allocVector(INTSXP, tuples_count));
 
         if (out_type == MDFSOutputType::AllTuples) {
             mdfs_output.copyAllTuples(INTEGER(Rout_vars), REAL(Rout_igs), INTEGER(Rout_tuples));
@@ -499,7 +498,7 @@ SEXP r_compute_all_matching_tuples_discrete(
             mdfs_output.copyMatchingTuples(INTEGER(Rout_vars), REAL(Rout_igs), INTEGER(Rout_tuples));
         }
 
-        SEXP Rout_result = PROTECT(allocVector(VECSXP, result_members_count));
+        SEXP Rout_result = PROTECT(Rf_allocVector(VECSXP, result_members_count));
         SET_VECTOR_ELT(Rout_result, 0, Rout_vars);
         SET_VECTOR_ELT(Rout_result, 1, Rout_tuples);
         SET_VECTOR_ELT(Rout_result, 2, Rout_igs);
@@ -519,12 +518,12 @@ SEXP r_discretize(
         SEXP Rin_seed,
         SEXP Rin_range)
 {
-    const R_len_t obj_count = length(Rin_variable);
-    const int discretization_nr = asInteger(Rin_discretization_nr);
-    const int variable_idx = asInteger(Rin_variable_idx);
-    const int divisions = asInteger(Rin_divisions);
-    const int seed = asInteger(Rin_seed);
-    double range = asReal(Rin_range);
+    const R_len_t obj_count = Rf_length(Rin_variable);
+    const int discretization_nr = Rf_asInteger(Rin_discretization_nr);
+    const int variable_idx = Rf_asInteger(Rin_variable_idx);
+    const int divisions = Rf_asInteger(Rin_divisions);
+    const int seed = Rf_asInteger(Rin_seed);
+    double range = Rf_asReal(Rin_range);
     double* variable = REAL(Rin_variable);
 
     std::vector<double> sorted_variable(variable, variable + obj_count);
@@ -532,7 +531,7 @@ SEXP r_discretize(
 
     uint8_t* discretized_variable = new uint8_t[obj_count];
     discretize(seed, discretization_nr, variable_idx, divisions, obj_count, variable, sorted_variable, discretized_variable, range);
-    SEXP Rout_result = PROTECT(allocVector(INTSXP, obj_count));
+    SEXP Rout_result = PROTECT(Rf_allocVector(INTSXP, obj_count));
     std::copy(discretized_variable, discretized_variable + obj_count, INTEGER(Rout_result));
     delete[] discretized_variable;
     UNPROTECT(1);
@@ -545,7 +544,7 @@ SEXP r_omp_set_num_threads(
         SEXP Rin_num_threads)
 {
     #ifdef _OPENMP
-    const int num_threads = asInteger(Rin_num_threads);
+    const int num_threads = Rf_asInteger(Rin_num_threads);
     omp_set_num_threads(num_threads);
     #endif
     return R_NilValue;
diff --git a/src/r_interface.h b/src/r_interface.h
index 4bace68..827f65d 100644
--- a/src/r_interface.h
+++ b/src/r_interface.h
@@ -1,3 +1,7 @@
+// to make CRAN happy as it will be the default in R 4.5.0
+#ifndef R_NO_REMAP
+# define R_NO_REMAP
+#endif
 #include <Rinternals.h>
 
 extern "C"