--- a/NEWS
+++ b/NEWS
@@ -74,9 +74,10 @@
 
  ** Other new functions added in 3.8.0:
 
-      colorcube   lines         splinefit
-      erfcinv     rgbplot       tetramesh
-      findfigs    shrinkfaces
+      betaincinv  lines         tetramesh
+      colorcube   rgbplot
+      erfcinv     shrinkfaces
+      findfigs    splinefit
       
  ** Deprecated functions.
 

--- a/doc/interpreter/arith.txi
+++ b/doc/interpreter/arith.txi
@@ -255,6 +255,8 @@
 
 @DOCSTRING(betainc)
 
+@DOCSTRING(betaincinv)
+
 @DOCSTRING(betaln)
 
 @DOCSTRING(bincoeff)

--- a/liboctave/lo-specfun.cc
+++ b/liboctave/lo-specfun.cc
@@ -615,7 +615,7 @@
 {
   FloatComplex retval;
 
-  float r = x.real (), i = x.imag();
+  float r = x.real (), i = x.imag ();
 
   if (fabs (r) < 0.5 && fabs (i) < 0.5)
     {
@@ -2158,6 +2158,28 @@
    d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
 }
 
+static void
+gripe_betaincinv_nonconformant (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2, octave_idx_type r3,
+                                octave_idx_type c3)
+{
+  (*current_liboctave_error_handler)
+   ("betaincinv: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)",
+     r1, c1, r2, c2, r3, c3);
+}
+
+static void
+gripe_betaincinv_nonconformant (const dim_vector& d1, const dim_vector& d2,
+                                const dim_vector& d3)
+{
+  std::string d1_str = d1.str ();
+  std::string d2_str = d2.str ();
+  std::string d3_str = d3.str ();
+
+  (*current_liboctave_error_handler)
+  ("betaincinv: nonconformant arguments (x is %s, a is %s, b is %s)",
+   d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
+}
+
 double
 betainc (double x, double a, double b)
 {
@@ -2858,7 +2880,7 @@
 
       (*current_liboctave_error_handler)
         ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
-         x_str.c_str (), a_str. c_str ());
+         x_str.c_str (), a_str.c_str ());
     }
 
  done:
@@ -3124,3 +3146,476 @@
 {
   return erfcx_impl (x);
 }
+
+//
+//  Incomplete Beta function ratio
+//
+//  Algorithm based on the one by John Burkardt.
+//  See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html
+//
+//  The original code is distributed under the GNU LGPL license.
+//
+//  Reference:
+//
+//    KL Majumder, GP Bhattacharjee,
+//    Algorithm AS 63:
+//    The incomplete Beta Integral,
+//    Applied Statistics,
+//    Volume 22, Number 3, 1973, pages 409-411.
+//
+double
+betain (double x, double p, double q, double beta, bool& err)
+{
+  double acu = 0.1E-14, ai, cx;
+  bool indx;
+  int ns;
+  double pp, psq, qq, rx, temp, term, value, xx;
+
+  value = x;
+  err = false;
+
+  //  Check the input arguments.
+
+  if ((p <= 0.0 || q <= 0.0) || (x < 0.0 || 1.0 < x))
+    {
+      err = true;
+      return value;
+    }
+
+  //  Special cases.
+
+  if (x == 0.0 || x == 1.0)
+    {
+      return value;
+    }
+
+  //  Change tail if necessary and determine S.
+
+  psq = p + q;
+  cx = 1.0 - x;
+
+  if (p < psq * x)
+    {
+      xx = cx;
+      cx = x;
+      pp = q;
+      qq = p;
+      indx = true;
+    }
+  else
+    {
+      xx = x;
+      pp = p;
+      qq = q;
+      indx = false;
+    }
+
+  term = 1.0;
+  ai = 1.0;
+  value = 1.0;
+  ns = (int) (qq + cx * psq);
+
+  //  Use the Soper reduction formula.
+
+  rx = xx / cx;
+  temp = qq - ai;
+  if (ns == 0)
+    {
+      rx = xx;
+    }
+
+  for ( ; ; )
+    {
+      term = term * temp * rx / (pp + ai);
+      value = value + term;
+      temp = fabs (term);
+
+      if (temp <= acu && temp <= acu * value)
+        {
+          value = value * exp (pp * log (xx)
+          + (qq - 1.0) * log (cx) - beta) / pp;
+
+          if (indx)
+            {
+              value = 1.0 - value;
+            }
+          break;
+        }
+
+      ai = ai + 1.0;
+      ns = ns - 1;
+
+      if (0 <= ns)
+        {
+          temp = qq - ai;
+          if (ns == 0)
+            {
+              rx = xx;
+            }
+        }
+      else
+        {
+          temp = psq;
+          psq = psq + 1.0;
+        }
+    }
+
+  return value;
+}
+
+//
+//  Inverse of the incomplete Beta function
+//
+//  Algorithm based on the one by John Burkardt.
+//  See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html
+//
+//  The original code is distributed under the GNU LGPL license.
+//
+//  Reference:
+//
+//    GW Cran, KJ Martin, GE Thomas,
+//    Remark AS R19 and Algorithm AS 109:
+//    A Remark on Algorithms AS 63: The Incomplete Beta Integral
+//    and AS 64: Inverse of the Incomplete Beta Integeral,
+//    Applied Statistics,
+//    Volume 26, Number 1, 1977, pages 111-114.
+//
+double
+betaincinv (double y, double p, double q) {
+  double a, acu, adj, fpu, g, h;
+  int iex;
+  bool indx;
+  double pp, prev, qq, r, s, sae = -37.0, sq, t, tx, value, w, xin, ycur, yprev;
+
+  double beta = lgamma (p) + lgamma (q) - lgamma (p + q);
+  bool err = false;
+  fpu = pow (10.0, sae);
+  value = y;
+
+  //  Test for admissibility of parameters.
+
+  if (p <= 0.0 || q <= 0.0)
+    {
+      (*current_liboctave_error_handler)
+        ("betaincinv: wrong parameters");
+    }
+
+  if (y < 0.0 || 1.0 < y)
+    {
+      (*current_liboctave_error_handler)
+        ("betaincinv: wrong parameter Y");
+    }
+
+  if (y == 0.0 || y == 1.0)
+    {
+      return value;
+    }
+
+  //  Change tail if necessary.
+
+  if (0.5 < y)
+    {
+      a = 1.0 - y;
+      pp = q;
+      qq = p;
+      indx = true;
+    }
+  else
+    {
+      a = y;
+      pp = p;
+      qq = q;
+      indx = false;
+    }
+
+  //  Calculate the initial approximation.
+
+  r = sqrt (- log (a * a));
+
+  ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r);
+
+  if (1.0 < pp && 1.0 < qq)
+    {
+      r = (ycur * ycur - 3.0) / 6.0;
+      s = 1.0 / (pp + pp - 1.0);
+      t = 1.0 / (qq + qq - 1.0);
+      h = 2.0 / (s + t);
+      w = ycur * sqrt (h + r) / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h));
+      value = pp / (pp + qq * exp (w + w));
+    }
+  else
+    {
+      r = qq + qq;
+      t = 1.0 / (9.0 * qq);
+      t = r * pow (1.0 - t + ycur * sqrt (t), 3);
+
+      if (t <= 0.0)
+        {
+          value = 1.0 - exp ((log ((1.0 - a) * qq) + beta) / qq);
+        }
+      else
+        {
+          t = (4.0 * pp + r - 2.0) / t;
+
+          if (t <= 1.0)
+            {
+              value = exp ((log (a * pp) + beta) / pp);
+            }
+          else
+            {
+              value = 1.0 - 2.0 / (t + 1.0);
+            }
+        }
+    }
+
+  //  Solve for X by a modified Newton-Raphson method,
+  //  using the function BETAIN.
+
+  r = 1.0 - pp;
+  t = 1.0 - qq;
+  yprev = 0.0;
+  sq = 1.0;
+  prev = 1.0;
+
+  if (value < 0.0001)
+    {
+      value = 0.0001;
+    }
+
+  if (0.9999 < value)
+    {
+      value = 0.9999;
+    }
+
+  iex = std::max (- 5.0 / pp / pp - 1.0 / pow (a, 0.2) - 13.0, sae);
+
+  acu = pow (10.0, iex);
+
+  for ( ; ; )
+    {
+      ycur = betain (value, pp, qq, beta, err);
+
+      if (err)
+        {
+          return value;
+        }
+
+      xin = value;
+      ycur = (ycur - a) * exp (beta + r * log (xin) + t * log (1.0 - xin));
+
+      if (ycur * yprev <= 0.0)
+        {
+          prev = std::max (sq, fpu);
+        }
+
+      g = 1.0;
+
+      for ( ; ; )
+        {
+          for ( ; ; )
+            {
+              adj = g * ycur;
+              sq = adj * adj;
+
+              if (sq < prev)
+                {
+                  tx = value - adj;
+
+                  if (0.0 <= tx && tx <= 1.0)
+                    {
+                      break;
+                    }
+                }
+              g = g / 3.0;
+            }
+
+          if (prev <= acu)
+            {
+              if (indx)
+                {
+                  value = 1.0 - value;
+                }
+              return value;
+            }
+
+          if (ycur * ycur <= acu)
+            {
+              if (indx)
+                {
+                  value = 1.0 - value;
+                }
+              return value;
+            }
+
+          if (tx != 0.0 && tx != 1.0)
+            {
+              break;
+            }
+
+          g = g / 3.0;
+        }
+
+      if (tx == value)
+        {
+          break;
+        }
+
+      value = tx;
+      yprev = ycur;
+    }
+
+  if (indx)
+    {
+      value = 1.0 - value;
+    }
+
+  return value;
+}
+
+Matrix
+betaincinv (double x, double a, const Matrix& b)
+{
+  octave_idx_type nr = b.rows ();
+  octave_idx_type nc = b.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      retval(i,j) = betaincinv (x, a, b(i,j));
+
+  return retval;
+}
+
+Matrix
+betaincinv (double x, const Matrix& a, double b)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      retval(i,j) = betaincinv (x, a(i,j), b);
+
+  return retval;
+}
+
+Matrix
+betaincinv (double x, const Matrix& a, const Matrix& b)
+{
+  Matrix retval;
+
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  if (a_nr == b_nr && a_nc == b_nc)
+    {
+      retval.resize (a_nr, a_nc);
+
+      for (octave_idx_type j = 0; j < a_nc; j++)
+        for (octave_idx_type i = 0; i < a_nr; i++)
+          retval(i,j) = betaincinv (x, a(i,j), b(i,j));
+    }
+  else
+    gripe_betaincinv_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc);
+
+  return retval;
+}
+
+Matrix
+betaincinv (const Matrix& x, double a, double b)
+{
+  octave_idx_type nr = x.rows ();
+  octave_idx_type nc = x.cols ();
+
+  Matrix retval (nr, nc);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type i = 0; i < nr; i++)
+      retval(i,j) = betaincinv (x(i,j), a, b);
+
+  return retval;
+}
+
+Matrix
+betaincinv (const Matrix& x, double a, const Matrix& b)
+{
+  Matrix retval;
+
+  octave_idx_type nr = x.rows ();
+  octave_idx_type nc = x.cols ();
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  if (nr == b_nr && nc == b_nc)
+    {
+      retval.resize (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          retval(i,j) = betaincinv (x(i,j), a, b(i,j));
+    }
+  else
+    gripe_betaincinv_nonconformant (nr, nc, 1, 1, b_nr, b_nc);
+
+  return retval;
+}
+
+Matrix
+betaincinv (const Matrix& x, const Matrix& a, double b)
+{
+  Matrix retval;
+
+  octave_idx_type nr = x.rows ();
+  octave_idx_type nc = x.cols ();
+
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (nr == a_nr && nc == a_nc)
+    {
+      retval.resize (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          retval(i,j) = betaincinv (x(i,j), a(i,j), b);
+    }
+  else
+    gripe_betaincinv_nonconformant (nr, nc, a_nr, a_nc, 1, 1);
+
+  return retval;
+}
+
+Matrix
+betaincinv (const Matrix& x, const Matrix& a, const Matrix& b)
+{
+  Matrix retval;
+
+  octave_idx_type nr = x.rows ();
+  octave_idx_type nc = x.cols ();
+
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc)
+    {
+      retval.resize (nr, nc);
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          retval(i,j) = betaincinv (x(i,j), a(i,j), b(i,j));
+    }
+  else
+    gripe_betaincinv_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc);
+
+  return retval;
+}

--- a/liboctave/lo-specfun.h
+++ b/liboctave/lo-specfun.h
@@ -581,4 +581,16 @@
 extern OCTAVE_API double erfcx (double x);
 extern OCTAVE_API float erfcx (float x);
 
+extern OCTAVE_API double betaincinv (double y, double p, double q);
+extern OCTAVE_API Matrix betaincinv (double x, double a, const Matrix& b);
+extern OCTAVE_API Matrix betaincinv (double x, const Matrix& a, double b);
+extern OCTAVE_API Matrix betaincinv (double x, const Matrix& a, const Matrix& b);
+
+extern OCTAVE_API Matrix betaincinv (const Matrix& x, double a, double b);
+extern OCTAVE_API Matrix betaincinv (const Matrix& x, double a, const Matrix& b);
+extern OCTAVE_API Matrix betaincinv (const Matrix& x, const Matrix& a, double b);
+extern OCTAVE_API Matrix betaincinv (const Matrix& x, const Matrix& a, const Matrix& b);
+
+extern OCTAVE_API double betain (double x, double p, double q, double beta, bool& err);
+
 #endif

--- a/scripts/help/unimplemented.m
+++ b/scripts/help/unimplemented.m
@@ -98,7 +98,6 @@
   "bar3",
   "bar3h",
   "bench",
-  "betaincinv",
   "bicgstabl",
   "brush",
   "builddocsearchdb",

--- a/src/mappers.cc
+++ b/src/mappers.cc
@@ -1945,6 +1945,256 @@
 %!error tanh (1, 2)
 */
 
+DEFUN (betaincinv, args, ,
+    "-*- texinfo -*-\n\
+@deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{z}, @var{w})\n\
+Compute the inverse of the incomplete Beta function, i.e., @var{x} such that\n\
+\n\
+@example\n\
+@var{y} == betainc (@var{x}, @var{z}, @var{w}) \n\
+@end example\n\
+@seealso{beta, betainc, betaln}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin == 3)
+    {
+      octave_value x_arg = args(0);
+      octave_value a_arg = args(1);
+      octave_value b_arg = args(2);
+
+      // FIXME Can we make a template version of the duplicated code below
+      if (x_arg.is_single_type () || a_arg.is_single_type () ||
+          b_arg.is_single_type ())
+        {
+          if (x_arg.is_scalar_type ())
+            {
+              float x = x_arg.float_value ();
+
+              if (a_arg.is_scalar_type ())
+                {
+                  float a = a_arg.float_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          float b = b_arg.float_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          FloatNDArray b = b_arg.float_array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+              else
+                {
+                  FloatNDArray a = a_arg.float_array_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          float b = b_arg.float_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          FloatNDArray b = b_arg.float_array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+            }
+          else
+            {
+              FloatNDArray x = x_arg.float_array_value ();
+
+              if (a_arg.is_scalar_type ())
+                {
+                  float a = a_arg.float_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          float b = b_arg.float_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          FloatNDArray b = b_arg.float_array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+              else
+                {
+                  FloatNDArray a = a_arg.float_array_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          float b = b_arg.float_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          FloatNDArray b = b_arg.float_array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+            }
+        }
+      else
+        {
+          if (x_arg.is_scalar_type ())
+            {
+              double x = x_arg.double_value ();
+
+              if (a_arg.is_scalar_type ())
+                {
+                  double a = a_arg.double_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          double b = b_arg.double_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          NDArray b = b_arg.array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+              else
+                {
+                  NDArray a = a_arg.array_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          double b = b_arg.double_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          NDArray b = b_arg.array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+            }
+          else
+            {
+              NDArray x = x_arg.array_value ();
+
+              if (a_arg.is_scalar_type ())
+                {
+                  double a = a_arg.double_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          double b = b_arg.double_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          NDArray b = b_arg.array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+              else
+                {
+                  NDArray a = a_arg.array_value ();
+
+                  if (! error_state)
+                    {
+                      if (b_arg.is_scalar_type ())
+                        {
+                          double b = b_arg.double_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                      else
+                        {
+                          NDArray b = b_arg.array_value ();
+
+                          if (! error_state)
+                            retval = betaincinv (x, a, b);
+                        }
+                    }
+                }
+            }
+        }
+    }
+  else
+    print_usage ();
+
+  return retval;
+}
+
+/*
+%!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps));
+%!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps));
+%!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps));
+%!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps));
+%!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps));
+%!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps));
+%!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps));
+%!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps));
+%!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps));
+
+%!assert (betaincinv (0, 42, 42), 0, sqrt (eps));
+%!assert (betaincinv (1, 42, 42), 1, sqrt (eps));
+
+%!error betaincinv()
+%!error betaincinv(1, 2)
+*/
+
 DEFUNX ("toascii", Ftoascii, args, ,
     "-*- texinfo -*-\n\
 @deftypefn {Mapping Function} {} toascii (@var{s})\n\