# HG changeset patch
# User Max Brister <max@2bass.com>
# Date 1341331516 18000
# Node ID 8b4606e3be32a896333fd72d4bf377670991e36b
# Parent  f649b66ef1afddb287d5c1edf733affb594cc279# Parent  619fedc6ea61b04428393eac68348cdac471d7ef
maint: periodic merge of default to jit

diff --git a/NEWS b/NEWS
--- 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.
 
diff --git a/doc/interpreter/arith.txi b/doc/interpreter/arith.txi
--- a/doc/interpreter/arith.txi
+++ b/doc/interpreter/arith.txi
@@ -255,6 +255,8 @@
 
 @DOCSTRING(betainc)
 
+@DOCSTRING(betaincinv)
+
 @DOCSTRING(betaln)
 
 @DOCSTRING(bincoeff)
diff --git a/doc/interpreter/oop.txi b/doc/interpreter/oop.txi
--- a/doc/interpreter/oop.txi
+++ b/doc/interpreter/oop.txi
@@ -586,7 +586,7 @@
 @item @tab a >= b @tab ge (a, b) @tab Greater than or equal to operator @tab
 @item @tab a == b @tab eq (a, b) @tab Equal to operator @tab
 @item @tab a != b @tab ne (a, b) @tab Not equal to operator @tab
-@item @tab a \& b @tab and (a, b) @tab Logical and operator @tab
+@item @tab a & b @tab and (a, b) @tab Logical and operator @tab
 @item @tab a | b @tab or (a, b) @tab Logical or operator @tab
 @item @tab ! b @tab not (a) @tab Logical not operator @tab
 @item @tab a' @tab ctranspose (a) @tab Complex conjugate transpose operator @tab
diff --git a/liboctave/lo-ieee.h b/liboctave/lo-ieee.h
--- a/liboctave/lo-ieee.h
+++ b/liboctave/lo-ieee.h
@@ -118,4 +118,26 @@
 #define lo_ieee_signbit(x) (sizeof (x) == sizeof (float) ? \
                           __lo_ieee_float_signbit (x) : __lo_ieee_signbit (x))
 
+#ifdef __cplusplus
+
+template <typename T>
+struct octave_numeric_limits
+{
+  static T NA (void) { return static_cast<T> (0); }
+};
+
+template <>
+struct octave_numeric_limits<double>
+{
+  static double NA (void) { return octave_NA; }
+};
+
+template <>
+struct octave_numeric_limits<float>
+{
+  static float NA (void) { return octave_Float_NA; }
+};
+
 #endif
+
+#endif
diff --git a/liboctave/lo-specfun.cc b/liboctave/lo-specfun.cc
--- 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)
 {
@@ -2166,93 +2188,42 @@
   return retval;
 }
 
-Matrix
-betainc (double x, double a, const Matrix& b)
+Array<double>
+betainc (double x, double a, const Array<double>& 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) = betainc (x, a, b(i,j));
-
-  return retval;
-}
-
-Matrix
-betainc (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) = betainc (x, a(i,j), b);
+  dim_vector dv = b.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<double> retval (dv);
+
+  double *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betainc (x, a, b(i));
 
   return retval;
 }
 
-Matrix
-betainc (double x, const Matrix& a, const Matrix& b)
+Array<double>
+betainc (double x, const Array<double>& a, double 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) = betainc (x, a(i,j), b(i,j));
-    }
-  else
-    gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc);
+  dim_vector dv = a.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<double> retval (dv);
+
+  double *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betainc (x, a(i), b);
 
   return retval;
 }
 
-NDArray
-betainc (double x, double a, const NDArray& b)
+Array<double>
+betainc (double x, const Array<double>& a, const Array<double>& b)
 {
-  dim_vector dv = b.dims ();
-  octave_idx_type nel = dv.numel ();
-
-  NDArray retval (dv);
-
-  for (octave_idx_type i = 0; i < nel; i++)
-    retval (i) = betainc (x, a, b(i));
-
-  return retval;
-}
-
-NDArray
-betainc (double x, const NDArray& a, double b)
-{
-  dim_vector dv = a.dims ();
-  octave_idx_type nel = dv.numel ();
-
-  NDArray retval (dv);
-
-  for (octave_idx_type i = 0; i < nel; i++)
-    retval (i) = betainc (x, a(i), b);
-
-  return retval;
-}
-
-NDArray
-betainc (double x, const NDArray& a, const NDArray& b)
-{
-  NDArray retval;
+  Array<double> retval;
   dim_vector dv = a.dims ();
 
   if (dv == b.dims ())
@@ -2261,8 +2232,10 @@
 
       retval.resize (dv);
 
+      double *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x, a(i), b(i));
+        *pretval++ = betainc (x, a(i), b(i));
     }
   else
     gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ());
@@ -2270,118 +2243,26 @@
   return retval;
 }
 
-
-Matrix
-betainc (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) = betainc (x(i,j), a, b);
-
-  return retval;
-}
-
-Matrix
-betainc (const Matrix& x, double a, const Matrix& b)
+Array<double>
+betainc (const Array<double>& x, double a, double 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) = betainc (x(i,j), a, b(i,j));
-    }
-  else
-    gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc);
+  dim_vector dv = x.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<double> retval (dv);
+
+  double *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betainc (x(i), a, b);
 
   return retval;
 }
 
-Matrix
-betainc (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) = betainc (x(i,j), a(i,j), b);
-    }
-  else
-    gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1);
-
-  return retval;
-}
-
-Matrix
-betainc (const Matrix& x, const Matrix& a, const Matrix& b)
+Array<double>
+betainc (const Array<double>& x, double a, const Array<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 ();
-
-  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) = betainc (x(i,j), a(i,j), b(i,j));
-    }
-  else
-    gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc);
-
-  return retval;
-}
-
-NDArray
-betainc (const NDArray& x, double a, double b)
-{
-  dim_vector dv = x.dims ();
-  octave_idx_type nel = dv.numel ();
-
-  NDArray retval (dv);
-
-  for (octave_idx_type i = 0; i < nel; i++)
-    retval (i) = betainc (x(i), a, b);
-
-  return retval;
-}
-
-NDArray
-betainc (const NDArray& x, double a, const NDArray& b)
-{
-  NDArray retval;
+  Array<double> retval;
   dim_vector dv = x.dims ();
 
   if (dv == b.dims ())
@@ -2390,8 +2271,10 @@
 
       retval.resize (dv);
 
+      double *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x(i), a, b(i));
+        *pretval++ = betainc (x(i), a, b(i));
     }
   else
     gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ());
@@ -2399,10 +2282,10 @@
   return retval;
 }
 
-NDArray
-betainc (const NDArray& x, const NDArray& a, double b)
+Array<double>
+betainc (const Array<double>& x, const Array<double>& a, double b)
 {
-  NDArray retval;
+  Array<double> retval;
   dim_vector dv = x.dims ();
 
   if (dv == a.dims ())
@@ -2411,8 +2294,10 @@
 
       retval.resize (dv);
 
+      double *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x(i), a(i), b);
+        *pretval++ = betainc (x(i), a(i), b);
     }
   else
     gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
@@ -2420,10 +2305,10 @@
   return retval;
 }
 
-NDArray
-betainc (const NDArray& x, const NDArray& a, const NDArray& b)
+Array<double>
+betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b)
 {
-  NDArray retval;
+  Array<double> retval;
   dim_vector dv = x.dims ();
 
   if (dv == a.dims () && dv == b.dims ())
@@ -2432,8 +2317,10 @@
 
       retval.resize (dv);
 
+      double *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x(i), a(i), b(i));
+        *pretval++ = betainc (x(i), a(i), b(i));
     }
   else
     gripe_betainc_nonconformant (dv, a.dims (), b.dims ());
@@ -2449,93 +2336,42 @@
   return retval;
 }
 
-FloatMatrix
-betainc (float x, float a, const FloatMatrix& b)
+Array<float>
+betainc (float x, float a, const Array<float>& b)
 {
-  octave_idx_type nr = b.rows ();
-  octave_idx_type nc = b.cols ();
-
-  FloatMatrix retval (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      retval(i,j) = betainc (x, a, b(i,j));
-
-  return retval;
-}
-
-FloatMatrix
-betainc (float x, const FloatMatrix& a, float b)
-{
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-
-  FloatMatrix retval (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      retval(i,j) = betainc (x, a(i,j), b);
+  dim_vector dv = b.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<float> retval (dv);
+
+  float *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betainc (x, a, b(i));
 
   return retval;
 }
 
-FloatMatrix
-betainc (float x, const FloatMatrix& a, const FloatMatrix& b)
+Array<float>
+betainc (float x, const Array<float>& a, float b)
 {
-  FloatMatrix 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) = betainc (x, a(i,j), b(i,j));
-    }
-  else
-    gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc);
+  dim_vector dv = a.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<float> retval (dv);
+
+  float *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betainc (x, a(i), b);
 
   return retval;
 }
 
-FloatNDArray
-betainc (float x, float a, const FloatNDArray& b)
+Array<float>
+betainc (float x, const Array<float>& a, const Array<float>& b)
 {
-  dim_vector dv = b.dims ();
-  octave_idx_type nel = dv.numel ();
-
-  FloatNDArray retval (dv);
-
-  for (octave_idx_type i = 0; i < nel; i++)
-    retval (i) = betainc (x, a, b(i));
-
-  return retval;
-}
-
-FloatNDArray
-betainc (float x, const FloatNDArray& a, float b)
-{
-  dim_vector dv = a.dims ();
-  octave_idx_type nel = dv.numel ();
-
-  FloatNDArray retval (dv);
-
-  for (octave_idx_type i = 0; i < nel; i++)
-    retval (i) = betainc (x, a(i), b);
-
-  return retval;
-}
-
-FloatNDArray
-betainc (float x, const FloatNDArray& a, const FloatNDArray& b)
-{
-  FloatNDArray retval;
+  Array<float> retval;
   dim_vector dv = a.dims ();
 
   if (dv == b.dims ())
@@ -2544,8 +2380,10 @@
 
       retval.resize (dv);
 
+      float *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x, a(i), b(i));
+        *pretval++ = betainc (x, a(i), b(i));
     }
   else
     gripe_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ());
@@ -2553,118 +2391,26 @@
   return retval;
 }
 
-
-FloatMatrix
-betainc (const FloatMatrix& x, float a, float b)
-{
-  octave_idx_type nr = x.rows ();
-  octave_idx_type nc = x.cols ();
-
-  FloatMatrix retval (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      retval(i,j) = betainc (x(i,j), a, b);
-
-  return retval;
-}
-
-FloatMatrix
-betainc (const FloatMatrix& x, float a, const FloatMatrix& b)
+Array<float>
+betainc (const Array<float>& x, float a, float b)
 {
-  FloatMatrix 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) = betainc (x(i,j), a, b(i,j));
-    }
-  else
-    gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc);
+  dim_vector dv = x.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<float> retval (dv);
+
+  float *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betainc (x(i), a, b);
 
   return retval;
 }
 
-FloatMatrix
-betainc (const FloatMatrix& x, const FloatMatrix& a, float b)
-{
-  FloatMatrix 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) = betainc (x(i,j), a(i,j), b);
-    }
-  else
-    gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1);
-
-  return retval;
-}
-
-FloatMatrix
-betainc (const FloatMatrix& x, const FloatMatrix& a, const FloatMatrix& b)
+Array<float>
+betainc (const Array<float>& x, float a, const Array<float>& b)
 {
-  FloatMatrix 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) = betainc (x(i,j), a(i,j), b(i,j));
-    }
-  else
-    gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc);
-
-  return retval;
-}
-
-FloatNDArray
-betainc (const FloatNDArray& x, float a, float b)
-{
-  dim_vector dv = x.dims ();
-  octave_idx_type nel = dv.numel ();
-
-  FloatNDArray retval (dv);
-
-  for (octave_idx_type i = 0; i < nel; i++)
-    retval (i) = betainc (x(i), a, b);
-
-  return retval;
-}
-
-FloatNDArray
-betainc (const FloatNDArray& x, float a, const FloatNDArray& b)
-{
-  FloatNDArray retval;
+  Array<float> retval;
   dim_vector dv = x.dims ();
 
   if (dv == b.dims ())
@@ -2673,8 +2419,10 @@
 
       retval.resize (dv);
 
+      float *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x(i), a, b(i));
+        *pretval++ = betainc (x(i), a, b(i));
     }
   else
     gripe_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ());
@@ -2682,10 +2430,10 @@
   return retval;
 }
 
-FloatNDArray
-betainc (const FloatNDArray& x, const FloatNDArray& a, float b)
+Array<float>
+betainc (const Array<float>& x, const Array<float>& a, float b)
 {
-  FloatNDArray retval;
+  Array<float> retval;
   dim_vector dv = x.dims ();
 
   if (dv == a.dims ())
@@ -2694,8 +2442,10 @@
 
       retval.resize (dv);
 
+      float *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x(i), a(i), b);
+        *pretval++ = betainc (x(i), a(i), b);
     }
   else
     gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
@@ -2703,10 +2453,10 @@
   return retval;
 }
 
-FloatNDArray
-betainc (const FloatNDArray& x, const FloatNDArray& a, const FloatNDArray& b)
+Array<float>
+betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b)
 {
-  FloatNDArray retval;
+  Array<float> retval;
   dim_vector dv = x.dims ();
 
   if (dv == a.dims () && dv == b.dims ())
@@ -2715,8 +2465,10 @@
 
       retval.resize (dv);
 
+      float *pretval = retval.fortran_vec ();
+
       for (octave_idx_type i = 0; i < nel; i++)
-        retval (i) = betainc (x(i), a(i), b(i));
+        *pretval++ = betainc (x(i), a(i), b(i));
     }
   else
     gripe_betainc_nonconformant (dv, a.dims (), b.dims ());
@@ -3128,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:
@@ -3394,3 +3146,468 @@
 {
   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 v3 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 v3 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;
+}
+
+Array<double>
+betaincinv (double x, double a, const Array<double>& b)
+{
+  dim_vector dv = b.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<double> retval (dv);
+
+  double *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betaincinv (x, a, b(i));
+
+  return retval;
+}
+
+Array<double>
+betaincinv (double x, const Array<double>& a, double b)
+{
+  dim_vector dv = a.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<double> retval (dv);
+
+  double *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betaincinv (x, a(i), b);
+
+  return retval;
+}
+
+Array<double>
+betaincinv (double x, const Array<double>& a, const Array<double>& b)
+{
+  Array<double> retval;
+  dim_vector dv = a.dims ();
+
+  if (dv == b.dims ())
+    {
+      octave_idx_type nel = dv.numel ();
+
+      retval.resize (dv);
+
+      double *pretval = retval.fortran_vec ();
+
+      for (octave_idx_type i = 0; i < nel; i++)
+        *pretval++ = betaincinv (x, a(i), b(i));
+    }
+  else
+    gripe_betaincinv_nonconformant (dim_vector (0, 0), dv, b.dims ());
+
+  return retval;
+}
+
+Array<double>
+betaincinv (const Array<double>& x, double a, double b)
+{
+  dim_vector dv = x.dims ();
+  octave_idx_type nel = dv.numel ();
+
+  Array<double> retval (dv);
+
+  double *pretval = retval.fortran_vec ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    *pretval++ = betaincinv (x(i), a, b);
+
+  return retval;
+}
+
+Array<double>
+betaincinv (const Array<double>& x, double a, const Array<double>& b)
+{
+  Array<double> retval;
+  dim_vector dv = x.dims ();
+
+  if (dv == b.dims ())
+    {
+      octave_idx_type nel = dv.numel ();
+
+      retval.resize (dv);
+
+      double *pretval = retval.fortran_vec ();
+
+      for (octave_idx_type i = 0; i < nel; i++)
+        *pretval++ = betaincinv (x(i), a, b(i));
+    }
+  else
+    gripe_betaincinv_nonconformant (dv, dim_vector (0, 0), b.dims ());
+
+  return retval;
+}
+
+Array<double>
+betaincinv (const Array<double>& x, const Array<double>& a, double b)
+{
+  Array<double> retval;
+  dim_vector dv = x.dims ();
+
+  if (dv == a.dims ())
+    {
+      octave_idx_type nel = dv.numel ();
+
+      retval.resize (dv);
+
+      double *pretval = retval.fortran_vec ();
+
+      for (octave_idx_type i = 0; i < nel; i++)
+        *pretval++ = betaincinv (x(i), a(i), b);
+    }
+  else
+    gripe_betaincinv_nonconformant (dv, a.dims (), dim_vector (0, 0));
+
+  return retval;
+}
+
+Array<double>
+betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b)
+{
+  Array<double> retval;
+  dim_vector dv = x.dims ();
+
+  if (dv == a.dims () && dv == b.dims ())
+    {
+      octave_idx_type nel = dv.numel ();
+
+      retval.resize (dv);
+
+      double *pretval = retval.fortran_vec ();
+
+      for (octave_idx_type i = 0; i < nel; i++)
+        *pretval++ = betaincinv (x(i), a(i), b(i));
+    }
+  else
+    gripe_betaincinv_nonconformant (dv, a.dims (), b.dims ());
+
+  return retval;
+}
diff --git a/liboctave/lo-specfun.h b/liboctave/lo-specfun.h
--- a/liboctave/lo-specfun.h
+++ b/liboctave/lo-specfun.h
@@ -520,42 +520,24 @@
 biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr);
 
 extern OCTAVE_API double betainc (double x, double a, double b);
-extern OCTAVE_API Matrix betainc (double x, double a, const Matrix& b);
-extern OCTAVE_API Matrix betainc (double x, const Matrix& a, double b);
-extern OCTAVE_API Matrix betainc (double x, const Matrix& a, const Matrix& b);
-
-extern OCTAVE_API NDArray betainc (double x, double a, const NDArray& b);
-extern OCTAVE_API NDArray betainc (double x, const NDArray& a, double b);
-extern OCTAVE_API NDArray betainc (double x, const NDArray& a, const NDArray& b);
-
-extern OCTAVE_API Matrix betainc (const Matrix& x, double a, double b);
-extern OCTAVE_API Matrix betainc (const Matrix& x, double a, const Matrix& b);
-extern OCTAVE_API Matrix betainc (const Matrix& x, const Matrix& a, double b);
-extern OCTAVE_API Matrix betainc (const Matrix& x, const Matrix& a, const Matrix& b);
-
-extern OCTAVE_API NDArray betainc (const NDArray& x, double a, double b);
-extern OCTAVE_API NDArray betainc (const NDArray& x, double a, const NDArray& b);
-extern OCTAVE_API NDArray betainc (const NDArray& x, const NDArray& a, double b);
-extern OCTAVE_API NDArray betainc (const NDArray& x, const NDArray& a, const NDArray& b);
+extern OCTAVE_API Array<double> betainc (double x, double a, const Array<double>& b);
+extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, double b);
+extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, const Array<double>& b);
+extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, double b);
+extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, double b);
+extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, const Array<double>& b);
+extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, double b);
+extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b);
 
 extern OCTAVE_API float betainc (float x, float a, float b);
-extern OCTAVE_API FloatMatrix betainc (float x, float a, const FloatMatrix& b);
-extern OCTAVE_API FloatMatrix betainc (float x, const FloatMatrix& a, float b);
-extern OCTAVE_API FloatMatrix betainc (float x, const FloatMatrix& a, const FloatMatrix& b);
-
-extern OCTAVE_API FloatNDArray betainc (float x, float a, const FloatNDArray& b);
-extern OCTAVE_API FloatNDArray betainc (float x, const FloatNDArray& a, float b);
-extern OCTAVE_API FloatNDArray betainc (float x, const FloatNDArray& a, const FloatNDArray& b);
-
-extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, float a, float b);
-extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, float a, const FloatMatrix& b);
-extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, const FloatMatrix& a, float b);
-extern OCTAVE_API FloatMatrix betainc (const FloatMatrix& x, const FloatMatrix& a, const FloatMatrix& b);
-
-extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, float a, float b);
-extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, float a, const FloatNDArray& b);
-extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, const FloatNDArray& a, float b);
-extern OCTAVE_API FloatNDArray betainc (const FloatNDArray& x, const FloatNDArray& a, const FloatNDArray& b);
+extern OCTAVE_API Array<float> betainc (float x, float a, const Array<float>& b);
+extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, float b);
+extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, const Array<float>& b);
+extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, float b);
+extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, float b);
+extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, const Array<float>& b);
+extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, float b);
+extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b);
 
 extern OCTAVE_API double gammainc (double x, double a, bool& err);
 extern OCTAVE_API Matrix gammainc (double x, const Matrix& a);
@@ -599,4 +581,14 @@
 extern OCTAVE_API double erfcx (double x);
 extern OCTAVE_API float erfcx (float x);
 
+extern OCTAVE_API double betaincinv (double x, double a, double b);
+extern OCTAVE_API Array<double> betaincinv (double x, double a, const Array<double>& b);
+extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, double b);
+extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b);
+extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b);
+extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b);
+extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, const Array<double>& b);
+extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, double b);
+extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b);
+
 #endif
diff --git a/liboctave/lo-utils.cc b/liboctave/lo-utils.cc
--- a/liboctave/lo-utils.cc
+++ b/liboctave/lo-utils.cc
@@ -195,10 +195,14 @@
   return retval;
 }
 
-static inline double
-read_inf_nan_na (std::istream& is, char c0, char sign = '+')
+// Note that the caller is responsible for repositioning the stream on
+// failure.
+
+template <typename T>
+T
+read_inf_nan_na (std::istream& is, char c0)
 {
-  double d = 0.0;
+  T val = 0.0;
 
   switch (c0)
     {
@@ -209,21 +213,12 @@
           {
             char c2 = is.get ();
             if (c2 == 'f' || c2 == 'F')
-              d = sign == '-' ? -octave_Inf : octave_Inf;
+              val = std::numeric_limits<T>::infinity ();
             else
-              {
-                is.putback (c2);
-                is.putback (c1);
-                is.putback (c0);
-                is.setstate (std::ios::failbit);
-              }
+              is.setstate (std::ios::failbit);
           }
         else
-          {
-            is.putback (c1);
-            is.putback (c0);
-            is.setstate (std::ios::failbit);
-          }
+          is.setstate (std::ios::failbit);
       }
       break;
 
@@ -234,19 +229,12 @@
           {
             char c2 = is.get ();
             if (c2 == 'n' || c2 == 'N')
-              d = octave_NaN;
+              val = std::numeric_limits<T>::quiet_NaN ();
             else
-              {
-                is.putback (c2);
-                d = octave_NA;
-              }
+              val = octave_numeric_limits<T>::NA ();
           }
         else
-          {
-            is.putback (c1);
-            is.putback (c0);
-            is.setstate (std::ios::failbit);
-          }
+          is.setstate (std::ios::failbit);
       }
       break;
 
@@ -254,74 +242,80 @@
       abort ();
     }
 
-  return d;
+  return val;
 }
 
 // Read a double value.  Discard any sign on NaN and NA.
 
-template <>
+template <typename T>
 double
-octave_read_value (std::istream& is)
+octave_read_fp_value (std::istream& is)
 {
-  double d = 0.0;
+  T val = 0.0;
+
+  // FIXME -- resetting stream position is likely to fail unless we are
+  // reading from a file.
+  std::ios::streampos pos = is.tellg ();
 
   char c1 = ' ';
 
   while (isspace (c1))
     c1 = is.get ();
 
+  bool neg = false;
+
   switch (c1)
     {
     case '-':
-      {
-        char c2 = 0;
-        c2 = is.get ();
-        if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N')
-          d = read_inf_nan_na (is, c2, c1);
-        else
-          {
-            is.putback (c2);
-            is.putback (c1);
-            is >> d;
-          }
-      }
-      break;
+      neg = true;
+      // fall through...
 
     case '+':
       {
         char c2 = 0;
         c2 = is.get ();
         if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N')
-          d = read_inf_nan_na (is, c2, c1);
+          val = read_inf_nan_na<T> (is, c2);
         else
           {
             is.putback (c2);
-            is.putback (c1);
-            is >> d;
+            is >> val;
           }
+
+        if (neg && ! is.fail ())
+          val = -val;
       }
       break;
 
     case 'i': case 'I':
     case 'n': case 'N':
-      d = read_inf_nan_na (is, c1);
+      val = read_inf_nan_na<T> (is, c1);
       break;
 
     default:
       is.putback (c1);
-      is >> d;
+      is >> val;
+      break;
     }
 
-  return d;
+  std::ios::iostate status = is.rdstate ();
+  if (status & std::ios::failbit)
+    {
+      is.clear ();
+      is.seekg (pos);
+      is.setstate (status);
+    }
+
+  return val;
 }
 
-template <>
-Complex
-octave_read_value (std::istream& is)
+template <typename T>
+std::complex<T>
+octave_read_cx_fp_value (std::istream& is)
 {
-  double re = 0.0, im = 0.0;
+  T re = 0.0, im = 0.0;
 
-  Complex cx = 0.0;
+  std::complex<T> cx = 0.0;
 
   char ch = ' ';
 
@@ -330,16 +324,16 @@
 
   if (ch == '(')
     {
-      re = octave_read_value<double> (is);
+      re = octave_read_value<T> (is);
       ch = is.get ();
 
       if (ch == ',')
         {
-          im = octave_read_value<double> (is);
+          im = octave_read_value<T> (is);
           ch = is.get ();
 
           if (ch == ')')
-            cx = Complex (re, im);
+            cx = std::complex<T> (re, im);
           else
             is.setstate (std::ios::failbit);
         }
@@ -355,170 +349,26 @@
     }
 
   return cx;
-
 }
 
-static inline float
-read_float_inf_nan_na (std::istream& is, char c0, char sign = '+')
+template <> OCTAVE_API double octave_read_value (std::istream& is)
 {
-  float d = 0.0;
-
-  switch (c0)
-    {
-    case 'i': case 'I':
-      {
-        char c1 = is.get ();
-        if (c1 == 'n' || c1 == 'N')
-          {
-            char c2 = is.get ();
-            if (c2 == 'f' || c2 == 'F')
-              d = sign == '-' ? -octave_Float_Inf : octave_Float_Inf;
-            else
-              {
-                is.putback (c2);
-                is.putback (c1);
-                is.putback (c0);
-                is.setstate (std::ios::failbit);
-              }
-          }
-        else
-          {
-            is.putback (c1);
-            is.putback (c0);
-            is.setstate (std::ios::failbit);
-          }
-      }
-      break;
-
-    case 'n': case 'N':
-      {
-        char c1 = is.get ();
-        if (c1 == 'a' || c1 == 'A')
-          {
-            char c2 = is.get ();
-            if (c2 == 'n' || c2 == 'N')
-              d = octave_Float_NaN;
-            else
-              {
-                is.putback (c2);
-                d = octave_Float_NA;
-              }
-          }
-        else
-          {
-            is.putback (c1);
-            is.putback (c0);
-            is.setstate (std::ios::failbit);
-          }
-      }
-      break;
-
-    default:
-      abort ();
-    }
-
-  return d;
+  return octave_read_fp_value<double> (is);
 }
 
-// Read a float value.  Discard any sign on NaN and NA.
-
-template <>
-float
-octave_read_value (std::istream& is)
+template <> OCTAVE_API Complex octave_read_value (std::istream& is)
 {
-  float d = 0.0;
-
-  char c1 = ' ';
-
-  while (isspace (c1))
-    c1 = is.get ();
-
-  switch (c1)
-    {
-    case '-':
-      {
-        char c2 = 0;
-        c2 = is.get ();
-        if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N')
-          d = read_float_inf_nan_na (is, c2, c1);
-        else
-          {
-            is.putback (c2);
-            is.putback (c1);
-            is >> d;
-          }
-      }
-      break;
-
-    case '+':
-      {
-        char c2 = 0;
-        c2 = is.get ();
-        if (c2 == 'i' || c2 == 'I' || c2 == 'n' || c2 == 'N')
-          d = read_float_inf_nan_na (is, c2, c1);
-        else
-          {
-            is.putback (c2);
-            is.putback (c1);
-            is >> d;
-          }
-      }
-      break;
-
-    case 'i': case 'I':
-    case 'n': case 'N':
-      d = read_float_inf_nan_na (is, c1);
-      break;
-
-    default:
-      is.putback (c1);
-      is >> d;
-    }
-
-  return d;
+  return octave_read_cx_fp_value<double> (is);
 }
 
-template <>
-FloatComplex
-octave_read_value (std::istream& is)
+template <> OCTAVE_API float octave_read_value (std::istream& is)
 {
-  float re = 0.0, im = 0.0;
-
-  FloatComplex cx = 0.0;
-
-  char ch = ' ';
-
-  while (isspace (ch))
-    ch = is.get ();
-
-  if (ch == '(')
-    {
-      re = octave_read_value<float> (is);
-      ch = is.get ();
+  return octave_read_fp_value<float> (is);
+}
 
-      if (ch == ',')
-        {
-          im = octave_read_value<float> (is);
-          ch = is.get ();
-
-          if (ch == ')')
-            cx = FloatComplex (re, im);
-          else
-            is.setstate (std::ios::failbit);
-        }
-      else if (ch == ')')
-        cx = re;
-      else
-        is.setstate (std::ios::failbit);
-    }
-  else
-    {
-      is.putback (ch);
-      cx = octave_read_value<float> (is);
-    }
-
-  return cx;
-
+template <> OCTAVE_API FloatComplex octave_read_value (std::istream& is)
+{
+  return octave_read_cx_fp_value<float> (is);
 }
 
 void
diff --git a/scripts/help/unimplemented.m b/scripts/help/unimplemented.m
--- a/scripts/help/unimplemented.m
+++ b/scripts/help/unimplemented.m
@@ -98,7 +98,6 @@
   "bar3",
   "bar3h",
   "bench",
-  "betaincinv",
   "bicgstabl",
   "brush",
   "builddocsearchdb",
diff --git a/scripts/io/strread.m b/scripts/io/strread.m
--- a/scripts/io/strread.m
+++ b/scripts/io/strread.m
@@ -301,6 +301,13 @@
     ## Format conversion specifiers following literals w/o space/delim
     ## in between are separate now.  Separate those w trailing literals
     idy2 = find (! cellfun ("isempty", strfind (fmt_words, "%")));
+
+    ## Check for unsupported format specifiers
+    errpat = '(\[.*\]|[cq]|[nfdu]8|[nfdu]16|[nfdu]32|[nfdu]64)';
+    if (! all (cellfun ("isempty", regexp (fmt_words(idy2), errpat))))
+      error ("strread: %q, %c, %[] or bit width format specifiers are not supported yet.");
+    endif
+
     a = strfind (fmt_words(idy2), "%");
     b = regexp (fmt_words(idy2), '[nfdus]', 'end');
     for jj = 1:numel (a)
@@ -931,3 +938,19 @@
 %! assert (isempty (b));
 %! assert (isempty (c));
 
+%% Unsupported format specifiers
+%!test
+%!error <format specifiers are not supported> strread ('a', '%c')
+%!error <format specifiers are not supported> strread ('a', '%*c %d')
+%!error <format specifiers are not supported> strread ('a', '%q')
+%!error <format specifiers are not supported> strread ('a', '%*q %d')
+%!error <format specifiers are not supported> strread ('a', '%[a]')
+%!error <format specifiers are not supported> strread ('a', '%*[a] %d')
+%!error <format specifiers are not supported> strread ('a', '%[^a]')
+%!error <format specifiers are not supported> strread ('a', '%*[â] %d')
+%!error <format specifiers are not supported> strread ('a', '%d8')
+%!error <format specifiers are not supported> strread ('a', '%*d8 %s')
+%!error <format specifiers are not supported> strread ('a', '%f64')
+%!error <format specifiers are not supported> strread ('a', '%*f64 %s')
+%!error <format specifiers are not supported> strread ('a', '%u32')
+%!error <format specifiers are not supported> strread ('a', '%*u32 %d')
diff --git a/scripts/miscellaneous/perl.m b/scripts/miscellaneous/perl.m
--- a/scripts/miscellaneous/perl.m
+++ b/scripts/miscellaneous/perl.m
@@ -17,13 +17,14 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {[@var{output}, @var{status}] =} perl (@var{scriptfile})
-## @deftypefnx {Function File} {[@var{output}, @var{status}] =} perl (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{})
-## Invoke Perl script @var{scriptfile} with possibly a list of
-## command line arguments.
-## Returns output in @var{output} and status
-## in @var{status}.
-## @seealso{system}
+## @deftypefn  {Function File} {@var{output} =} perl (@var{scriptfile})
+## @deftypefnx {Function File} {@var{output} =} perl (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{})
+## @deftypefnx {Function File} {[@var{output}, @var{status}] =} perl (@dots{})
+## Invoke Perl script @var{scriptfile}, possibly with a list of command line
+## arguments.  Return output in @var{output} and optional status in
+## @var{status}.  If @var{scriptfile} is not an absolute file name it is
+## is searched for in the current directory and then in the Octave loadpath.
+## @seealso{system, python}
 ## @end deftypefn
 
 function [output, status] = perl (scriptfile = "-e ''", varargin)
@@ -33,8 +34,13 @@
   ## initial value in the argument list of the function definition.
 
   if (ischar (scriptfile)
-      && ((nargin != 1 && iscellstr (varargin))
-          || (nargin == 1 && ! isempty (scriptfile))))
+      && (   (nargin == 1 && ! isempty (scriptfile))
+          || (nargin != 1 && iscellstr (varargin))))
+    if (! strcmp (scriptfile(1:2), "-e"))
+      ## Attempt to find file in loadpath.  No effect for absolute filenames.
+      scriptfile = file_in_loadpath (scriptfile);
+    endif
+
     [status, output] = system (cstrcat ("perl ", scriptfile,
                                         sprintf (" %s", varargin{:})));
   else
diff --git a/scripts/miscellaneous/python.m b/scripts/miscellaneous/python.m
--- a/scripts/miscellaneous/python.m
+++ b/scripts/miscellaneous/python.m
@@ -18,26 +18,28 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {[@var{output}, @var{status}] =} python (@var{scriptfile})
-## @deftypefnx {Function File} {[@var{output}, @var{status}] =} python (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{})
-## Invoke python script @var{scriptfile} with possibly a list of
-## command line arguments.
-## Returns output in @var{output} and status
-## in @var{status}.
-## @seealso{system}
+## @deftypefn  {Function File} {@var{output} =} python (@var{scriptfile})
+## @deftypefnx {Function File} {@var{output} =} python (@var{scriptfile}, @var{argument1}, @var{argument2}, @dots{})
+## @deftypefnx {Function File} {[@var{output}, @var{status}] =} python (@dots{})
+## Invoke Python script @var{scriptfile}, possibly with a list of command line
+## arguments.  Return output in @var{output} and optional status in
+## @var{status}.  If @var{scriptfile} is not an absolute file name it is
+## is searched for in the current directory and then in the Octave loadpath.
+## @seealso{system, perl}
 ## @end deftypefn
 
 ## Author: Carnë Draug <carandraug+dev@gmail.com>
 
 function [output, status] = python (scriptfile = "-c ''", varargin)
 
-  ## VARARGIN is intialized to {}(1x0) if no additional arguments are
-  ## supplied, so there is no need to check for it, or provide an
-  ## initial value in the argument list of the function definition.
+  if (ischar (scriptfile)
+      && (   (nargin == 1 && ! isempty (scriptfile))
+          || (nargin != 1 && iscellstr (varargin))))
+    if (! strcmp (scriptfile(1:2), "-c"))
+      ## Attempt to find file in loadpath.  No effect for absolute filenames.
+      scriptfile = file_in_loadpath (scriptfile);
+    endif
 
-  if (ischar (scriptfile)
-      && ((nargin != 1 && iscellstr (varargin))
-          || (nargin == 1 && ! isempty (scriptfile))))
     [status, output] = system (cstrcat ("python ", scriptfile,
                                         sprintf (" %s", varargin{:})));
   else
diff --git a/scripts/polynomial/ppval.m b/scripts/polynomial/ppval.m
--- a/scripts/polynomial/ppval.m
+++ b/scripts/polynomial/ppval.m
@@ -68,7 +68,7 @@
   ndv = length (dimvec);
 
   ## Offsets.
-  dx = (xi - x(idx));
+  dx = (xi - x(idx))(:)';
   dx = repmat (dx, [prod(d), 1]);
   dx = reshape (dx, dimvec);
   dx = shiftdim (dx, ndv - 1);
@@ -119,4 +119,11 @@
 %!assert (ppval (pp2, xi), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr)
 %!assert (ppval (pp2, xi'), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr)
 %!assert (size (ppval (pp2, [xi;xi])), [2 2 4])
-
+%!test
+%! breaks = [0, 1, 2, 3];
+%! coefs = rand (6, 4);
+%! pp = mkpp (breaks, coefs, 2);
+%! ret = zeros (2, 4, 2);
+%! ret(:,:,1) = ppval (pp, breaks');
+%! ret(:,:,2) = ppval (pp, breaks');
+%! assert (ppval (pp, [breaks',breaks']), ret)
diff --git a/scripts/sparse/gmres.m b/scripts/sparse/gmres.m
--- a/scripts/sparse/gmres.m
+++ b/scripts/sparse/gmres.m
@@ -174,7 +174,7 @@
 
     x = x_old + V(:, 1:restart_it) * Y(1:restart_it);
 
-    resvec(iter) = presn;
+    resvec(iter+1) = presn;
     if (norm (x - x_old, inf) <= eps)
       flag = 3;  # Stagnation: no change between iterations
       break;
@@ -191,8 +191,7 @@
       flag = 0;  # Converged to solution within tolerance
     endif
 
-    resvec = resvec(1:iter-1);
-    it = [ceil(iter / restart), rem(iter, restart)];
+    it = [floor(iter/restart), restart_it-1];
   endif
 
 endfunction
diff --git a/scripts/specfun/beta.m b/scripts/specfun/beta.m
--- a/scripts/specfun/beta.m
+++ b/scripts/specfun/beta.m
@@ -31,6 +31,7 @@
 ## @end example
 ##
 ## @end ifnottex
+## @seealso{betaln, betainc}
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
diff --git a/scripts/strings/strcat.m b/scripts/strings/strcat.m
--- a/scripts/strings/strcat.m
+++ b/scripts/strings/strcat.m
@@ -23,11 +23,35 @@
 ## horizontally.  If the arguments are cells strings,  @code{strcat}
 ## returns a cell string with the individual cells concatenated.
 ## For numerical input, each element is converted to the
-## corresponding ASCII character.  Trailing white space is eliminated.
+## corresponding ASCII character.  Trailing white space for each of
+## the inputs (@var{s1}, @var{S2}, @dots{}) is eliminated before they
+## are concatenated.
+##
 ## For example:
 ##
 ## @example
 ## @group
+## strcat ("|", " leading space is preserved", "|")
+##     @result{} | leading space is perserved|
+## @end group
+## @end example
+##
+## @example
+## @group
+## strcat ("|", "trailing space is eliminated ", "|")
+##     @result{} |trailing space is eliminated|
+## @end group
+## @end example
+##
+## @example
+## @group
+## strcat ("homogeneous space |", "  ", "| is also eliminated")
+##     @result{} homogeneous space || is also eliminated
+## @end group
+## @end example
+##
+## @example
+## @group
 ## s = [ "ab"; "cde" ];
 ## strcat (s, s, s)
 ##     @result{}
diff --git a/src/DLD-FUNCTIONS/betainc.cc b/src/DLD-FUNCTIONS/betainc.cc
--- a/src/DLD-FUNCTIONS/betainc.cc
+++ b/src/DLD-FUNCTIONS/betainc.cc
@@ -32,6 +32,9 @@
 #include "oct-obj.h"
 #include "utils.h"
 
+// FIXME: These functions do not need to be dynamically loaded.  They should
+//        be placed elsewhere in the Octave code hierarchy.
+
 DEFUN_DLD (betainc, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Mapping Function} {} betainc (@var{x}, @var{a}, @var{b})\n\
@@ -94,7 +97,7 @@
                         }
                       else
                         {
-                          FloatNDArray b = b_arg.float_array_value ();
+                          Array<float> b = b_arg.float_array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -103,7 +106,7 @@
                 }
               else
                 {
-                  FloatNDArray a = a_arg.float_array_value ();
+                  Array<float> a = a_arg.float_array_value ();
 
                   if (! error_state)
                     {
@@ -116,7 +119,7 @@
                         }
                       else
                         {
-                          FloatNDArray b = b_arg.float_array_value ();
+                          Array<float> b = b_arg.float_array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -126,7 +129,7 @@
             }
           else
             {
-              FloatNDArray x = x_arg.float_array_value ();
+              Array<float> x = x_arg.float_array_value ();
 
               if (a_arg.is_scalar_type ())
                 {
@@ -143,7 +146,7 @@
                         }
                       else
                         {
-                          FloatNDArray b = b_arg.float_array_value ();
+                          Array<float> b = b_arg.float_array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -152,7 +155,7 @@
                 }
               else
                 {
-                  FloatNDArray a = a_arg.float_array_value ();
+                  Array<float> a = a_arg.float_array_value ();
 
                   if (! error_state)
                     {
@@ -165,7 +168,7 @@
                         }
                       else
                         {
-                          FloatNDArray b = b_arg.float_array_value ();
+                          Array<float> b = b_arg.float_array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -195,7 +198,7 @@
                         }
                       else
                         {
-                          NDArray b = b_arg.array_value ();
+                          Array<double> b = b_arg.array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -204,7 +207,7 @@
                 }
               else
                 {
-                  NDArray a = a_arg.array_value ();
+                  Array<double> a = a_arg.array_value ();
 
                   if (! error_state)
                     {
@@ -217,7 +220,7 @@
                         }
                       else
                         {
-                          NDArray b = b_arg.array_value ();
+                          Array<double> b = b_arg.array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -227,7 +230,7 @@
             }
           else
             {
-              NDArray x = x_arg.array_value ();
+              Array<double> x = x_arg.array_value ();
 
               if (a_arg.is_scalar_type ())
                 {
@@ -244,7 +247,7 @@
                         }
                       else
                         {
-                          NDArray b = b_arg.array_value ();
+                          Array<double> b = b_arg.array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -253,7 +256,7 @@
                 }
               else
                 {
-                  NDArray a = a_arg.array_value ();
+                  Array<double> a = a_arg.array_value ();
 
                   if (! error_state)
                     {
@@ -266,7 +269,7 @@
                         }
                       else
                         {
-                          NDArray b = b_arg.array_value ();
+                          Array<double> b = b_arg.array_value ();
 
                           if (! error_state)
                             retval = betainc (x, a, b);
@@ -283,7 +286,7 @@
 }
 
 /*
-%% test/octave.test/arith/betainc-1.m
+## Double precision
 %!test
 %! a = [1, 1.5, 2, 3];
 %! b = [4, 3, 2, 1];
@@ -295,7 +298,7 @@
 %! assert (v1, v2, sqrt (eps));
 %! assert (v3, v4, sqrt (eps));
 
-%% Single precision
+## Single precision
 %!test
 %! a = single ([1, 1.5, 2, 3]);
 %! b = single ([4, 3, 2, 1]);
@@ -307,7 +310,7 @@
 %! assert (v1, v2, sqrt (eps ("single")));
 %! assert (v3, v4, sqrt (eps ("single")));
 
-%% Mixed double/single precision
+## Mixed double/single precision
 %!test
 %! a = single ([1, 1.5, 2, 3]);
 %! b = [4, 3, 2, 1];
@@ -316,15 +319,172 @@
 %! x = [.2, .4, .6, .8];
 %! v3 = betainc (x, a, b);
 %! v4 = 1-betainc (1.-x, b, a);
-%! assert (v1, v2, sqrt (eps ('single')));
-%! assert (v3, v4, sqrt (eps ('single')));
+%! assert (v1, v2, sqrt (eps ("single")));
+%! assert (v3, v4, sqrt (eps ("single")));
+
+%!error betainc ()
+%!error betainc (1)
+%!error betainc (1,2)
+%!error betainc (1,2,3,4)
+*/
+
+DEFUN_DLD (betaincinv, args, ,
+    "-*- texinfo -*-\n\
+@deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{a}, @var{b})\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{a}, @var{b}) \n\
+@end example\n\
+@seealso{betainc, beta, 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);
+
+      if (x_arg.is_scalar_type ())
+        {
+          double x = x_arg.double_value ();
+
+          if (a_arg.is_scalar_type ())
+            {
+              double a = a_arg.double_value ();
 
-%% test/octave.test/arith/betainc-2.m
-%!error betainc ()
+              if (! error_state)
+                {
+                  if (b_arg.is_scalar_type ())
+                    {
+                      double b = b_arg.double_value ();
+
+                      if (! error_state)
+                        retval = betaincinv (x, a, b);
+                    }
+                  else
+                    {
+                      Array<double> b = b_arg.array_value ();
+
+                      if (! error_state)
+                        retval = betaincinv (x, a, b);
+                    }
+                }
+            }
+          else
+            {
+              Array<double> 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
+                    {
+                      Array<double> b = b_arg.array_value ();
+
+                      if (! error_state)
+                        retval = betaincinv (x, a, b);
+                    }
+                }
+            }
+        }
+      else
+        {
+          Array<double> x = x_arg.array_value ();
 
-%% test/octave.test/arith/betainc-3.m
-%!error betainc> betainc (1)
+          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
+                    {
+                      Array<double> b = b_arg.array_value ();
+
+                      if (! error_state)
+                        retval = betaincinv (x, a, b);
+                    }
+                }
+            }
+          else
+            {
+              Array<double> 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
+                    {
+                      Array<double> b = b_arg.array_value ();
 
-%% test/octave.test/arith/betainc-4.m
-%!error betainc> betainc (1,2)
+                      if (! error_state)
+                        retval = betaincinv (x, a, b);
+                    }
+                }
+            }
+        }
+
+      // FIXME: It would be better to have an algorithm for betaincinv which accepted
+      //        float inputs and returned float outputs.  As it is, we do extra work
+      //        to calculate betaincinv to double precision and then throw that precision
+      //        away.
+      if (x_arg.is_single_type () || a_arg.is_single_type () ||
+          b_arg.is_single_type ())
+        {
+          retval = Array<float> (retval.array_value ());
+        }
+    }
+  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))
+
+## Test class single as well
+%!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single")))
+%!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single")))
+%!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single")))
+
+## Extreme values
+%!assert (betaincinv (0, 42, 42), 0, sqrt (eps))
+%!assert (betaincinv (1, 42, 42), 1, sqrt (eps))
+
+%!error betaincinv ()
+%!error betaincinv (1, 2)
 */
+
diff --git a/src/mappers.cc b/src/mappers.cc
diff --git a/src/syscalls.cc b/src/syscalls.cc
--- a/src/syscalls.cc
+++ b/src/syscalls.cc
@@ -1594,8 +1594,10 @@
 
 DEFUNX ("canonicalize_file_name", Fcanonicalize_file_name, args, ,
   "-*- texinfo -*-\n\
-@deftypefn {Built-in Function} {[@var{cname}, @var{status}, @var{msg}]} canonicalize_file_name (@var{name})\n\
-Return the canonical name of file @var{name}.\n\
+@deftypefn {Built-in Function} {[@var{cname}, @var{status}, @var{msg}] =} canonicalize_file_name (@var{fname})\n\
+Return the canonical name of file @var{fname}.  If the file does not exist\n\
+the empty string (\"\") is returned.\n\
+@seealso{make_absolute_filename, is_absolute_filename, is_rooted_relative_filename}\n\
 @end deftypefn")
 {
   octave_value_list retval;
diff --git a/src/utils.cc b/src/utils.cc
--- a/src/utils.cc
+++ b/src/utils.cc
@@ -862,8 +862,9 @@
 DEFUN (make_absolute_filename, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} make_absolute_filename (@var{file})\n\
-Return the full name of @var{file}, relative to the current directory.\n\
-@seealso{is_absolute_filename, is_rooted_relative_filename, isdir}\n\
+Return the full name of @var{file} beginning from the root of the file\n\
+system.  No check is done for the existence of @var{file}.\n\
+@seealso{canonicalize_file_name, is_absolute_filename, is_rooted_relative_filename, isdir}\n\
 @end deftypefn")
 {
   octave_value retval = std::string ();