Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision b966a96fbdc2825411d2c2a16c6f9e47f38e9ba5)
+++ sasmodels/kernel_header.c	(revision 1557a1ec73531576042f730e22cd319338baafc3)
@@ -7,5 +7,15 @@
 // If opencl is not available, then we are compiling a C function
 // Note: if using a C++ compiler, then define kernel as extern "C"
-#ifndef USE_OPENCL
+#ifdef USE_OPENCL
+   typedef int int32_t;
+#  if defined(USE_SINCOS)
+#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
+#  else
+#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
+#  endif
+   // OpenCL only has type generic math
+   #define erff erf
+   #define erfcf erfc
+#else // !USE_OPENCL
 // Use SAS_DOUBLE to force the use of double even for float kernels
 #  define SAS_DOUBLE dou ## ble
@@ -26,4 +36,5 @@
          #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN
          #define INFINITY (std::numeric_limits<double>::infinity())
+         #define NEED_ERF
          #define NEED_EXPM1
          #define NEED_TGAMMA
@@ -33,5 +44,5 @@
      #endif
      inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); }
-#  else
+#  else // !__cplusplus
      #include <inttypes.h>  // C99 guarantees that int32_t types is here
      #include <stdio.h>
@@ -47,4 +58,5 @@
          inline SAS_DOUBLE fmin(SAS_DOUBLE x, SAS_DOUBLE y) { return x>y ? y : x; }
          inline SAS_DOUBLE fmax(SAS_DOUBLE x, SAS_DOUBLE y) { return x<y ? y : x; }
+         #define NEED_ERF
          #define NEED_EXPM1
          #define NEED_TGAMMA
@@ -55,5 +67,5 @@
      #define kernel
      #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
-#  endif
+#  endif  // !__cplusplus
 #  define global
 #  define local
@@ -63,12 +75,5 @@
 #  define powr(a,b) pow(a,b)
 #  define pown(a,b) pow(a,b)
-#else
-   typedef int int32_t;
-#  if defined(USE_SINCOS)
-#    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar)
-#  else
-#    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0)
-#  endif
-#endif
+#endif // !USE_OPENCL
 
 #if defined(NEED_EXPM1)
Index: sasmodels/models/lib/sas_erf.c
===================================================================
--- sasmodels/models/lib/sas_erf.c	(revision 1557a1ec73531576042f730e22cd319338baafc3)
+++ sasmodels/models/lib/sas_erf.c	(revision 1557a1ec73531576042f730e22cd319338baafc3)
@@ -0,0 +1,352 @@
+
+/*
+ * Cephes Math Library Release 2.2:  June, 1992
+ * Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ */
+
+/*
+ *
+ *	Error function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, erf();
+ *
+ * y = erf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * The integral is
+ *
+ *                           x
+ *                            -
+ *                 2         | |          2
+ *   erf(x)  =  --------     |    exp( - t  ) dt.
+ *              sqrt(pi)   | |
+ *                          -
+ *                           0
+ *
+ * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
+ * erf(x) = 1 - erfc(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,1         30000       3.7e-16     1.0e-16
+ *
+ */
+
+/*
+ *
+ *	Complementary error function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, erfc();
+ *
+ * y = erfc( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *
+ *  1 - erf(x) =
+ *
+ *                           inf.
+ *                             -
+ *                  2         | |          2
+ *   erfc(x)  =  --------     |    exp( - t  ) dt
+ *               sqrt(pi)   | |
+ *                           -
+ *                            x
+ *
+ *
+ * For small x, erfc(x) = 1 - erf(x); otherwise rational
+ * approximations are computed.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,26.6417   30000       5.7e-14     1.5e-14
+ */
+
+
+#ifdef NEED_ERF
+
+#if FLOAT_SIZE>4  // DOUBLE_PRECISION
+double erf(double x);
+double erfc(double a);
+
+constant double PD[] = {
+    2.46196981473530512524E-10,
+    5.64189564831068821977E-1,
+    7.46321056442269912687E0,
+    4.86371970985681366614E1,
+    1.96520832956077098242E2,
+    5.26445194995477358631E2,
+    9.34528527171957607540E2,
+    1.02755188689515710272E3,
+    5.57535335369399327526E2
+};
+
+constant double QD[] = {
+    /* 1.00000000000000000000E0, */
+    1.32281951154744992508E1,
+    8.67072140885989742329E1,
+    3.54937778887819891062E2,
+    9.75708501743205489753E2,
+    1.82390916687909736289E3,
+    2.24633760818710981792E3,
+    1.65666309194161350182E3,
+    5.57535340817727675546E2
+};
+
+constant double RD[] = {
+    5.64189583547755073984E-1,
+    1.27536670759978104416E0,
+    5.01905042251180477414E0,
+    6.16021097993053585195E0,
+    7.40974269950448939160E0,
+    2.97886665372100240670E0
+};
+
+constant double SD[] = {
+    /* 1.00000000000000000000E0, */
+    2.26052863220117276590E0,
+    9.39603524938001434673E0,
+    1.20489539808096656605E1,
+    1.70814450747565897222E1,
+    9.60896809063285878198E0,
+    3.36907645100081516050E0
+};
+
+constant double TD[] = {
+    9.60497373987051638749E0,
+    9.00260197203842689217E1,
+    2.23200534594684319226E3,
+    7.00332514112805075473E3,
+    5.55923013010394962768E4
+};
+
+constant double UD[] = {
+    /* 1.00000000000000000000E0, */
+    3.35617141647503099647E1,
+    5.21357949780152679795E2,
+    4.59432382970980127987E3,
+    2.26290000613890934246E4,
+    4.92673942608635921086E4
+};
+
+double erfc(double a)
+{
+    double MAXLOG = 88.72283905206835;
+    double p, q, x, y, z;
+
+
+    /*if (a < 0.0)
+        x = -a;
+    else
+        x = a;*/
+
+    x = fabs(a);
+
+
+    if (x < 1.0) {
+        //The line bellow is a troublemaker for GPU, so sas_erf function
+        //is explicit here for the case < 1.0
+        //return (1.0 - sas_erf(a));
+        z = x * x;
+        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
+
+        return y;
+    }
+
+    z = -a * a;
+
+    if (z < -MAXLOG) {
+        if (a < 0)
+            return (2.0);
+        else
+            return (0.0);
+    }
+
+    z = exp(z);
+
+
+    if (x < 8.0) {
+        p = polevl(x, PD, 8);
+        q = p1evl(x, QD, 8);
+    }
+    else {
+        p = polevl(x, RD, 5);
+        q = p1evl(x, SD, 6);
+    }
+    y = (z * p) / q;
+
+    if (a < 0)
+        y = 2.0 - y;
+
+    if (y == 0.0) {
+        if (a < 0)
+            return (2.0);
+        else
+            return (0.0);
+    }
+    return y;
+}
+
+
+double erf(double x)
+{
+    double y, z;
+
+    if (fabs(x) > 1.0)
+        return (1.0 - erfc(x));
+
+    z = x * x;
+    #if FLOAT_SIZE>4
+        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
+    #else
+        y = x * polevl( z, TF, 6 );
+    #endif
+
+    return y;
+}
+
+#else // SINGLE PRECISION
+
+double erff(double x);
+double erfcf(double a);
+
+/* erfc(x) = exp(-x^2) P(1/x), 1 < x < 2 */
+constant double PF[] = {
+    2.326819970068386E-002,
+    -1.387039388740657E-001,
+    3.687424674597105E-001,
+    -5.824733027278666E-001,
+    6.210004621745983E-001,
+    -4.944515323274145E-001,
+    3.404879937665872E-001,
+    -2.741127028184656E-001,
+    5.638259427386472E-001
+};
+
+/* erfc(x) = exp(-x^2) 1/x P(1/x^2), 2 < x < 14 */
+constant double RF[] = {
+    -1.047766399936249E+001,
+    1.297719955372516E+001,
+    -7.495518717768503E+000,
+    2.921019019210786E+000,
+    -1.015265279202700E+000,
+    4.218463358204948E-001,
+    -2.820767439740514E-001,
+    5.641895067754075E-001
+};
+
+/* erf(x) = x P(x^2), 0 < x < 1 */
+ constant double TF[] = {
+    7.853861353153693E-005,
+    -8.010193625184903E-004,
+    5.188327685732524E-003,
+    -2.685381193529856E-002,
+    1.128358514861418E-001,
+    -3.761262582423300E-001,
+    1.128379165726710E+000
+};
+
+
+float erfcf(float a)
+{
+    float MAXLOG = 88.72283905206835;
+    float p, q, x, y, z;
+
+
+    /*if (a < 0.0)
+        x = -a;
+    else
+        x = a;*/
+
+    x = fabsf(a);
+
+
+    if (x < 1.0) {
+        //The line below is a troublemaker for GPU, so sas_erf function
+        //is explicit here for the case < 1.0
+        //return (1.0 - sas_erf(a));
+        z = x * x;
+        y = x * polevl( z, TF, 6 );
+
+        return y;
+    }
+
+    z = -a * a;
+
+    if (z < -MAXLOG) {
+        if (a < 0)
+            return (2.0);
+        else
+            return (0.0);
+    }
+
+    z = expf(z);
+
+
+    q=1.0/x;
+    y=q*q;
+    if( x < 2.0 ) {
+        p = polevl( y, PF, 8 );
+    } else {
+        p = polevl( y, RF, 7 );
+    }
+    y = z * q * p;
+
+    if (a < 0)
+        y = 2.0 - y;
+
+    if (y == 0.0) {
+        if (a < 0)
+            return (2.0);
+        else
+            return (0.0);
+    }
+    return y;
+}
+
+
+float erff(float x)
+{
+    float y, z;
+
+    if (fabsf(x) > 1.0)
+        return (1.0 - erfcf(x));
+
+    z = x * x;
+    y = x * polevl( z, TF, 6 );
+
+    return y;
+}
+
+#endif // SINGLE_PRECISION
+#endif // NEED_ERF
+
+#if FLOAT_SIZE>4
+#define sas_erf erf
+#define sas_erfc erfc
+#else
+#define sas_erf erff
+#define sas_erfc erfcf
+#endif
