Yet another DIY RPN calculator

[Albert Chan]

erfinv code had a hidden gem
With minor changes, adding a tail = 1 - |x|, we get *accurate* inverse erfc for free.

Code:

static double erfinv(double x, double tail) {
  
  if (x + tail != 1) tail = 1-x;    // assume x >= 0
  if (tail <= 0) return tail ? NAN : INFINITY;
  ...
  double r = sqrt(M_LN2 - log(tail));
  ...

Calling C from Lua, code is trivial

Code:

static int Lierf(lua_State *L)
{
  double x = luaL_checknumber(L,1);
  double tail = lua_tonumber(L,2);
  lua_pushnumber(L, copysign(erfinv(fabs(x), tail), x));
  return 1;
}

Note that tail is suggestive, throw away if x + tail ≠ 1
This make ierfc(x) easy to implement, even if x > 1

lua> for k, v in pairs(require'mathx') do _G[k]=v end -- get erf, erfc, ierf
lua> ierfc = function(x) return ierf(1-x, x) end

lua> x = erfc(5)
lua> x, ierfc(x)
1.5374597944280351e-012     4.999999999999999

lua> x = erfc(10)
lua> x, ierfc(x)
2.088487583762545e-045      10

lua> x = erfc(15)
lua> x, ierfc(x)
7.212994172451206e-100      15

lua> x = erfc(20)
lua> x, ierfc(x)
5.3958656116079005e-176    20