erfcx: "Compilation ERROR, function(s)/method(s) not created!"

I’m trying to make a model with an exponentially modified Gaussian, with a scale parameter (not rate, as in the built-in version) that can vary between positive, zero, or negative.

In order to do this, I tried to implement the erfcx function, which seems to be missing from stan. Based on the example code at https://stackoverflow.com/questions/39777360/accurate-computation-of-scaled-complementary-error-function-erfcx , I’ve written the following functions:

functions {
  real fmaf(real x, real y, real z) {
    return x * y + z;
  }
  real erfcx(real x) {
     real a;
     real d;
     real e;
     real m;
     real p;
     real q;
     real r;
     real s;
     real t;

       a = abs(x); // NaN-preserving absolute value computation

       /* Compute q = (a-2)/(a+2) accurately. [0,INF) -> [-1,1] */
       m = a - 2.0;
       p = a + 2.0;
       r = 1.0 / p;
       q = m * r;
       t = fmaf (q + 1.0, -2.0, a);
       e = fmaf (q, -a, t);
       q = fmaf (r, e, q);

       /* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1,1]
       p =              0x1.f10000p-15f;  //  5.92470169e-5
       p = fmaf (p, q,  0x1.521cc6p-13f); //  1.61224554e-4
       p = fmaf (p, q, -0x1.6b4ffep-12f); // -3.46481771e-4
       p = fmaf (p, q, -0x1.6e2a7cp-10f); // -1.39681227e-3
       p = fmaf (p, q,  0x1.3c1d7ep-10f); //  1.20588380e-3
       p = fmaf (p, q,  0x1.1cc236p-07f); //  8.69014394e-3
       p = fmaf (p, q, -0x1.069940p-07f); // -8.01387429e-3
       p = fmaf (p, q, -0x1.bc1b6cp-05f); // -5.42122945e-2
       p = fmaf (p, q,  0x1.4ff8acp-03f); //  1.64048523e-1
       p = fmaf (p, q, -0x1.54081ap-03f); // -1.66031078e-1
       p = fmaf (p, q, -0x1.7bf5cep-04f); // -9.27637145e-2
       p = fmaf (p, q,  0x1.1ba03ap-02f); //  2.76978403e-1*/

       p =              5.92470169e-5;
       p = fmaf (p, q,  1.61224554e-4);
       p = fmaf (p, q, -3.46481771e-4);
       p = fmaf (p, q, -1.39681227e-3);
       p = fmaf (p, q,  1.20588380e-3);
       p = fmaf (p, q,  8.69014394e-3);
       p = fmaf (p, q, -8.01387429e-3);
       p = fmaf (p, q, -5.42122945e-2);
       p = fmaf (p, q,  1.64048523e-1);
       p = fmaf (p, q, -1.66031078e-1);
       p = fmaf (p, q, -9.27637145e-2);
       p = fmaf (p, q,  2.76978403e-1);

       /* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
       d = a + 0.5;
       r = 1.0 / d;
       r = r * 0.5;
       q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
       t = q + q;
       e = (p - q) + fmaf (t, -a, 1.0); // residual: (p+1)-q*(1+2*a)
       r = fmaf (e, r, q);

       //if (a > 0x1.fffffep127f) r = 0.0f; // 3.40282347e+38 // handle INF argument

       /* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
       if (x < 0.0) {
           s = x * x;
           d = fmaf (x, x, -s);
           e = exp (s);
           r = e - r;
           r = fmaf (e, d + d, r);
           r = r + e;
       }
       return r;
    }
}

But when I try to compile this model, I get “Compilation ERROR, function(s)/method(s) not created!”. That’s a pretty opaque error. Can anybody help me debug it?

That should be fabs to avoid the ambiguous overload error messages that should have popped up.

Yes, I’ve fixed the fabs thing. But it still doesn’t compile.

Can somebody check if this function block works on some other machine? Ie, whether it’s something local like my compiler environment, or something about the function block itself? Thanks a lot!

Yes, it compiles for me if I use fabs instead of abs.

That’s a bug. I just filed a GitHub issue for stan-dev/math.