# HG changeset patch
# User Michele Ginesi <michele.ginesi@gmail.com>
# Date 1504793596 -7200
# Node ID bd89440407aa52ed157da464e41d10830293e0dc
# Parent  451f4f291f460cab6db7cee469556edc88514820
Incomplete beta function moved to a .m file, fixing accuracy and
input validation (see bugs #51157 and #34405). It's inverse also
has been rewritten as .m file.
--
added libinterp/corefcn/__betainc_lentz__.cc
added scripts/specfun/betainc.m
added scripts/specfun/betaincinv.m
changed libinterp/corefcn/module.mk
changed liboctave/external/slatec-fn/module.mk
changed liboctave/numeric/lo-specfun.cc
changed liboctave/numeric/lo-specfun.h
changed scripts/specfun/module.mk
changed scripts/statistics/distributions/betainv.m
changed scripts/statistics/distributions/binocdf.m
removed libinterp/corefcn/betainc.cc
removed liboctave/external/slatec-fn/betai.f
removed liboctave/external/slatec-fn/dbetai.f
removed liboctave/external/slatec-fn/xbetai.f
removed liboctave/external/slatec-fn/xdbetai.f

diff -r 451f4f291f46 -r bd89440407aa libinterp/corefcn/__betainc_lentz__.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__betainc_lentz__.cc	Thu Sep 07 16:13:16 2017 +0200
@@ -0,0 +1,106 @@
+// Copyright (C) 2017 Michele Ginesi
+//
+// This file is part of Octave.
+//
+// Octave is free software; you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 3 of the License, or
+// (at your option) any later version.
+//
+// Octave is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with Octave; see the file COPYING.  If not, see
+// <http://www.gnu.org/licenses/>.
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include "defun.h"
+
+DEFUN (__betainc_lentz__, args, , "Continued fraction for incomplete gamma function")
+{
+  int nargin = args.length ();
+
+  if (nargin != 3)
+    print_usage ();
+  else
+    {
+      octave_value x_arg = args(0);
+      octave_value a_arg = args(1);
+      octave_value b_arg = args(2);
+      if (x_arg.is_single_type () || a_arg.is_single_type () || b_arg.is_single_type ())
+        {
+          const float x = args(0).float_value ();
+          const float a = args(1).float_value ();
+          const float b = args(2).float_value ();
+          static const float tiny = pow (2, -50);
+          static const float eps = std::numeric_limits<float>::epsilon();
+          float f = tiny;
+          float C = f;
+          float D = 0;
+          float alpha_m = 1;
+          float beta_m = a - (a * (a+b)) / (a + 1) * x;
+          float x2 = x * x;
+          float Delta = 0;
+          int m = 1;
+          int maxit = 100;
+          while((std::abs(Delta - 1) > eps) & (m < maxit))
+            {
+               D = beta_m + alpha_m * D;
+               if (D == 0)
+                 D = tiny;
+               C = beta_m + alpha_m / C;
+               if (C == 0)
+                 C = tiny;
+               D = 1 / D;
+               Delta = C * D;
+               f *= Delta;
+               alpha_m = ((a + m - 1) * (a + b + m - 1) * (b - m) * m) / ((a + 2 * m - 1) * (a + 2 * m - 1)) * x2;
+               beta_m = a + 2 * m + ((m * (b - m)) / (a + 2 * m - 1) - ((a + m) * (a + b + m)) / (a + 2 * m + 1)) * x;
+               m++;
+             }
+           if (! error_state)
+            return octave_value (f);
+        }
+      else
+        {
+          const double x = args(0).double_value ();
+          const double a = args(1).double_value ();
+          const double b = args(2).double_value ();
+          static const double tiny = pow (2, -100);
+          static const double eps = std::numeric_limits<double>::epsilon();
+          double f = tiny;
+          double C = f;
+          double D = 0;
+          double alpha_m = 1;
+          double beta_m = a - (a * (a+b)) / (a + 1) * x;
+          double x2 = x * x;
+          double Delta = 0;
+          int m = 1;
+          int maxit = 200;
+          while((std::abs(Delta - 1) > eps) & (m < maxit))
+            {
+               D = beta_m + alpha_m * D;
+               if (D == 0)
+                 D = tiny;
+               C = beta_m + alpha_m / C;
+               if (C == 0)
+                 C = tiny;
+               D = 1 / D;
+               Delta = C * D;
+               f *= Delta;
+               alpha_m = ((a + m - 1) * (a + b + m - 1) * (b - m) * m) / ((a + 2 * m - 1) * (a + 2 * m - 1)) * x2;
+               beta_m = a + 2 * m + ((m * (b - m)) / (a + 2 * m - 1) - ((a + m) * (a + b + m)) / (a + 2 * m + 1)) * x;
+               m++;
+             }
+           if (! error_state)
+            return octave_value (f);
+        }
+      }
+  return octave_value_list ();
+}
diff -r 451f4f291f46 -r bd89440407aa libinterp/corefcn/betainc.cc
--- a/libinterp/corefcn/betainc.cc	Thu Sep 07 15:58:26 2017 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,424 +0,0 @@
-/*
-
-Copyright (C) 1997-2017 John W. Eaton
-
-This file is part of Octave.
-
-Octave is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-Octave is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<https://www.gnu.org/licenses/>.
-
-*/
-
-#if defined (HAVE_CONFIG_H)
-#  include "config.h"
-#endif
-
-#include "lo-specfun.h"
-
-#include "defun.h"
-#include "error.h"
-#include "errwarn.h"
-#include "ovl.h"
-#include "utils.h"
-
-
-DEFUN (betainc, args, ,
-       doc: /* -*- texinfo -*-
-@deftypefn {} {} betainc (@var{x}, @var{a}, @var{b})
-Compute the regularized incomplete Beta function.
-
-The regularized incomplete Beta function is defined by
-@tex
-$$
- I (x, a, b) = {1 \over {B (a, b)}} \int_0^x t^{(a-z)} (1-t)^{(b-1)} dt.
-$$
-@end tex
-@ifnottex
-@c Set example in small font to prevent overfull line
-
-@smallexample
-@group
-                                   x
-                          1       /
-betainc (x, a, b) = -----------   | t^(a-1) (1-t)^(b-1) dt.
-                    beta (a, b)   /
-                               t=0
-@end group
-@end smallexample
-
-@end ifnottex
-
-If @var{x} has more than one component, both @var{a} and @var{b} must be
-scalars.  If @var{x} is a scalar, @var{a} and @var{b} must be of
-compatible dimensions.
-@seealso{betaincinv, beta, betaln}
-@end deftypefn */)
-{
-  if (args.length () != 3)
-    print_usage ();
-
-  octave_value retval;
-
-  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 (b_arg.is_scalar_type ())
-                {
-                  float b = b_arg.float_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<float> b = b_arg.float_array_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-            }
-          else
-            {
-              Array<float> a = a_arg.float_array_value ();
-
-              if (b_arg.is_scalar_type ())
-                {
-                  float b = b_arg.float_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<float> b = b_arg.float_array_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-            }
-        }
-      else
-        {
-          Array<float> x = x_arg.float_array_value ();
-
-          if (a_arg.is_scalar_type ())
-            {
-              float a = a_arg.float_value ();
-
-              if (b_arg.is_scalar_type ())
-                {
-                  float b = b_arg.float_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<float> b = b_arg.float_array_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-            }
-          else
-            {
-              Array<float> a = a_arg.float_array_value ();
-
-              if (b_arg.is_scalar_type ())
-                {
-                  float b = b_arg.float_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<float> b = b_arg.float_array_value ();
-
-                  retval = octave::math::betainc (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 (b_arg.is_scalar_type ())
-                {
-                  double b = b_arg.double_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<double> b = b_arg.array_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-            }
-          else
-            {
-              Array<double> a = a_arg.array_value ();
-
-              if (b_arg.is_scalar_type ())
-                {
-                  double b = b_arg.double_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<double> b = b_arg.array_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-            }
-        }
-      else
-        {
-          Array<double> x = x_arg.array_value ();
-
-          if (a_arg.is_scalar_type ())
-            {
-              double a = a_arg.double_value ();
-
-              if (b_arg.is_scalar_type ())
-                {
-                  double b = b_arg.double_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<double> b = b_arg.array_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-            }
-          else
-            {
-              Array<double> a = a_arg.array_value ();
-
-              if (b_arg.is_scalar_type ())
-                {
-                  double b = b_arg.double_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-              else
-                {
-                  Array<double> b = b_arg.array_value ();
-
-                  retval = octave::math::betainc (x, a, b);
-                }
-            }
-        }
-    }
-
-  return retval;
-}
-
-/*
-## Double precision
-%!test
-%! a = [1, 1.5, 2, 3];
-%! b = [4, 3, 2, 1];
-%! v1 = betainc (1,a,b);
-%! v2 = [1,1,1,1];
-%! x = [.2, .4, .6, .8];
-%! v3 = betainc (x, a, b);
-%! v4 = 1 - betainc (1.-x, b, a);
-%! assert (v1, v2, sqrt (eps));
-%! assert (v3, v4, sqrt (eps));
-
-## Single precision
-%!test
-%! a = single ([1, 1.5, 2, 3]);
-%! b = single ([4, 3, 2, 1]);
-%! v1 = betainc (1,a,b);
-%! v2 = single ([1,1,1,1]);
-%! x = single ([.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")));
-
-## Mixed double/single precision
-%!test
-%! a = single ([1, 1.5, 2, 3]);
-%! b = [4, 3, 2, 1];
-%! v1 = betainc (1,a,b);
-%! v2 = single ([1,1,1,1]);
-%! 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")));
-
-%!error betainc ()
-%!error betainc (1)
-%!error betainc (1,2)
-%!error betainc (1,2,3,4)
-*/
-
-DEFUN (betaincinv, args, ,
-       doc: /* -*- texinfo -*-
-@deftypefn {} {} betaincinv (@var{y}, @var{a}, @var{b})
-Compute the inverse of the incomplete Beta function.
-
-The inverse is the value @var{x} such that
-
-@example
-@var{y} == betainc (@var{x}, @var{a}, @var{b})
-@end example
-@seealso{betainc, beta, betaln}
-@end deftypefn */)
-{
-  if (args.length () != 3)
-    print_usage ();
-
-  octave_value retval;
-
-  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 (b_arg.is_scalar_type ())
-            {
-              double b = b_arg.double_value ();
-
-              retval = octave::math::betaincinv (x, a, b);
-            }
-          else
-            {
-              Array<double> b = b_arg.array_value ();
-
-              retval = octave::math::betaincinv (x, a, b);
-            }
-        }
-      else
-        {
-          Array<double> a = a_arg.array_value ();
-
-          if (b_arg.is_scalar_type ())
-            {
-              double b = b_arg.double_value ();
-
-              retval = octave::math::betaincinv (x, a, b);
-            }
-          else
-            {
-              Array<double> b = b_arg.array_value ();
-
-              retval = octave::math::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 (b_arg.is_scalar_type ())
-            {
-              double b = b_arg.double_value ();
-
-              retval = octave::math::betaincinv (x, a, b);
-            }
-          else
-            {
-              Array<double> b = b_arg.array_value ();
-
-              retval = octave::math::betaincinv (x, a, b);
-            }
-        }
-      else
-        {
-          Array<double> a = a_arg.array_value ();
-
-          if (b_arg.is_scalar_type ())
-            {
-              double b = b_arg.double_value ();
-
-              retval = octave::math::betaincinv (x, a, b);
-            }
-          else
-            {
-              Array<double> b = b_arg.array_value ();
-
-              retval = octave::math::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 ());
-    }
-
-  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 -r 451f4f291f46 -r bd89440407aa libinterp/corefcn/module.mk
--- a/libinterp/corefcn/module.mk	Thu Sep 07 15:58:26 2017 +0200
+++ b/libinterp/corefcn/module.mk	Thu Sep 07 16:13:16 2017 +0200
@@ -107,6 +107,7 @@
 
 COREFCN_SRC = \
   %reldir%/Cell.cc \
+  %reldir%/__betainc_lentz__.cc \
   %reldir%/__contourc__.cc \
   %reldir%/__dsearchn__.cc \
   %reldir%/__gammainc_lentz__.cc \
@@ -119,7 +120,6 @@
   %reldir%/__qp__.cc \
   %reldir%/balance.cc \
   %reldir%/besselj.cc \
-  %reldir%/betainc.cc \
   %reldir%/bitfcns.cc \
   %reldir%/bsxfun.cc \
   %reldir%/c-file-ptr-stream.cc \
diff -r 451f4f291f46 -r bd89440407aa liboctave/external/slatec-fn/betai.f
--- a/liboctave/external/slatec-fn/betai.f	Thu Sep 07 15:58:26 2017 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,118 +0,0 @@
-*DECK BETAI
-      REAL FUNCTION BETAI (X, PIN, QIN)
-C***BEGIN PROLOGUE  BETAI
-C***PURPOSE  Calculate the incomplete Beta function.
-C***LIBRARY   SLATEC (FNLIB)
-C***CATEGORY  C7F
-C***TYPE      SINGLE PRECISION (BETAI-S, DBETAI-D)
-C***KEYWORDS  FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS
-C***AUTHOR  Fullerton, W., (LANL)
-C***DESCRIPTION
-C
-C   BETAI calculates the REAL incomplete beta function.
-C
-C   The incomplete beta function ratio is the probability that a
-C   random variable from a beta distribution having parameters PIN and
-C   QIN will be less than or equal to X.
-C
-C     -- Input Arguments -- All arguments are REAL.
-C   X      upper limit of integration.  X must be in (0,1) inclusive.
-C   PIN    first beta distribution parameter.  PIN must be .GT. 0.0.
-C   QIN    second beta distribution parameter.  QIN must be .GT. 0.0.
-C
-C***REFERENCES  Nancy E. Bosten and E. L. Battiste, Remark on Algorithm
-C                 179, Communications of the ACM 17, 3 (March 1974),
-C                 pp. 156.
-C***ROUTINES CALLED  ALBETA, R1MACH, XERMSG
-C***REVISION HISTORY  (YYMMDD)
-C   770401  DATE WRITTEN
-C   890531  Changed all specific intrinsics to generic.  (WRB)
-C   890531  REVISION DATE from Version 3.2
-C   891214  Prologue converted to Version 4.0 format.  (BAB)
-C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
-C   900326  Removed duplicate information from DESCRIPTION section.
-C           (WRB)
-C   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
-C***END PROLOGUE  BETAI
-      LOGICAL FIRST
-      SAVE EPS, ALNEPS, SML, ALNSML, FIRST
-      DATA FIRST /.TRUE./
-C***FIRST EXECUTABLE STATEMENT  BETAI
-      IF (FIRST) THEN
-         EPS = R1MACH(3)
-         ALNEPS = LOG(EPS)
-         SML = R1MACH(1)
-         ALNSML = LOG(SML)
-      ENDIF
-      FIRST = .FALSE.
-C
-      IF (X .LT. 0. .OR. X .GT. 1.0) CALL XERMSG ('SLATEC', 'BETAI',
-     +   'X IS NOT IN THE RANGE (0,1)', 1, 2)
-      IF (PIN .LE. 0. .OR. QIN .LE. 0.) CALL XERMSG ('SLATEC', 'BETAI',
-     +   'P AND/OR Q IS LE ZERO', 2, 2)
-C
-      Y = X
-      P = PIN
-      Q = QIN
-      IF (Q.LE.P .AND. X.LT.0.8) GO TO 20
-      IF (X.LT.0.2) GO TO 20
-      Y = 1.0 - Y
-      P = QIN
-      Q = PIN
-C
- 20   IF ((P+Q)*Y/(P+1.).LT.EPS) GO TO 80
-C
-C EVALUATE THE INFINITE SUM FIRST.
-C TERM WILL EQUAL Y**P/BETA(PS,P) * (1.-PS)I * Y**I / FAC(I)
-C
-      PS = Q - AINT(Q)
-      IF (PS.EQ.0.) PS = 1.0
-      XB = P*LOG(Y) -  ALBETA(PS, P) - LOG(P)
-      BETAI = 0.0
-      IF (XB.LT.ALNSML) GO TO 40
-C
-      BETAI = EXP (XB)
-      TERM = BETAI*P
-      IF (PS.EQ.1.0) GO TO 40
-C
-      N = MAX (ALNEPS/LOG(Y), 4.0E0)
-      DO 30 I=1,N
-        TERM = TERM*(I-PS)*Y/I
-        BETAI = BETAI + TERM/(P+I)
- 30   CONTINUE
-C
-C NOW EVALUATE THE FINITE SUM, MAYBE.
-C
- 40   IF (Q.LE.1.0) GO TO 70
-C
-      XB = P*LOG(Y) + Q*LOG(1.0-Y) - ALBETA(P,Q) - LOG(Q)
-      IB = MAX (XB/ALNSML, 0.0E0)
-      TERM = EXP (XB - IB*ALNSML)
-      C = 1.0/(1.0-Y)
-      P1 = Q*C/(P+Q-1.)
-C
-      FINSUM = 0.0
-      N = Q
-      IF (Q.EQ.REAL(N)) N = N - 1
-      DO 50 I=1,N
-        IF (P1.LE.1.0 .AND. TERM/EPS.LE.FINSUM) GO TO 60
-        TERM = (Q-I+1)*C*TERM/(P+Q-I)
-C
-        IF (TERM.GT.1.0) IB = IB - 1
-        IF (TERM.GT.1.0) TERM = TERM*SML
-C
-        IF (IB.EQ.0) FINSUM = FINSUM + TERM
- 50   CONTINUE
-C
- 60   BETAI = BETAI + FINSUM
- 70   IF (Y.NE.X .OR. P.NE.PIN) BETAI = 1.0 - BETAI
-      BETAI = MAX (MIN (BETAI, 1.0), 0.0)
-      RETURN
-C
- 80   BETAI = 0.0
-      XB = P*LOG(MAX(Y,SML)) - LOG(P) - ALBETA(P,Q)
-      IF (XB.GT.ALNSML .AND. Y.NE.0.) BETAI = EXP (XB)
-      IF (Y.NE.X .OR. P.NE.PIN) BETAI = 1.0 - BETAI
-      RETURN
-C
-      END
diff -r 451f4f291f46 -r bd89440407aa liboctave/external/slatec-fn/dbetai.f
--- a/liboctave/external/slatec-fn/dbetai.f	Thu Sep 07 15:58:26 2017 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,121 +0,0 @@
-
-*DECK DBETAI
-      DOUBLE PRECISION FUNCTION DBETAI (X, PIN, QIN)
-C***BEGIN PROLOGUE  DBETAI
-C***PURPOSE  Calculate the incomplete Beta function.
-C***LIBRARY   SLATEC (FNLIB)
-C***CATEGORY  C7F
-C***TYPE      DOUBLE PRECISION (BETAI-S, DBETAI-D)
-C***KEYWORDS  FNLIB, INCOMPLETE BETA FUNCTION, SPECIAL FUNCTIONS
-C***AUTHOR  Fullerton, W., (LANL)
-C***DESCRIPTION
-C
-C   DBETAI calculates the DOUBLE PRECISION incomplete beta function.
-C
-C   The incomplete beta function ratio is the probability that a
-C   random variable from a beta distribution having parameters PIN and
-C   QIN will be less than or equal to X.
-C
-C     -- Input Arguments -- All arguments are DOUBLE PRECISION.
-C   X      upper limit of integration.  X must be in (0,1) inclusive.
-C   PIN    first beta distribution parameter.  PIN must be .GT. 0.0.
-C   QIN    second beta distribution parameter.  QIN must be .GT. 0.0.
-C
-C***REFERENCES  Nancy E. Bosten and E. L. Battiste, Remark on Algorithm
-C                 179, Communications of the ACM 17, 3 (March 1974),
-C                 pp. 156.
-C***ROUTINES CALLED  D1MACH, DLBETA, XERMSG
-C***REVISION HISTORY  (YYMMDD)
-C   770701  DATE WRITTEN
-C   890531  Changed all specific intrinsics to generic.  (WRB)
-C   890911  Removed unnecessary intrinsics.  (WRB)
-C   890911  REVISION DATE from Version 3.2
-C   891214  Prologue converted to Version 4.0 format.  (BAB)
-C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
-C   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
-C***END PROLOGUE  DBETAI
-      DOUBLE PRECISION X, PIN, QIN, ALNEPS, ALNSML, C, EPS, FINSUM, P,
-     1  PS, Q, SML, TERM, XB, XI, Y, D1MACH, DLBETA, P1
-      LOGICAL FIRST
-      SAVE EPS, ALNEPS, SML, ALNSML, FIRST
-      DATA FIRST /.TRUE./
-C***FIRST EXECUTABLE STATEMENT  DBETAI
-      IF (FIRST) THEN
-         EPS = D1MACH(3)
-         ALNEPS = LOG (EPS)
-         SML = D1MACH(1)
-         ALNSML = LOG (SML)
-      ENDIF
-      FIRST = .FALSE.
-C
-      IF (X .LT. 0.D0 .OR. X .GT. 1.D0) CALL XERMSG ('SLATEC', 'DBETAI',
-     +   'X IS NOT IN THE RANGE (0,1)', 1, 2)
-      IF (PIN .LE. 0.D0 .OR. QIN .LE. 0.D0) CALL XERMSG ('SLATEC',
-     +   'DBETAI', 'P AND/OR Q IS LE ZERO', 2, 2)
-C
-      Y = X
-      P = PIN
-      Q = QIN
-      IF (Q.LE.P .AND. X.LT.0.8D0) GO TO 20
-      IF (X.LT.0.2D0) GO TO 20
-      Y = 1.0D0 - Y
-      P = QIN
-      Q = PIN
-C
- 20   IF ((P+Q)*Y/(P+1.D0).LT.EPS) GO TO 80
-C
-C EVALUATE THE INFINITE SUM FIRST.  TERM WILL EQUAL
-C Y**P/BETA(PS,P) * (1.-PS)-SUB-I * Y**I / FAC(I) .
-C
-      PS = Q - AINT(Q)
-      IF (PS.EQ.0.D0) PS = 1.0D0
-      XB = P*LOG(Y) - DLBETA(PS,P) - LOG(P)
-      DBETAI = 0.0D0
-      IF (XB.LT.ALNSML) GO TO 40
-C
-      DBETAI = EXP (XB)
-      TERM = DBETAI*P
-      IF (PS.EQ.1.0D0) GO TO 40
-      N = MAX (ALNEPS/LOG(Y), 4.0D0)
-      DO 30 I=1,N
-        XI = I
-        TERM = TERM * (XI-PS)*Y/XI
-        DBETAI = DBETAI + TERM/(P+XI)
- 30   CONTINUE
-C
-C NOW EVALUATE THE FINITE SUM, MAYBE.
-C
- 40   IF (Q.LE.1.0D0) GO TO 70
-C
-      XB = P*LOG(Y) + Q*LOG(1.0D0-Y) - DLBETA(P,Q) - LOG(Q)
-      IB = MAX (XB/ALNSML, 0.0D0)
-      TERM = EXP(XB - IB*ALNSML)
-      C = 1.0D0/(1.D0-Y)
-      P1 = Q*C/(P+Q-1.D0)
-C
-      FINSUM = 0.0D0
-      N = Q
-      IF (Q.EQ.DBLE(N)) N = N - 1
-      DO 50 I=1,N
-        IF (P1.LE.1.0D0 .AND. TERM/EPS.LE.FINSUM) GO TO 60
-        XI = I
-        TERM = (Q-XI+1.0D0)*C*TERM/(P+Q-XI)
-C
-        IF (TERM.GT.1.0D0) IB = IB - 1
-        IF (TERM.GT.1.0D0) TERM = TERM*SML
-C
-        IF (IB.EQ.0) FINSUM = FINSUM + TERM
- 50   CONTINUE
-C
- 60   DBETAI = DBETAI + FINSUM
- 70   IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI
-      DBETAI = MAX (MIN (DBETAI, 1.0D0), 0.0D0)
-      RETURN
-C
- 80   DBETAI = 0.0D0
-      XB = P*LOG(MAX(Y,SML)) - LOG(P) - DLBETA(P,Q)
-      IF (XB.GT.ALNSML .AND. Y.NE.0.0D0) DBETAI = EXP(XB)
-      IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI
-C
-      RETURN
-      END
diff -r 451f4f291f46 -r bd89440407aa liboctave/external/slatec-fn/module.mk
--- a/liboctave/external/slatec-fn/module.mk	Thu Sep 07 15:58:26 2017 +0200
+++ b/liboctave/external/slatec-fn/module.mk	Thu Sep 07 16:13:16 2017 +0200
@@ -3,13 +3,11 @@
   %reldir%/alngam.f \
   %reldir%/alnrel.f \
   %reldir%/algams.f \
-  %reldir%/betai.f \
   %reldir%/csevl.f \
   %reldir%/d9gmit.f \
   %reldir%/d9lgic.f \
   %reldir%/d9lgit.f \
   %reldir%/d9lgmc.f \
-  %reldir%/dbetai.f \
   %reldir%/dcsevl.f \
   %reldir%/dgamlm.f \
   %reldir%/dgamma.f \
@@ -32,8 +30,6 @@
   %reldir%/r9lgmc.f \
   %reldir%/r9lgit.f \
   %reldir%/r9gmit.f \
-  %reldir%/r9lgic.f \
-  %reldir%/xdbetai.f \
-  %reldir%/xbetai.f
+  %reldir%/r9lgic.f
 
 DIRSTAMP_FILES += %reldir%/$(octave_dirstamp)
diff -r 451f4f291f46 -r bd89440407aa liboctave/external/slatec-fn/xbetai.f
--- a/liboctave/external/slatec-fn/xbetai.f	Thu Sep 07 15:58:26 2017 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-      subroutine xbetai (x, a, b, result)
-      external betai
-      real x, a, b, result, betai
-      result = betai (x, a, b)
-      return
-      end
diff -r 451f4f291f46 -r bd89440407aa liboctave/external/slatec-fn/xdbetai.f
--- a/liboctave/external/slatec-fn/xdbetai.f	Thu Sep 07 15:58:26 2017 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-      subroutine xdbetai (x, a, b, result)
-      external dbetai
-      double precision x, a, b, result, dbetai
-      result = dbetai (x, a, b)
-      return
-      end
diff -r 451f4f291f46 -r bd89440407aa liboctave/numeric/lo-specfun.cc
--- a/liboctave/numeric/lo-specfun.cc	Thu Sep 07 15:58:26 2017 +0200
+++ b/liboctave/numeric/lo-specfun.cc	Thu Sep 07 16:13:16 2017 +0200
@@ -1374,765 +1374,6 @@
 #undef NN_BESSEL
 #undef RC_BESSEL
 
-    OCTAVE_NORETURN static void
-    err_betainc_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)
-        ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
-         d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
-    }
-
-    OCTAVE_NORETURN static void
-    err_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)
-    {
-      double retval;
-      F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval));
-      return retval;
-    }
-
-    Array<double>
-    betainc (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++ = betainc (x, a, b(i));
-
-      return retval;
-    }
-
-    Array<double>
-    betainc (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++ = betainc (x, a(i), b);
-
-      return retval;
-    }
-
-    Array<double>
-    betainc (double x, const Array<double>& a, const Array<double>& b)
-    {
-      Array<double> retval;
-      dim_vector dv = a.dims ();
-
-      if (dv != b.dims ())
-        err_betainc_nonconformant (dim_vector (0, 0), 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++ = betainc (x, a(i), b(i));
-
-      return retval;
-    }
-
-    Array<double>
-    betainc (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++ = betainc (x(i), a, b);
-
-      return retval;
-    }
-
-    Array<double>
-    betainc (const Array<double>& x, double a, const Array<double>& b)
-    {
-      Array<double> retval;
-      dim_vector dv = x.dims ();
-
-      if (dv != b.dims ())
-        err_betainc_nonconformant (dv, dim_vector (0, 0), 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++ = betainc (x(i), a, b(i));
-
-      return retval;
-    }
-
-    Array<double>
-    betainc (const Array<double>& x, const Array<double>& a, double b)
-    {
-      Array<double> retval;
-      dim_vector dv = x.dims ();
-
-      if (dv != a.dims ())
-        err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
-
-      octave_idx_type nel = dv.numel ();
-
-      retval.resize (dv);
-
-      double *pretval = retval.fortran_vec ();
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        *pretval++ = betainc (x(i), a(i), b);
-
-      return retval;
-    }
-
-    Array<double>
-    betainc (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 ())
-        err_betainc_nonconformant (dv, a.dims (), 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++ = betainc (x(i), a(i), b(i));
-
-      return retval;
-    }
-
-    float
-    betainc (float x, float a, float b)
-    {
-      float retval;
-      F77_XFCN (xbetai, XBETAI, (x, a, b, retval));
-      return retval;
-    }
-
-    Array<float>
-    betainc (float x, float a, const Array<float>& 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;
-    }
-
-    Array<float>
-    betainc (float x, const Array<float>& a, float b)
-    {
-      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;
-    }
-
-    Array<float>
-    betainc (float x, const Array<float>& a, const Array<float>& b)
-    {
-      Array<float> retval;
-      dim_vector dv = a.dims ();
-
-      if (dv != b.dims ())
-        err_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ());
-
-      octave_idx_type nel = dv.numel ();
-
-      retval.resize (dv);
-
-      float *pretval = retval.fortran_vec ();
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        *pretval++ = betainc (x, a(i), b(i));
-
-      return retval;
-    }
-
-    Array<float>
-    betainc (const Array<float>& x, float a, float b)
-    {
-      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;
-    }
-
-    Array<float>
-    betainc (const Array<float>& x, float a, const Array<float>& b)
-    {
-      Array<float> retval;
-      dim_vector dv = x.dims ();
-
-      if (dv != b.dims ())
-        err_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ());
-
-      octave_idx_type nel = dv.numel ();
-
-      retval.resize (dv);
-
-      float *pretval = retval.fortran_vec ();
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        *pretval++ = betainc (x(i), a, b(i));
-
-      return retval;
-    }
-
-    Array<float>
-    betainc (const Array<float>& x, const Array<float>& a, float b)
-    {
-      Array<float> retval;
-      dim_vector dv = x.dims ();
-
-      if (dv != a.dims ())
-        err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
-
-      octave_idx_type nel = dv.numel ();
-
-      retval.resize (dv);
-
-      float *pretval = retval.fortran_vec ();
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        *pretval++ = betainc (x(i), a(i), b);
-
-      return retval;
-    }
-
-    Array<float>
-    betainc (const Array<float>& x, const Array<float>& a,
-             const Array<float>& b)
-    {
-      Array<float> retval;
-      dim_vector dv = x.dims ();
-
-      if (dv != a.dims () || dv != b.dims ())
-        err_betainc_nonconformant (dv, a.dims (), b.dims ());
-
-      octave_idx_type nel = dv.numel ();
-
-      retval.resize (dv);
-
-      float *pretval = retval.fortran_vec ();
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        *pretval++ = betainc (x(i), a(i), b(i));
-
-      return retval;
-    }
-
-    //
-    //  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 = static_cast<int> (qq + cx * psq);
-
-      //  Use the Soper reduction formula.
-
-      rx = xx / cx;
-      temp = qq - ai;
-      if (ns == 0)
-        {
-          rx = xx;
-        }
-
-      for ( ; ; )
-        {
-          term *= temp * rx / (pp + ai);
-          value += term;
-          temp = fabs (term);
-
-          if (temp <= acu && temp <= acu * value)
-            {
-              value *= exp (pp * std::log (xx)
-                            + (qq - 1.0) * std::log (cx) - beta) / pp;
-
-              if (indx)
-                {
-                  value = 1.0 - value;
-                }
-              break;
-            }
-
-          ai += 1.0;
-          ns -= 1;
-
-          if (0 <= ns)
-            {
-              temp = qq - ai;
-              if (ns == 0)
-                {
-                  rx = xx;
-                }
-            }
-          else
-            {
-              temp = 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 = std::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 = std::sqrt (- std::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 * std::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 * std::pow (1.0 - t + ycur * std::sqrt (t), 3);
-
-          if (t <= 0.0)
-            {
-              value = 1.0 - exp ((std::log ((1.0 - a) * qq) + beta) / qq);
-            }
-          else
-            {
-              t = (4.0 * pp + r - 2.0) / t;
-
-              if (t <= 1.0)
-                {
-                  value = exp ((std::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 / std::pow (a, 0.2) - 13.0, sae);
-
-      acu = std::pow (10.0, iex);
-
-      for ( ; ; )
-        {
-          ycur = betain (value, pp, qq, beta, err);
-
-          if (err)
-            {
-              return value;
-            }
-
-          xin = value;
-          ycur = (ycur - a) * exp (beta + r * std::log (xin)
-                                   + t * std::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 /= 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 /= 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 ())
-        err_betaincinv_nonconformant (dim_vector (0, 0), 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));
-
-      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 ())
-        err_betaincinv_nonconformant (dv, dim_vector (0, 0), 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));
-
-      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 ())
-        err_betaincinv_nonconformant (dv, a.dims (), dim_vector (0, 0));
-
-      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);
-
-      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 ())
-        err_betaincinv_nonconformant (dv, a.dims (), 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));
-
-      return retval;
-    }
-
     Complex
     biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr)
     {
@@ -3132,6 +2373,7 @@
 FloatComplexNDArray airy (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::airy (z, deriv, scaled, ierr); }
 FloatComplexNDArray biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::biry (z, deriv, scaled, ierr); }
 
+<<<<<<< dest
 Array<double> betainc (double x, double a, const Array<double>& b) { return octave::math::betainc (x, a, b); }
 Array<double> betainc (double x, const Array<double>& a, double b) { return octave::math::betainc (x, a, b); }
 Array<double> betainc (double x, const Array<double>& a, const Array<double>& b) { return octave::math::betainc (x, a, b); }
diff -r 451f4f291f46 -r bd89440407aa liboctave/numeric/lo-specfun.h
--- a/liboctave/numeric/lo-specfun.h	Thu Sep 07 15:58:26 2017 +0200
+++ b/liboctave/numeric/lo-specfun.h	Thu Sep 07 16:13:16 2017 +0200
@@ -282,55 +282,6 @@
     extern OCTAVE_API FloatComplexMatrix besselh2 (const FloatRowVector& alpha,
                                                    const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr);
 
-    extern OCTAVE_API double betainc (double x, double a, double 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,
-                                             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 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,
-                                            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 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,
-                                                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);
-
     extern OCTAVE_API Complex biry (const Complex& z, bool deriv, bool scaled,
                                     octave_idx_type& ierr);
     extern OCTAVE_API ComplexMatrix biry (const ComplexMatrix& z, bool deriv,
@@ -771,42 +722,6 @@
 OCTAVE_DEPRECATED (4.2, "use 'octave::math::biry' instead")
 extern OCTAVE_API FloatComplexNDArray biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr);
 
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-inline double betainc (double x, double a, double b) { return octave::math::betainc (x, a, b); }
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<double> betainc (double x, double a, const Array<double>& b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, double b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<double> betainc (double x, const Array<double>& a, const Array<double>& b);
-
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-inline float betainc (float x, float a, float b) { return octave::math::betainc (x, a, b); }
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<double> betainc (const Array<double>& x, double a, const Array<double>& b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, double b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<double> betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b);
-
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API float betainc (float x, float a, float b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<float> betainc (float x, float a, const Array<float>& b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, float b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<float> betainc (float x, const Array<float>& a, const Array<float>& b);
-
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, float b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<float> betainc (const Array<float>& x, float a, const Array<float>& b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, float b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betainc' instead")
-extern OCTAVE_API Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b);
-
 OCTAVE_DEPRECATED (4.2, "use 'octave::math::gammainc' instead")
 inline double gammainc (double x, double a, bool& err) { return octave::math::gammainc (x, a, err); }
 OCTAVE_DEPRECATED (4.2, "use 'octave::math::gammainc' instead")
@@ -887,27 +802,6 @@
 OCTAVE_DEPRECATED (4.2, "use 'octave::math::dawson' instead")
 inline FloatComplex dawson (const FloatComplex& x) { return octave::math::dawson (x); }
 
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-inline double betaincinv (double x, double a, double b) { return octave::math::betaincinv (x, a, b); }
-
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API double betaincinv (double x, double a, double b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API Array<double> betaincinv (double x, double a, const Array<double>& b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, double b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b);
-
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, const Array<double>& b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, double b);
-OCTAVE_DEPRECATED (4.2, "use 'octave::math::betaincinv' instead")
-extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b);
-
 OCTAVE_DEPRECATED (4.2, "use 'octave::math::ellipj' instead")
 inline void ellipj (double u, double m, double& sn, double& cn, double& dn, double& err) { octave::math::ellipj (u, m, sn, cn, dn, err); }
 OCTAVE_DEPRECATED (4.2, "use 'octave::math::ellipj' instead")
diff -r 451f4f291f46 -r bd89440407aa scripts/specfun/betainc.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/specfun/betainc.m	Thu Sep 07 16:13:16 2017 +0200
@@ -0,0 +1,231 @@
+## Copyright (C) 2017 Michele Ginesi
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## Authors: Michele Ginesi <michele.ginesi@gmail.com>
+
+## -*- texinfo -*-
+## @deftypefn {} {} betainc (@var{x}, @var{a}, @var{b})
+## @deftypefnx {} {} betainc (@var{x}, @var{a}, @var{b}, @var{tail})
+## Compute the incomplete beta function ratio.
+##
+## This is defined as
+## @tex
+## $$
+## I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt}
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##                x
+##               /
+##              |
+## I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt
+##              |
+##              /
+##             0
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## with @var{x} real in [0,1], @var{a} and @var{b} real and strictly positive.
+## If one of the input has more than one components, then the others must be
+## scalar or of compatible dimensions.
+##
+## By default or if @var{tail} is @qcode{"lower"} the incomplete beta function
+## ratio integrated from 0 to @var{x} is computed. If @var{tail} is
+## @qcode{"upper"} then the complementary function integrated from @var{x} to
+## 1 is calculated. The two choices are related as
+##
+## betainc (@var{x}, @var{a}, @var{b}, @qcode{"lower"}) = 
+## 1 - betainc (@var{x}, @var{a}, @var{b}, @qcode{"upper"}).
+##
+## Reference
+##
+## @nospell{A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland, W.B. Jones}
+## @cite{Handbook of Continued Fractions for Special Functions},
+## ch. 18.
+##
+## @seealso{beta, betaincinv, betaln}
+##
+## @end deftypefn
+
+function [y] = betainc (x, a, b, tail = "lower")
+
+  if (nargin > 4 || nargin < 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (x) || ! isscalar (a) || ! isscalar (b))
+    [err, x, a, b] = common_size (x, a, b);
+    if (err > 0)
+      error ("betainc: x, a and b must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (x) || iscomplex (a) || iscomplex (b))
+    error ("betainc: inputs must be real or integer");
+  endif
+
+  if (any (a <= 0))
+    error ("betainc: a must be strictly positive");
+  endif
+
+    if (any (b <= 0))
+    error ("betainc: b must be strictly positive");
+  endif
+
+  if (any (x > 1 | x < 0))
+    error ("betainc: x must be between 0 and 1");
+  endif
+
+  if (isinteger (x))
+    y = double (x);
+  endif
+
+  if (isinteger (a))
+    a = double (a);
+  endif
+
+  if (isinteger (b))
+    b = double (b);
+  endif
+
+  y = zeros (size (x));
+
+  # If any of the argument is single, also the output should be
+
+  if (strcmpi (class (y), "single") || strcmpi (class (a), "single") || strcmpi (class (b), "single"))
+    a = single (a);
+    b = single (b);
+    x = single (x);
+    y = single (y);
+  endif
+
+  # In the following, we use the fact that the continued fraction we will
+  # use is more efficient when x <= a / (a+b).
+  # Moreover, to compute the upper version, which is defined as
+  # I_x(a,b,"upper") = 1 - I_x(a,b) we used the property
+  # I_x(a,b) + I_(1-x) (b,a) = 1.
+
+  if (strcmpi (tail, "upper"))
+    flag = (x < a./(a+b));
+    x(!flag) = 1 - x(!flag);
+    [a(!flag), b(!flag)] = deal (b(!flag), a(!flag));
+  elseif (strcmpi (tail, "lower"))
+    flag = (x > a./(a+b));
+    x (flag) = 1 - x(flag);
+    [a(flag), b(flag)] = deal (b(flag), a(flag));
+  else
+    error ("betainc: invalid value for flag")
+  endif
+
+  f = zeros (size (x), class (x));
+
+  ## Continued fractions: CPVWJ, formula 18.5.20, modified Lentz algorithm
+  ## implemented in a separate .cc file. This particular continued fraction
+  ## gives (B(a,b) * I_x(a,b)) / (x^a * (1-x)^b).
+
+  for i = 1 : numel (x)
+    f(i) = __betainc_lentz__ (x(i), a(i), b(i));
+  endfor
+  # We divide the continued fraction by B(a,b) / (x^a * (1-x)^b)
+  # to obtain I_x(a,b).
+  y = a .* log (x) + b .* log1p (-x) + gammaln (a + b) - ...
+    gammaln (a) - gammaln (b) + log (f);
+  y = real (exp (y));
+  y(flag) = 1 - y(flag);
+
+endfunction
+
+## Tests from betainc.cc
+
+# Double precision
+%!test
+%! a = [1, 1.5, 2, 3];
+%! b = [4, 3, 2, 1];
+%! v1 = betainc (1,a,b);
+%! v2 = [1,1,1,1];
+%! x = [.2, .4, .6, .8];
+%! v3 = betainc (x, a, b);
+%! v4 = 1 - betainc (1.-x, b, a);
+%! assert (v1, v2, sqrt (eps));
+%! assert (v3, v4, sqrt (eps));
+
+# Single precision
+%!test
+%! a = single ([1, 1.5, 2, 3]);
+%! b = single ([4, 3, 2, 1]);
+%! v1 = betainc (1,a,b);
+%! v2 = single ([1,1,1,1]);
+%! x = single ([.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")));
+
+# Mixed double/single precision
+%!test
+%! a = single ([1, 1.5, 2, 3]);
+%! b = [4, 3, 2, 1];
+%! v1 = betainc (1,a,b);
+%! v2 = single ([1,1,1,1]);
+%! 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")));
+
+## New test
+
+%!test #<51157>
+%! y = betainc([0.00780;0.00782;0.00784],250.005,49750.995);
+%! y_ex = [0.999999999999989; 0.999999999999992; 0.999999999999995];
+%! assert (y, y_ex, -1e-14);
+
+%!assert (betainc (0.001, 20, 30), 2.750687665855991e-47, -3e-14);
+%!assert (betainc (0.0001, 20, 30), 2.819953178893307e-67, -3e-14);
+%!assert (betainc (0.99, 20, 30, "upper"), 1.5671643161872703e-47, -3e-14);
+%!assert (betainc (0.999, 20, 30, "upper"), 1.850806276141535e-77, -3e-14);
+%!assert (betainc (0.5, 200, 300), 0.9999964565197356, -1e-15);
+%!assert (betainc (0.5, 200, 300, "upper"), 3.54348026439253e-06, -1e-13);
+
+# Test trivial values
+
+%!test
+%! [a,b] = ndgrid (linspace(1e-4, 100, 20), linspace(1e-4, 100, 20));
+%! assert (betainc (0,a,b), zeros(20));
+%! assert (betainc (1,a,b), ones(20));
+
+## Test input validation
+
+%!error betainc ()
+%!error betainc (1)
+%!error betainc (1,2)
+%!error betainc (1.1,1,1)
+%!error betainc (-0.1,1,1)
+%!error betainc (0.5,0,1)
+%!error betainc (0.5,1,0)
+%!error betainc (1,2,3,4)
+%!error betainc (1,2)
+%!error betainc (1,2,3,4,5)
+%!error betainc (0.5i, 1, 2)
+%!error betainc (0, 1i, 1)
+%!error betainc (0, 1, 1i)
diff -r 451f4f291f46 -r bd89440407aa scripts/specfun/betaincinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/specfun/betaincinv.m	Thu Sep 07 16:13:16 2017 +0200
@@ -0,0 +1,283 @@
+## Copyright (C) 2017 Michele Ginesi
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## Author: Michele Ginesi <michele.ginesi@gmail.com>
+
+## -*- texinfo -*-
+## @deftypefn  {} {} betaincinv (@var{y}, @var{a}, @var{b})
+## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "lower")
+## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "upper")
+## Compute the inverse of the normalized incomplete beta function.
+##
+## The normalized incomplete beta function is defined as
+## @tex
+## $$
+##  I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt}
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##                x
+##               /
+##              |
+## I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt
+##              |
+##              /
+##             0
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+## If two inputs are scalar, then @code{betaincinv (@var{y}, @var{a}, @var{b})}
+## is returned for each of the other inputs.
+##
+## If two or more inputs are not scalar, the sizes of them must agree, and
+## @code{betaincinv} is applied element-by-element.
+##
+## The variable @var{y} must be in the interval [0,1], while @var{a} and @var{b}
+## must be real and strictly positive.
+##
+## By default the inverse of the incomplete beta function integrated from 0
+## to @var{x} is computed.  If @qcode{"upper"} is given then the complementary
+## function integrated from @var{x} to 1 is inverted.
+##
+## The function is computed by standard Newton's method, by solving
+## @tex
+## $$
+##  y - I_x (a, b) = 0
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##
+## @var{y} - betainc (@var{x}, @var{a}, @var{b}) = 0
+##
+## @end group
+## @end example
+##
+## @end ifnottex
+##
+##
+## @seealso{beta, betainc, betaln}
+## @end deftypefn
+
+function [x] = betaincinv (y, a, b, tail = "lower")
+
+  if (nargin > 4 || nargin < 3)
+    print_usage ();
+  endif
+
+  if (! isscalar (y) || ! isscalar (a) || ! isscalar (b))
+    [err, y, a, b] = common_size (y, a, b);
+    if (err > 0)
+      error ("betaincinv: y, a must be of common size or scalars");
+    endif
+  endif
+
+  if (iscomplex (y) || iscomplex (a) || iscomplex (b))
+    error ("betaincinv: inputs must be real or integer");
+  endif
+
+  if (any (a <= 0) | any (b <= 0))
+    error ("betaincinv: a must be strictly positive");
+  endif
+
+  if (any (y > 1 | y < 0))
+    error ("betaincinv: y must be between 0 and 1");
+  endif
+
+  if (isinteger (y))
+    y = double (y);
+  endif
+
+  if (isinteger (a))
+    a = double (a);
+  endif
+
+  if (isinteger (b))
+    b = double (b);
+  endif
+
+
+  # Extract the size.
+  sz = size (y);
+  # Write the inputs as two column vectors.
+  y = y(:);
+  a = a(:);
+  b = b(:);
+  l = length (y);
+  # Initialise the output.
+  x = zeros (l, 1);
+
+  # If one of the inputs is of single type, also the output should be
+
+  if (strcmpi (class (y), "single") || strcmpi (class (a), "single") || strcmpi (class (b), "single"))
+    a = single (a);
+    b = single (b);
+    y = single (y);
+    x = single (x);
+  endif
+
+  # Parameters of the Newton method
+  maxit = 20;
+  tol = eps (class (y));
+
+  if (strcmpi (tail, "upper"))
+    p = 1 - y;
+    q = y;
+    x(y == 0) = 1;
+    x(y == 1) = 0;
+  elseif (strcmpi (tail, "lower"))
+    p = y;
+    q = 1 - y;
+    x(y == 0) = 0;
+    x(y == 1) = 1;
+  else
+    error ("betaincinv: invalid value for tail")
+  endif
+
+  # Trivial values have been already computed.
+  i_miss = ((y != 0) & (y != 1));
+
+  # We will invert the lower version for p < 0.5 and the upper otherwise.
+  i_low = (p < 0.5);
+  i_upp = (!i_low);
+
+  len_low = nnz (i_miss & i_low);
+  len_upp = nnz (i_miss & i_upp);
+
+  if (any (i_miss & i_low));
+    # Function and derivative of the lower version.
+    F = @(x, a, b, y) y - betainc (x, a, b);
+    JF = @(x, a, b) - real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ...
+        gammaln (a+b) - gammaln (a) - gammaln (b)));
+    # We compute the initial guess with a bisection method of 10 steps.
+    x0 = bisection_method (F, zeros (len_low,1), ones (len_low,1),...
+        a(i_low), b(i_low), p(i_low), 10);
+    x(i_low) = newton_method (F, JF, x0, a(i_low), b(i_low), p(i_low), ...
+        tol, maxit);
+  endif
+
+  if (any (i_miss & i_upp));
+    # Function and derivative of the lower version.
+    F = @(x, a, b, y) y - betainc (x, a, b, "upper");
+    JF = @(x, a, b) real (exp ((a-1) .* log (x) + (b-1) .* log1p (-x) + ...
+        gammaln (a+b) - gammaln (a) - gammaln (b)));
+    # We compute the initial guess with a bisection method of 10 steps.
+    x0 = bisection_method (F, zeros (len_upp,1), ones (len_upp,1),...
+        a(i_upp), b(i_upp), q(i_upp), 10);
+    x(i_upp) = newton_method (F, JF, x0, a(i_upp), b(i_upp), q(i_upp), ...
+        tol, maxit);
+  endif
+
+  ## Re-organize the output.
+
+  x = reshape (x, sz);
+
+endfunction
+
+
+## Subfunctions: Bisection and Newton Methods
+
+function xc = bisection_method (F, xl, xr, a, b, y, maxit)
+    F_l = feval (F, xl, a, b, y);
+    F_r = feval (F, xr, a, b, y);
+   for it = 1:maxit
+    xc = (xl + xr) / 2;
+    F_c = feval (F, xc, a, b, y);
+    flag_l = ((F_c .* F_r) < 0);
+    flag_r = ((F_c .* F_l) < 0);
+    flag_c = (F_c == 0);
+    xl(flag_l) = xc(flag_l);
+    xr(flag_r) = xc(flag_r);
+    xl(flag_c) = xr(flag_c) = xc(flag_c);
+    F_l(flag_l) = F_c(flag_l);
+    F_r(flag_r) = F_c(flag_r);
+    F_l(flag_c) = F_r(flag_c) = 0;
+   endfor
+endfunction
+
+function x = newton_method (F, JF, x0, a, b, y, tol, maxit);
+  l = length (y);
+  res = -feval (F, x0, a, b, y) ./ feval (JF, x0, a, b);
+  i_miss = (abs(res) >= tol * abs (x0));
+  x = x0;
+  it = 0;
+  while (any (i_miss) && (it < maxit))
+    it++;
+    x(i_miss) += res(i_miss);
+    res(i_miss) = - feval (F, x(i_miss), a(i_miss), b(i_miss), ...
+        y(i_miss)) ./ feval (JF, x(i_miss), a(i_miss), b(i_miss));
+    i_miss = (abs(res) >= tol * abs (x));
+  endwhile
+  x += res;
+endfunction
+
+
+## Test
+
+%!test
+%! x = linspace(0.1, 0.9, 11);
+%! a = [2, 3, 4];
+%! [x,a,b] = ndgrid (x,a,a);
+%! xx = betaincinv (betainc (x, a, b), a, b);
+%! assert (xx, x, 3e-15);
+
+%!test
+%! x = linspace(0.1, 0.9, 11);
+%! a = [2, 3, 4];
+%! [x,a,b] = ndgrid (x,a,a);
+%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper");
+%! assert (xx, x, 3e-15);
+
+%!test
+%! x = linspace(0.1, 0.9, 11);
+%! a = [0.1:0.1:1];
+%! [x,a,b] = ndgrid (x,a,a);
+%! xx = betaincinv (betainc (x, a, b), a, b);
+%! assert (xx, x, 3e-15);
+
+%!test
+%! x = linspace(0.1, 0.9, 11);
+%! a = [0.1:0.1:1];
+%! [x,a,b] = ndgrid (x,a,a);
+%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper");
+%! assert (xx, x, 3e-15);
+
+## Test on the conservation of the class
+%!assert (class (betaincinv (0.5, 1, 1)), "double")
+%!assert (class (betaincinv (single (0.5), 1, 1)), "single")
+%!assert (class (betaincinv (0.5, single (1), 1)), "single")
+%!assert (class (betaincinv (int8 (0), 1, 1)), "double")
+%!assert (class (betaincinv (0.5, int8 (1), 1)), "double")
+%!assert (class (betaincinv (int8 (0), single (1), 1)), "single")
+%!assert (class (betaincinv (single (0.5), int8 (1), 1)), "single")
+
+## Test input validation
+%!error betaincinv ()
+%!error betaincinv (1)
+%!error betaincinv (1,2,3,4)
+%!error betaincinv (1, "2")
+%!error betaincinv (0.5i, 1, 2)
+%!error betaincinv (0, 1i, 1)
+%!error betaincinv (0, 1, 1i)
diff -r 451f4f291f46 -r bd89440407aa scripts/specfun/module.mk
--- a/scripts/specfun/module.mk	Thu Sep 07 15:58:26 2017 +0200
+++ b/scripts/specfun/module.mk	Thu Sep 07 16:13:16 2017 +0200
@@ -2,6 +2,8 @@
 
 %canon_reldir%_FCN_FILES = \
   %reldir%/beta.m \
+  %reldir%/betainc.m \
+  %reldir%/betaincinv.m \
   %reldir%/betaln.m \
   %reldir%/ellipke.m \
   %reldir%/expint.m \