diff --git a/NEWS.md b/NEWS.md
index e50b0b2..007f7e2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,8 @@
   the `rxode2` internals use the function pointers, while the models
   built by `rxode2` use the binary linkage.
 
+* As requested by `CRAN` use strict headers.
+
 # PreciseSums 0.6
 
 * Updated linkage for C.
diff --git a/src/Makevars b/src/Makevars
index 39fa148..e3b80c0 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,2 +1,3 @@
+PKG_CFLAGS  = -DSTRICT_R_HEADERS=1
 SOURCES_C = sum.c prod.c init.c
-OBJECTS = $(SOURCES_C:.c=.o)
\ No newline at end of file
+OBJECTS = $(SOURCES_C:.c=.o)
diff --git a/src/PreciseSums.h b/src/PreciseSums.h
index c4458e5..c16c169 100644
--- a/src/PreciseSums.h
+++ b/src/PreciseSums.h
@@ -11,6 +11,7 @@
 #define PS_Prod 2
 #define PS_LogifyProd 3
 
+#define R_NO_REMAP
 #include <R.h>
 #include <Rinternals.h>
 #include <Rdefines.h>
diff --git a/src/PreciseSumsPtr.h b/src/PreciseSumsPtr.h
index 56da7d1..689d10b 100644
--- a/src/PreciseSumsPtr.h
+++ b/src/PreciseSumsPtr.h
@@ -11,6 +11,7 @@
 #define PS_Prod 2
 #define PS_LogifyProd 3
 
+#define R_NO_REMAP
 #include <R.h>
 #include <Rinternals.h>
 #include <Rdefines.h>
diff --git a/src/init.c b/src/init.c
index 9913c47..b3627c0 100644
--- a/src/init.c
+++ b/src/init.c
@@ -1,5 +1,6 @@
 #include <stdio.h>
 #include <stdlib.h>
+#define R_NO_REMAP
 #include <R.h>
 #include <Rinternals.h>
 #include <Rmath.h> //Rmath includes math.
@@ -40,7 +41,7 @@ SEXP _PreciseSumsPtr(void) {
   SEXP PreciseSums_sum_get_ptr = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&PreciseSums_sum_get, R_NilValue, R_NilValue)); pro++;
   SEXP PreciseSums_prod_get_ptr = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&PreciseSums_prod_get, R_NilValue, R_NilValue)); pro++;
 
-  SEXP ret = PROTECT(allocVector(VECSXP, 6)); pro++;
+  SEXP ret = PROTECT(Rf_allocVector(VECSXP, 6)); pro++;
   SET_VECTOR_ELT(ret, 0, PreciseSums_sum_ptr);
   SET_VECTOR_ELT(ret, 1, PreciseSums_prod_ptr);
   SET_VECTOR_ELT(ret, 2, PreciseSums_sum_r_ptr);
@@ -48,15 +49,15 @@ SEXP _PreciseSumsPtr(void) {
   SET_VECTOR_ELT(ret, 4, PreciseSums_sum_get_ptr);
   SET_VECTOR_ELT(ret, 5, PreciseSums_prod_get_ptr);
 
-  SEXP retn = PROTECT(allocVector(STRSXP, 6)); pro++;
-  SET_STRING_ELT(retn, 0, mkChar("PreciseSums_sum"));
-  SET_STRING_ELT(retn, 1, mkChar("PreciseSums_prod"));
-  SET_STRING_ELT(retn, 2, mkChar("PreciseSums_sum_r"));
-  SET_STRING_ELT(retn, 3, mkChar("PreciseSums_prod_r"));
-  SET_STRING_ELT(retn, 4, mkChar("PreciseSums_sum_get"));
-  SET_STRING_ELT(retn, 5, mkChar("PreciseSums_prod_get"));
+  SEXP retn = PROTECT(Rf_allocVector(STRSXP, 6)); pro++;
+  SET_STRING_ELT(retn, 0, Rf_mkChar("PreciseSums_sum"));
+  SET_STRING_ELT(retn, 1, Rf_mkChar("PreciseSums_prod"));
+  SET_STRING_ELT(retn, 2, Rf_mkChar("PreciseSums_sum_r"));
+  SET_STRING_ELT(retn, 3, Rf_mkChar("PreciseSums_prod_r"));
+  SET_STRING_ELT(retn, 4, Rf_mkChar("PreciseSums_sum_get"));
+  SET_STRING_ELT(retn, 5, Rf_mkChar("PreciseSums_prod_get"));
 
-  setAttrib(ret, R_NamesSymbol, retn);
+  Rf_setAttrib(ret, R_NamesSymbol, retn);
 
   UNPROTECT(pro);
   return ret;
diff --git a/src/prod.c b/src/prod.c
index 88a5d35..6a3e13b 100644
--- a/src/prod.c
+++ b/src/prod.c
@@ -1,10 +1,14 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdarg.h>
+#define R_NO_REMAP
 #include <R.h>
 #include <Rinternals.h>
 #include <Rmath.h> //Rmath includes math.
 
+// https://github.com/wch/r-source/blob/41ca75a4ffa1a4ec66f7697820fcf229b6a64f8e/src/nmath/qgamma.c#L122
+#define LN_EPS -36.043653389117156
+
 extern double PreciseSums_sum (double *input, int n);
 
 extern double PreciseSums_pairwise_add_DOUBLE(double *a, int n);
@@ -12,10 +16,11 @@ extern double PreciseSums_pairwise_add_DOUBLE(double *a, int n);
 double PreciseSums_safe_log(double x){
   if (x <= 0){
     // Warning?
-    return log(DBL_EPSILON);
+    return LN_EPS;
   } else {
     return log(x);
   }
+  return LN_EPS;
 }
 
 int PreciseSums_prod_type = 3;
@@ -38,7 +43,7 @@ extern double PreciseSums_prod_ld(double *input, int n){
   long double p = 1;
   for  (int i = 0; i < n; i++){
     if (input[i] == 0){
-      return 0.0; 
+      return 0.0;
     }
     p *= input[i];
   }
@@ -49,7 +54,7 @@ extern double PreciseSums_prod_d(double *input, int n){
   double p = 1;
   for  (int i = 0; i < n; i++){
     if (input[i] == 0){
-      return 0.0; 
+      return 0.0;
     }
     p *= input[i];
   }
@@ -70,9 +75,9 @@ extern double PreciseSums_prod_logify_r(double *input, double *p, int n){
 }
 
 extern double PreciseSums_prod_logify(double *input, int n){
-  double *p = Calloc(n,double);
+  double *p = R_Calloc(n,double);
   double ret= PreciseSums_prod_logify_r(input, p, n);
-  Free(p);
+  R_Free(p);
   return ret;
 }
 
@@ -110,9 +115,9 @@ extern double PreciseSums_prod(double *input, int n){
 
 
 SEXP _psProd(SEXP input){
-  int len = length(input);
+  int len = Rf_length(input);
   double *dinput = REAL(input);
-  SEXP rets = PROTECT(allocVector(REALSXP,1));
+  SEXP rets = PROTECT(Rf_allocVector(REALSXP,1));
   REAL(rets)[0] = PreciseSums_prod(dinput, len);
   UNPROTECT(1);
   return rets;
@@ -121,13 +126,13 @@ SEXP _psProd(SEXP input){
 extern double PreciseSums_prodV(int n, ...){
   va_list valist;
   va_start(valist, n);
-  double *p = Calloc(n, double);
+  double *p = R_Calloc(n, double);
   for (int i = 0; i < n; i++){
     p[i] = va_arg(valist, double);
   }
   va_end(valist);
   double s = PreciseSums_prod(p, n);
-  Free(p);
+  R_Free(p);
   return s;
 }
 
diff --git a/src/sum.c b/src/sum.c
index 3503bdb..9d90ef0 100644
--- a/src/sum.c
+++ b/src/sum.c
@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdarg.h>
+#define R_NO_REMAP
 #include <R.h>
 #include <Rinternals.h>
 #include <Rmath.h> //Rmath includes math.
@@ -87,9 +88,9 @@ PreciseSums_pairwise_add_DOUBLE(double *a, int n)
 }
 
 SEXP _psPairwiseSum(SEXP input){
-  int len = length(input);
+  int len = Rf_length(input);
   double *dinput = REAL(input);
-  SEXP rets = PROTECT(allocVector(REALSXP,1));
+  SEXP rets = PROTECT(Rf_allocVector(REALSXP,1));
   REAL(rets)[0] = PreciseSums_pairwise_add_DOUBLE(dinput, len);
   UNPROTECT(1);
   return rets;
@@ -102,7 +103,7 @@ extern double PreciseSums_KahanSum(double *input, int len){
   int i;
 
   for (i = 0; i < len; i++){
-    y = (double)input[i] - c; 
+    y = (double)input[i] - c;
     t = sum + y;
     c = (t - sum) - y;       // (t - sum) cancels the high-order part of y; subtracting y recovers negative (low part of y)
     sum = t;                 // Algebraically, c should always be zero. Beware overly-aggressive optimizing compilers!
@@ -111,9 +112,9 @@ extern double PreciseSums_KahanSum(double *input, int len){
 }
 
 SEXP _psKahanSum(SEXP input){
-  int len = length(input);
+  int len = Rf_length(input);
   double *dinput = REAL(input);
-  SEXP rets = PROTECT(allocVector(REALSXP,1));
+  SEXP rets = PROTECT(Rf_allocVector(REALSXP,1));
   REAL(rets)[0] = PreciseSums_KahanSum(dinput, len);
   UNPROTECT(1);
   return rets;
@@ -159,9 +160,9 @@ extern double PreciseSums_KleinSum(double *input, int len){
 }
 
 SEXP _psNeumaierSum(SEXP input){
-  int len = length(input);
+  int len = Rf_length(input);
   double *dinput = REAL(input);
-  SEXP rets = PROTECT(allocVector(REALSXP,1));
+  SEXP rets = PROTECT(Rf_allocVector(REALSXP,1));
   REAL(rets)[0] = PreciseSums_NeumaierSum(dinput, len);
   UNPROTECT(1);
   return rets;
@@ -194,8 +195,8 @@ extern double PreciseSums_Python_fsum_r(double *iterable, int iterable_len, doub
         p[i++] = lo;
       x = hi;
     }
-    
-    n = i; 
+
+    n = i;
     if (x != 0.0) {
       if (!R_FINITE(x)) {
         /* a nonfinite x could arise either as
@@ -203,7 +204,7 @@ extern double PreciseSums_Python_fsum_r(double *iterable, int iterable_len, doub
            as a result of a nan or inf in the
            summands */
         if (R_FINITE(xsave) || ISNAN(xsave)) {
-	  if (m > 0) Free(p);
+	  if (m > 0) R_Free(p);
 	  return PreciseSums_KleinSum(iterable, iterable_len);
           /* error("intermediate overflow in fsum"); */
         } else {
@@ -217,9 +218,9 @@ extern double PreciseSums_Python_fsum_r(double *iterable, int iterable_len, doub
           //&& _fsum_realloc(&p, n, ps, &m)
           // Doubles the size of array.
           m += m;
-          p = Realloc(p, m, double);
+          p = R_Realloc(p, m, double);
         } else if (m < 0 && n >= -m){
-	  /* if (m > 0) Free(p); */
+	  /* if (m > 0) R_Free(p); */
 	  return PreciseSums_KleinSum(iterable, iterable_len);
 	  /* error("The size of the saved partials is too small to calculate the sum."); */
 	}
@@ -229,8 +230,8 @@ extern double PreciseSums_Python_fsum_r(double *iterable, int iterable_len, doub
   }
   if (special_sum != 0.0) {
     if (ISNAN(inf_sum)){
-      if (m > 0) Free(p);
-      error("-inf + inf in fsum");
+      if (m > 0) R_Free(p);
+      Rf_error("-inf + inf in fsum");
     }
     sum = special_sum;
     return sum;
@@ -246,7 +247,7 @@ extern double PreciseSums_Python_fsum_r(double *iterable, int iterable_len, doub
       x = hi;
       y = p[--n];
       if (fabs(y) >= fabs(x)){
-	if (m > 0) Free(p);
+	if (m > 0) R_Free(p);
 	return PreciseSums_KleinSum(iterable, iterable_len);
         /* Rprintf("Partial Sums:\n"); */
         /* for (i = 0; i < j; i++){ */
@@ -254,7 +255,7 @@ extern double PreciseSums_Python_fsum_r(double *iterable, int iterable_len, doub
         /* } */
         /* Rprintf("Assertion Error:\n"); */
         /* Rprintf("fabs(y) >= fabs(x) or %f >= %f\n",fabs(y),fabs(x)); */
-        /* if (m > 0) Free(p); */
+        /* if (m > 0) R_Free(p); */
         /* error("Error in parital sums."); */
       }
       hi = x + y;
@@ -282,17 +283,17 @@ extern double PreciseSums_Python_fsum_r(double *iterable, int iterable_len, doub
 }
 
 extern double PreciseSums_Python_fsum(double *iterable, int iterable_len){
-  double *p = Calloc(NUM_PARTIALS, double);
+  double *p = R_Calloc(NUM_PARTIALS, double);
   int m = NUM_PARTIALS;
   double ret = PreciseSums_Python_fsum_r(&iterable[0], iterable_len, p, m);
-  Free(p);
+  R_Free(p);
   return ret;
 }
 
 SEXP _psPythonSum(SEXP input){
-  int len = length(input);
+  int len = Rf_length(input);
   double *dinput = REAL(input);
-  SEXP rets = PROTECT(allocVector(REALSXP,1));
+  SEXP rets = PROTECT(Rf_allocVector(REALSXP,1));
   REAL(rets)[0] = PreciseSums_Python_fsum(dinput, len);
   UNPROTECT(1);
   return rets;
@@ -321,7 +322,7 @@ extern double PreciseSums_sum (double *input, int n){
   }
   //PreciseSums_KahanSum;
   // PreciseSums_NeumaierSum
-  error("Unknown sum type.");
+  Rf_error("Unknown sum type.");
   return 0;
 }
 
@@ -347,7 +348,7 @@ extern double PreciseSums_sum_r(double *input, int n, double *p, int m, int type
   }
   //PreciseSums_KahanSum;
   // PreciseSums_NeumaierSum
-  error("Unknown sum type.");
+  Rf_error("Unknown sum type.");
   return 0;
 }
 
@@ -366,9 +367,9 @@ SEXP _psSetSum(SEXP input){
 }
 
 SEXP _psSum(SEXP input){
-  int len = length(input);
+  int len = Rf_length(input);
   double *dinput = REAL(input);
-  SEXP rets = PROTECT(allocVector(REALSXP,1));
+  SEXP rets = PROTECT(Rf_allocVector(REALSXP,1));
   REAL(rets)[0] = PreciseSums_sum(dinput, len);
   UNPROTECT(1);
   return rets;
@@ -377,13 +378,13 @@ SEXP _psSum(SEXP input){
 extern double PreciseSums_sumV(int n, ...){
   va_list valist;
   va_start(valist, n);
-  double *p = Calloc(n, double);
+  double *p = R_Calloc(n, double);
   for (unsigned int i = 0; i < n; i++){
     p[i] = va_arg(valist, double);
   }
   va_end(valist);
   double s = PreciseSums_sum(p, n);
-  Free(p);
+  R_Free(p);
   return s;
 }