--- a/src/DLD-FUNCTIONS/betainc.cc	Fri Jun 22 16:51:31 2012 +0200
+++ b/src/DLD-FUNCTIONS/betainc.cc	Thu Jun 28 10:18:38 2012 -0700
@@ -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\
@@ -324,3 +327,164 @@
 %!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 ();
+
+              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 ();
+
+          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 ();
+
+                      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)
+*/
+