# HG changeset patch
# User Michele Ginesi <michele.ginesi@gmail.com>
# Date 1486251018 -3600
#      Sun Feb 05 00:30:18 2017 +0100
# Node ID 4ab8551eb58e398795a8774b2c53009effef725f
# Parent  d05f688ae836dceddb765e670999849830129d0e
New strategy to compute gammainc (see bug #47800).

* scripts/specfun/gammainc.m: new file using seven differents algorithms
depending on the values of a and x
* scripts/specfun/module.mk: added gammainc.m
* libinterp/corefcn/__gammainc_lentz__.cc: new file using continued fractions
* libinterp/corefcn/gammainc.cc: removed
* liboctave/cruft/slatec-fn/dgami.f,dgamit.f,gami.f,gamit.f,xdgami.f,
xdgamit.f,xgmainc.f,xsgmainc.f: removed
* libinterp/corefcn/module.mk: removed gammainc.cc, added __gammainc_lentz__.cc
* liboctave/numeric/lo-specfun.cc: removed all lines regarding gammainc

diff -r d05f688ae836 -r 4ab8551eb58e libinterp/corefcn/__gammainc_lentz__.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__gammainc_lentz__.cc	Sun Feb 05 00:30:18 2017 +0100
@@ -0,0 +1,89 @@
+// Copyright (C) 2017 Nir Krakauer
+//
+// 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 (__gammainc_lentz__, args, , "Continued fraction for incomplete gamma function")
+{
+  int nargin = args.length ();
+
+  if (nargin != 2)
+    print_usage ();
+  else
+    {
+      octave_value x_arg = args(0);
+      octave_value a_arg = args(1);
+      if (x_arg.is_single_type () || a_arg.is_single_type ())
+        {
+          const float x = args(0).float_value ();
+          const float a = args(1).float_value ();
+          static const float tiny = pow (2, -50);
+          static const float eps = std::numeric_limits<float>::epsilon();
+          float y = tiny;
+          float Cj = y;
+          float Dj = 0;
+          float bj = x - a + 1;
+          float aj = a;
+          float Deltaj = 0;
+          int j = 1;
+          while(std::abs((Deltaj - 1) / y)  > eps)
+            {
+               Cj = bj + aj/Cj;
+               Dj = 1 / (bj + aj*Dj);
+               Deltaj = Cj * Dj;
+               y *= Deltaj;
+               bj += 2;
+               aj = j * (a - j);
+               j++;
+             }
+           if (! error_state)
+            return octave_value (y);
+        }
+      else
+        {
+          const double x = args(0).double_value ();
+          const double a = args(1).double_value ();
+          static const double tiny = pow (2, -100);
+          static const double eps = std::numeric_limits<double>::epsilon();
+          double y = tiny;
+          double Cj = y;
+          double Dj = 0;
+          double bj = x - a + 1;
+          double aj = a;
+          double Deltaj = 0;
+          int j = 1;
+          while(std::abs((Deltaj - 1) / y)  > eps)
+            {
+              Cj = bj + aj/Cj;
+              Dj = 1 / (bj + aj*Dj);
+              Deltaj = Cj * Dj;
+              y *= Deltaj;
+              bj += 2;
+              aj = j * (a - j);
+              j++;
+            }
+          if (! error_state)
+            return octave_value (y);
+        }
+      }
+  return octave_value_list ();
+}
diff -r d05f688ae836 -r 4ab8551eb58e libinterp/corefcn/gammainc.cc
--- a/libinterp/corefcn/gammainc.cc	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,213 +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
-<http://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 (gammainc, args, ,
-       doc: /* -*- texinfo -*-
-@deftypefn  {} {} gammainc (@var{x}, @var{a})
-@deftypefnx {} {} gammainc (@var{x}, @var{a}, "lower")
-@deftypefnx {} {} gammainc (@var{x}, @var{a}, "upper")
-Compute the normalized incomplete gamma function.
-
-This is defined as
-@tex
-$$
- \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt}
-$$
-@end tex
-@ifnottex
-
-@example
-@group
-                                x
-                       1       /
-gammainc (x, a) = ---------    | exp (-t) t^(a-1) dt
-                  gamma (a)    /
-                            t=0
-@end group
-@end example
-
-@end ifnottex
-with the limiting value of 1 as @var{x} approaches infinity.
-The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and
-@nospell{Stegun} (6.5.1).
-
-If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned
-for each element of @var{x} and vice versa.
-
-If neither @var{x} nor @var{a} is scalar, the sizes of @var{x} and
-@var{a} must agree, and @code{gammainc} is applied element-by-element.
-
-By default the incomplete gamma function integrated from 0 to @var{x} is
-computed.  If @qcode{"upper"} is given then the complementary function
-integrated from @var{x} to infinity is calculated.  It should be noted that
-
-@example
-gammainc (@var{x}, @var{a}) @equiv{} 1 - gammainc (@var{x}, @var{a}, "upper")
-@end example
-@seealso{gamma, gammaln}
-@end deftypefn */)
-{
-  int nargin = args.length ();
-
-  if (nargin < 2 || nargin > 3)
-    print_usage ();
-
-  bool lower = true;
-  if (nargin == 3)
-    {
-      std::string s = args(2).xstring_value ("gammainc: third argument must be \"lower\" or \"upper\"");
-
-      std::transform (s.begin (), s.end (), s.begin (), tolower);
-
-      if (s == "upper")
-        lower = false;
-      else if (s == "lower")
-        lower = true;
-      else
-        error ("gammainc: third argument must be \"lower\" or \"upper\"");
-    }
-
-  octave_value retval;
-
-  octave_value x_arg = args(0);
-  octave_value a_arg = args(1);
-
-  // FIXME: Can we make a template version of the duplicated code below
-  if (x_arg.is_single_type () || a_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 ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0f - octave::math::gammainc (x, a);
-            }
-          else
-            {
-              FloatNDArray a = a_arg.float_array_value ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0f - octave::math::gammainc (x, a);
-            }
-        }
-      else
-        {
-          FloatNDArray x = x_arg.float_array_value ();
-
-          if (a_arg.is_scalar_type ())
-            {
-              float a = a_arg.float_value ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0f - octave::math::gammainc (x, a);
-            }
-          else
-            {
-              FloatNDArray a = a_arg.float_array_value ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0f - octave::math::gammainc (x, a);
-            }
-        }
-    }
-  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 ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0 - octave::math::gammainc (x, a);
-            }
-          else
-            {
-              NDArray a = a_arg.array_value ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0 - octave::math::gammainc (x, a);
-            }
-        }
-      else
-        {
-          NDArray x = x_arg.array_value ();
-
-          if (a_arg.is_scalar_type ())
-            {
-              double a = a_arg.double_value ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0 - octave::math::gammainc (x, a);
-            }
-          else
-            {
-              NDArray a = a_arg.array_value ();
-
-              retval = lower ? octave::math::gammainc (x, a)
-                             : 1.0 - octave::math::gammainc (x, a);
-            }
-        }
-    }
-
-  return retval;
-}
-
-/*
-%!test
-%! a = [.5 .5 .5 .5 .5];
-%! x = [0 1 2 3 4];
-%! v1 = sqrt (pi)*erf (x)./gamma (a);
-%! v3 = gammainc (x.*x, a);
-%! assert (v1, v3, sqrt (eps));
-
-%!assert (gammainc (0:4,0.5, "upper"), 1-gammainc (0:4,0.5), 1e-10)
-
-%!test
-%! a = single ([.5 .5 .5 .5 .5]);
-%! x = single ([0 1 2 3 4]);
-%! v1 = sqrt (pi ("single"))*erf (x)./gamma (a);
-%! v3 = gammainc (x.*x, a);
-%! assert (v1, v3, sqrt (eps ("single")));
-
-%!assert (gammainc (single (0:4), single (0.5), "upper"),
-%!        single (1)-gammainc (single (0:4), single (0.5)),
-%!        single (1e-7))
-*/
diff -r d05f688ae836 -r 4ab8551eb58e libinterp/corefcn/module.mk
--- a/libinterp/corefcn/module.mk	Mon Mar 20 15:32:08 2017 -0700
+++ b/libinterp/corefcn/module.mk	Sun Feb 05 00:30:18 2017 +0100
@@ -103,6 +103,7 @@ COREFCN_SRC = \
   libinterp/corefcn/Cell.cc \
   libinterp/corefcn/__contourc__.cc \
   libinterp/corefcn/__dsearchn__.cc \
+  libinterp/corefcn/__gammainc_lentz__.cc \
   libinterp/corefcn/__ichol__.cc \
   libinterp/corefcn/__ilu__.cc \
   libinterp/corefcn/__lin_interpn__.cc \
@@ -148,7 +149,6 @@ COREFCN_SRC = \
   libinterp/corefcn/filter.cc \
   libinterp/corefcn/find.cc \
   libinterp/corefcn/ft-text-renderer.cc \
-  libinterp/corefcn/gammainc.cc \
   libinterp/corefcn/gcd.cc \
   libinterp/corefcn/getgrent.cc \
   libinterp/corefcn/getpwent.cc \
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/dgami.f
--- a/liboctave/cruft/slatec-fn/dgami.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-
-*DECK DGAMI
-      DOUBLE PRECISION FUNCTION DGAMI (A, X)
-C***BEGIN PROLOGUE  DGAMI
-C***PURPOSE  Evaluate the incomplete Gamma function.
-C***LIBRARY   SLATEC (FNLIB)
-C***CATEGORY  C7E
-C***TYPE      DOUBLE PRECISION (GAMI-S, DGAMI-D)
-C***KEYWORDS  FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS
-C***AUTHOR  Fullerton, W., (LANL)
-C***DESCRIPTION
-C
-C Evaluate the incomplete gamma function defined by
-C
-C DGAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) .
-C
-C DGAMI is evaluated for positive values of A and non-negative values
-C of X.  A slight deterioration of 2 or 3 digits accuracy will occur
-C when DGAMI is very large or very small, because logarithmic variables
-C are used.  The function and both arguments are double precision.
-C
-C***REFERENCES  (NONE)
-C***ROUTINES CALLED  DGAMIT, DLNGAM, XERMSG
-C***REVISION HISTORY  (YYMMDD)
-C   770701  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***END PROLOGUE  DGAMI
-      DOUBLE PRECISION A, X, FACTOR, DLNGAM, DGAMIT
-C***FIRST EXECUTABLE STATEMENT  DGAMI
-      IF (A .LE. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI',
-     +   'A MUST BE GT ZERO', 1, 2)
-      IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI',
-     +   'X MUST BE GE ZERO', 2, 2)
-C
-      DGAMI = 0.D0
-      IF (X.EQ.0.0D0) RETURN
-C
-C THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW.
-      FACTOR = EXP (DLNGAM(A) + A*LOG(X))
-C
-      DGAMI = FACTOR * DGAMIT (A, X)
-C
-      RETURN
-      END
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/dgamit.f
--- a/liboctave/cruft/slatec-fn/dgamit.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,119 +0,0 @@
-*DECK DGAMIT
-      DOUBLE PRECISION FUNCTION DGAMIT (A, X)
-C***BEGIN PROLOGUE  DGAMIT
-C***PURPOSE  Calculate Tricomi's form of the incomplete Gamma function.
-C***LIBRARY   SLATEC (FNLIB)
-C***CATEGORY  C7E
-C***TYPE      DOUBLE PRECISION (GAMIT-S, DGAMIT-D)
-C***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
-C             SPECIAL FUNCTIONS, TRICOMI
-C***AUTHOR  Fullerton, W., (LANL)
-C***DESCRIPTION
-C
-C   Evaluate Tricomi's incomplete Gamma function defined by
-C
-C   DGAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
-C              T**(A-1.)
-C
-C   for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
-C   GAMMA(X) is the complete gamma function of X.
-C
-C   DGAMIT is evaluated for arbitrary real values of A and for non-
-C   negative values of X (even though DGAMIT is defined for X .LT.
-C   0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite,
-C   which is a fatal error.
-C
-C   The function and both arguments are DOUBLE PRECISION.
-C
-C   A slight deterioration of 2 or 3 digits accuracy will occur when
-C   DGAMIT is very large or very small in absolute value, because log-
-C   arithmic variables are used.  Also, if the parameter  A  is very
-C   close to a negative integer (but not a negative integer), there is
-C   a loss of accuracy, which is reported if the result is less than
-C   half machine precision.
-C
-C***REFERENCES  W. Gautschi, A computational procedure for incomplete
-C                 gamma functions, ACM Transactions on Mathematical
-C                 Software 5, 4 (December 1979), pp. 466-481.
-C               W. Gautschi, Incomplete gamma functions, Algorithm 542,
-C                 ACM Transactions on Mathematical Software 5, 4
-C                 (December 1979), pp. 482-489.
-C***ROUTINES CALLED  D1MACH, D9GMIT, D9LGIC, D9LGIT, DGAMR, DLGAMS,
-C                    DLNGAM, XERCLR, XERMSG
-C***REVISION HISTORY  (YYMMDD)
-C   770701  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   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
-C***END PROLOGUE  DGAMIT
-      DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
-     1  BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, D9LGIT,
-     2  DLNGAM, D9LGIC
-      LOGICAL FIRST
-      SAVE ALNEPS, SQEPS, BOT, FIRST
-      DATA FIRST /.TRUE./
-C***FIRST EXECUTABLE STATEMENT  DGAMIT
-      IF (FIRST) THEN
-         ALNEPS = -LOG (D1MACH(3))
-         SQEPS = SQRT(D1MACH(4))
-         BOT = LOG (D1MACH(1))
-      ENDIF
-      FIRST = .FALSE.
-C
-      IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMIT', 'X IS NEGATIVE'
-     +   , 2, 2)
-C
-      IF (X.NE.0.D0) ALX = LOG (X)
-      SGA = 1.0D0
-      IF (A.NE.0.D0) SGA = SIGN (1.0D0, A)
-      AINTA = AINT (A + 0.5D0*SGA)
-      AEPS = A - AINTA
-C
-      IF (X.GT.0.D0) GO TO 20
-      DGAMIT = 0.0D0
-      IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
-      RETURN
-C
- 20   IF (X.GT.1.D0) GO TO 30
-      IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1,
-     1  SGNGAM)
-      DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
-      RETURN
-C
- 30   IF (A.LT.X) GO TO 40
-      T = D9LGIT (A, X, DLNGAM(A+1.0D0))
-      IF (T.LT.BOT) CALL XERCLR
-      DGAMIT = EXP (T)
-      RETURN
-C
- 40   ALNG = D9LGIC (A, X, ALX)
-C
-C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
-C
-      H = 1.0D0
-      IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50
-C
-      CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM)
-      T = LOG (ABS(A)) + ALNG - ALGAP1
-      IF (T.GT.ALNEPS) GO TO 60
-C
-      IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T)
-      IF (ABS(H).GT.SQEPS) GO TO 50
-C
-      CALL XERCLR
-      CALL XERMSG ('SLATEC', 'DGAMIT', 'RESULT LT HALF PRECISION', 1,
-     +   1)
-C
- 50   T = -A*ALX + LOG(ABS(H))
-      IF (T.LT.BOT) CALL XERCLR
-      DGAMIT = SIGN (EXP(T), H)
-      RETURN
-C
- 60   T = T - A*ALX
-      IF (T.LT.BOT) CALL XERCLR
-      DGAMIT = -SGA * SGNGAM * EXP(T)
-      RETURN
-C
-      END
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/gami.f
--- a/liboctave/cruft/slatec-fn/gami.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,45 +0,0 @@
-*DECK GAMI
-      FUNCTION GAMI (A, X)
-C***BEGIN PROLOGUE  GAMI
-C***PURPOSE  Evaluate the incomplete Gamma function.
-C***LIBRARY   SLATEC (FNLIB)
-C***CATEGORY  C7E
-C***TYPE      SINGLE PRECISION (GAMI-S, DGAMI-D)
-C***KEYWORDS  FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS
-C***AUTHOR  Fullerton, W., (LANL)
-C***DESCRIPTION
-C
-C Evaluate the incomplete gamma function defined by
-C
-C GAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) .
-C
-C GAMI is evaluated for positive values of A and non-negative values
-C of X.  A slight deterioration of 2 or 3 digits accuracy will occur
-C when GAMI is very large or very small, because logarithmic variables
-C are used.  GAMI, A, and X are single precision.
-C
-C***REFERENCES  (NONE)
-C***ROUTINES CALLED  ALNGAM, GAMIT, XERMSG
-C***REVISION HISTORY  (YYMMDD)
-C   770701  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***END PROLOGUE  GAMI
-C***FIRST EXECUTABLE STATEMENT  GAMI
-      IF (A .LE. 0.0) CALL XERMSG ('SLATEC', 'GAMI',
-     +   'A MUST BE GT ZERO', 1, 2)
-      IF (X .LT. 0.0) CALL XERMSG ('SLATEC', 'GAMI',
-     +   'X MUST BE GE ZERO', 2, 2)
-C
-      GAMI = 0.0
-      IF (X.EQ.0.0) RETURN
-C
-C THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW.
-      FACTOR = EXP (ALNGAM(A) + A*LOG(X) )
-C
-      GAMI = FACTOR * GAMIT(A, X)
-C
-      RETURN
-      END
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/gamit.f
--- a/liboctave/cruft/slatec-fn/gamit.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,112 +0,0 @@
-*DECK GAMIT
-      REAL FUNCTION GAMIT (A, X)
-C***BEGIN PROLOGUE  GAMIT
-C***PURPOSE  Calculate Tricomi's form of the incomplete Gamma function.
-C***LIBRARY   SLATEC (FNLIB)
-C***CATEGORY  C7E
-C***TYPE      SINGLE PRECISION (GAMIT-S, DGAMIT-D)
-C***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB,
-C             SPECIAL FUNCTIONS, TRICOMI
-C***AUTHOR  Fullerton, W., (LANL)
-C***DESCRIPTION
-C
-C   Evaluate Tricomi's incomplete gamma function defined by
-C
-C   GAMIT = X**(-A)/GAMMA(A) * integral from 0 to X of EXP(-T) *
-C             T**(A-1.)
-C
-C   for A .GT. 0.0 and by analytic continuation for A .LE. 0.0.
-C   GAMMA(X) is the complete gamma function of X.
-C
-C   GAMIT is evaluated for arbitrary real values of A and for non-
-C   negative values of X (even though GAMIT is defined for X .LT.
-C   0.0), except that for X = 0 and A .LE. 0.0, GAMIT is infinite,
-C   which is a fatal error.
-C
-C   The function and both arguments are REAL.
-C
-C   A slight deterioration of 2 or 3 digits accuracy will occur when
-C   GAMIT is very large or very small in absolute value, because log-
-C   arithmic variables are used.  Also, if the parameter  A  is very
-C   close to a negative integer (but not a negative integer), there is
-C   a loss of accuracy, which is reported if the result is less than
-C   half machine precision.
-C
-C***REFERENCES  W. Gautschi, A computational procedure for incomplete
-C                 gamma functions, ACM Transactions on Mathematical
-C                 Software 5, 4 (December 1979), pp. 466-481.
-C               W. Gautschi, Incomplete gamma functions, Algorithm 542,
-C                 ACM Transactions on Mathematical Software 5, 4
-C                 (December 1979), pp. 482-489.
-C***ROUTINES CALLED  ALGAMS, ALNGAM, GAMR, R1MACH, R9GMIT, R9LGIC,
-C                    R9LGIT, XERCLR, XERMSG
-C***REVISION HISTORY  (YYMMDD)
-C   770701  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   920528  DESCRIPTION and REFERENCES sections revised.  (WRB)
-C***END PROLOGUE  GAMIT
-      LOGICAL FIRST
-      SAVE ALNEPS, SQEPS, BOT, FIRST
-      DATA FIRST /.TRUE./
-C***FIRST EXECUTABLE STATEMENT  GAMIT
-      IF (FIRST) THEN
-         ALNEPS = -LOG(R1MACH(3))
-         SQEPS = SQRT(R1MACH(4))
-         BOT = LOG(R1MACH(1))
-      ENDIF
-      FIRST = .FALSE.
-C
-      IF (X .LT. 0.0) CALL XERMSG ('SLATEC', 'GAMIT', 'X IS NEGATIVE',
-     +   2, 2)
-C
-      IF (X.NE.0.0) ALX = LOG(X)
-      SGA = 1.0
-      IF (A.NE.0.0) SGA = SIGN (1.0, A)
-      AINTA = AINT (A+0.5*SGA)
-      AEPS = A - AINTA
-C
-      IF (X.GT.0.0) GO TO 20
-      GAMIT = 0.0
-      IF (AINTA.GT.0.0 .OR. AEPS.NE.0.0) GAMIT = GAMR(A+1.0)
-      RETURN
-C
- 20   IF (X.GT.1.0) GO TO 40
-      IF (A.GE.(-0.5) .OR. AEPS.NE.0.0) CALL ALGAMS (A+1.0, ALGAP1,
-     1  SGNGAM)
-      GAMIT = R9GMIT (A, X, ALGAP1, SGNGAM, ALX)
-      RETURN
-C
- 40   IF (A.LT.X) GO TO 50
-      T = R9LGIT (A, X, ALNGAM(A+1.0))
-      IF (T.LT.BOT) CALL XERCLR
-      GAMIT = EXP(T)
-      RETURN
-C
- 50   ALNG = R9LGIC (A, X, ALX)
-C
-C EVALUATE GAMIT IN TERMS OF LOG(GAMIC(A,X))
-C
-      H = 1.0
-      IF (AEPS.EQ.0.0 .AND. AINTA.LE.0.0) GO TO 60
-      CALL ALGAMS (A+1.0, ALGAP1, SGNGAM)
-      T = LOG(ABS(A)) + ALNG - ALGAP1
-      IF (T.GT.ALNEPS) GO TO 70
-      IF (T.GT.(-ALNEPS)) H = 1.0 - SGA*SGNGAM*EXP(T)
-      IF (ABS(H).GT.SQEPS) GO TO 60
-      CALL XERCLR
-      CALL XERMSG ('SLATEC', 'GAMIT', 'RESULT LT HALF PRECISION', 1, 1)
-C
- 60   T = -A*ALX + LOG(ABS(H))
-      IF (T.LT.BOT) CALL XERCLR
-      GAMIT = SIGN (EXP(T), H)
-      RETURN
-C
- 70   T = T - A*ALX
-      IF (T.LT.BOT) CALL XERCLR
-      GAMIT = -SGA*SGNGAM*EXP(T)
-      RETURN
-C
-      END
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/module.mk
--- a/liboctave/cruft/slatec-fn/module.mk	Mon Mar 20 15:32:08 2017 -0700
+++ b/liboctave/cruft/slatec-fn/module.mk	Sun Feb 05 00:30:18 2017 +0100
@@ -18,8 +18,6 @@ CRUFT_SOURCES += \
   liboctave/cruft/slatec-fn/dbetai.f \
   liboctave/cruft/slatec-fn/dcsevl.f \
   liboctave/cruft/slatec-fn/derf.f \
-  liboctave/cruft/slatec-fn/dgami.f \
-  liboctave/cruft/slatec-fn/dgamit.f \
   liboctave/cruft/slatec-fn/dgamlm.f \
   liboctave/cruft/slatec-fn/dgamma.f \
   liboctave/cruft/slatec-fn/dgamr.f \
@@ -31,8 +29,6 @@ CRUFT_SOURCES += \
   liboctave/cruft/slatec-fn/dpchst.f \
   liboctave/cruft/slatec-fn/dpsifn.f \
   liboctave/cruft/slatec-fn/erf.f \
-  liboctave/cruft/slatec-fn/gami.f \
-  liboctave/cruft/slatec-fn/gamit.f \
   liboctave/cruft/slatec-fn/gamlim.f \
   liboctave/cruft/slatec-fn/gamma.f \
   liboctave/cruft/slatec-fn/gamr.f \
@@ -51,16 +47,12 @@ CRUFT_SOURCES += \
   liboctave/cruft/slatec-fn/xdbetai.f \
   liboctave/cruft/slatec-fn/xderf.f \
   liboctave/cruft/slatec-fn/xderfc.f \
-  liboctave/cruft/slatec-fn/xdgami.f \
-  liboctave/cruft/slatec-fn/xdgamit.f \
   liboctave/cruft/slatec-fn/xdgamma.f \
-  liboctave/cruft/slatec-fn/xgmainc.f \
   liboctave/cruft/slatec-fn/xacosh.f \
   liboctave/cruft/slatec-fn/xasinh.f \
   liboctave/cruft/slatec-fn/xatanh.f \
   liboctave/cruft/slatec-fn/xerf.f \
   liboctave/cruft/slatec-fn/xerfc.f \
-  liboctave/cruft/slatec-fn/xsgmainc.f \
   liboctave/cruft/slatec-fn/xgamma.f \
   liboctave/cruft/slatec-fn/xbetai.f
 
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/xdgami.f
--- a/liboctave/cruft/slatec-fn/xdgami.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-      subroutine xdgami (a, x, result)
-      external dgami
-      double precision a, x, result, dgami
-      result = dgami (a, x)
-      return
-      end
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/xdgamit.f
--- a/liboctave/cruft/slatec-fn/xdgamit.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-      subroutine xdgamit (a, x, result)
-      external dgamit
-      double precision a, x, result, dgamit
-      result = dgamit (a, x)
-      return
-      end
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/xgmainc.f
--- a/liboctave/cruft/slatec-fn/xgmainc.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,100 +0,0 @@
-      subroutine xgammainc (a, x, result)
-
-c -- jwe, based on DGAMIT.
-c
-c -- Do a better job than dgami for large values of x.
-
-      double precision a, x, result
-      intrinsic exp, log, sqrt, sign, aint
-      external dgami, dlngam, d9lgit, d9lgic, d9gmit
-
-C     external dgamr
-C     DOUBLE PRECISION DGAMR
-
-      DOUBLE PRECISION AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
-     $     BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, D9GMIT,
-     $     D9LGIC, D9LGIT, DLNGAM, DGAMI
-
-      LOGICAL FIRST
-
-      SAVE ALNEPS, SQEPS, BOT, FIRST
-
-      DATA FIRST /.TRUE./
-
-      if (x .eq. 0.0d0) then
-
-        if (a .eq. 0.0d0) then
-          result = 1.0d0
-        else
-          result = 0.0d0
-        endif
-
-      else
-
-      IF (FIRST) THEN
-         ALNEPS = -LOG (D1MACH(3))
-         SQEPS = SQRT(D1MACH(4))
-         BOT = LOG (D1MACH(1))
-      ENDIF
-      FIRST = .FALSE.
-C
-      IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'XGMAINC', 'X IS NEGATIVE'
-     +   , 2, 2)
-C
-      IF (X.NE.0.D0) ALX = LOG (X)
-      SGA = 1.0D0
-      IF (A.NE.0.D0) SGA = SIGN (1.0D0, A)
-      AINTA = AINT (A + 0.5D0*SGA)
-      AEPS = A - AINTA
-C
-C      IF (X.GT.0.D0) GO TO 20
-C      DGAMIT = 0.0D0
-C      IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0)
-C      RETURN
-C
- 20   IF (X.GT.1.D0) GO TO 30
-      IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1,
-     1  SGNGAM)
-C      DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX)
-      result = exp (a*alx + log (D9GMIT (A, X, ALGAP1, SGNGAM, ALX)))
-      RETURN
-C
- 30   IF (A.LT.X) GO TO 40
-      T = D9LGIT (A, X, DLNGAM(A+1.0D0))
-      IF (T.LT.BOT) CALL XERCLR
-C      DGAMIT = EXP (T)
-      result = EXP (a*alx + T)
-      RETURN
-C
- 40   ALNG = D9LGIC (A, X, ALX)
-C
-C EVALUATE DGAMIT IN TERMS OF LOG (DGAMIC (A, X))
-C
-      H = 1.0D0
-      IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50
-C
-      CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM)
-      T = LOG (ABS(A)) + ALNG - ALGAP1
-      IF (T.GT.ALNEPS) GO TO 60
-C
-      IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * EXP(T)
-      IF (ABS(H).GT.SQEPS) GO TO 50
-C
-      CALL XERCLR
-      CALL XERMSG ('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1,
-     +   1)
-C
-C 50   T = -A*ALX + LOG(ABS(H))
-C      IF (T.LT.BOT) CALL XERCLR
-C      DGAMIT = SIGN (EXP(T), H)
- 50   result = H
-      RETURN
-C
-C 60   T = T - A*ALX
- 60   IF (T.LT.BOT) CALL XERCLR
-      result = -SGA * SGNGAM * EXP(T)
-      RETURN
-
-      endif
-      return
-      end
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/cruft/slatec-fn/xsgmainc.f
--- a/liboctave/cruft/slatec-fn/xsgmainc.f	Mon Mar 20 15:32:08 2017 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,100 +0,0 @@
-      subroutine xsgammainc (a, x, result)
-
-c -- jwe, based on GAMIT.
-c
-c -- Do a better job than gami for large values of x.
-
-      real a, x, result
-      intrinsic exp, log, sqrt, sign, aint
-      external gami, alngam, r9lgit, r9lgic, r9gmit
-
-C     external gamr
-C     real GAMR
-
-      REAL AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX,
-     $     BOT, H, SGA, SGNGAM, SQEPS, T, R1MACH, R9GMIT,
-     $     R9LGIC, R9LGIT, ALNGAM, GAMI
-
-      LOGICAL FIRST
-
-      SAVE ALNEPS, SQEPS, BOT, FIRST
-
-      DATA FIRST /.TRUE./
-
-      if (x .eq. 0.0e0) then
-
-        if (a .eq. 0.0e0) then
-          result = 1.0e0
-        else
-          result = 0.0e0
-        endif
-
-      else
-
-      IF (FIRST) THEN
-         ALNEPS = -LOG (R1MACH(3))
-         SQEPS = SQRT(R1MACH(4))
-         BOT = LOG (R1MACH(1))
-      ENDIF
-      FIRST = .FALSE.
-C
-      IF (X .LT. 0.E0) CALL XERMSG ('SLATEC', 'XGMAINC', 'X IS NEGATIVE'
-     +   , 2, 2)
-C
-      IF (X.NE.0.E0) ALX = LOG (X)
-      SGA = 1.0E0
-      IF (A.NE.0.E0) SGA = SIGN (1.0E0, A)
-      AINTA = AINT (A + 0.5E0*SGA)
-      AEPS = A - AINTA
-C
-C      IF (X.GT.0.E0) GO TO 20
-C      GAMIT = 0.0E0
-C      IF (AINTA.GT.0.E0 .OR. AEPS.NE.0.E0) GAMIT = GAMR(A+1.0E0)
-C      RETURN
-C
- 20   IF (X.GT.1.E0) GO TO 30
-      IF (A.GE.(-0.5E0) .OR. AEPS.NE.0.E0) CALL ALGAMS (A+1.0E0, ALGAP1,
-     1  SGNGAM)
-C      GAMIT = R9GMIT (A, X, ALGAP1, SGNGAM, ALX)
-      result = exp (a*alx + log (R9GMIT (A, X, ALGAP1, SGNGAM, ALX)))
-      RETURN
-C
- 30   IF (A.LT.X) GO TO 40
-      T = R9LGIT (A, X, ALNGAM(A+1.0E0))
-      IF (T.LT.BOT) CALL XERCLR
-C      GAMIT = EXP (T)
-      result = EXP (a*alx + T)
-      RETURN
-C
- 40   ALNG = R9LGIC (A, X, ALX)
-C
-C EVALUATE GAMIT IN TERMS OF LOG (DGAMIC (A, X))
-C
-      H = 1.0E0
-      IF (AEPS.EQ.0.E0 .AND. AINTA.LE.0.E0) GO TO 50
-C
-      CALL ALGAMS (A+1.0E0, ALGAP1, SGNGAM)
-      T = LOG (ABS(A)) + ALNG - ALGAP1
-      IF (T.GT.ALNEPS) GO TO 60
-C
-      IF (T.GT.(-ALNEPS)) H = 1.0E0 - SGA * SGNGAM * EXP(T)
-      IF (ABS(H).GT.SQEPS) GO TO 50
-C
-      CALL XERCLR
-      CALL XERMSG ('SLATEC', 'XGMAINC', 'RESULT LT HALF PRECISION', 1,
-     +   1)
-C
-C 50   T = -A*ALX + LOG(ABS(H))
-C      IF (T.LT.BOT) CALL XERCLR
-C      GAMIT = SIGN (EXP(T), H)
- 50   result = H
-      RETURN
-C
-C 60   T = T - A*ALX
- 60   IF (T.LT.BOT) CALL XERCLR
-      result = -SGA * SGNGAM * EXP(T)
-      RETURN
-
-      endif
-      return
-      end
diff -r d05f688ae836 -r 4ab8551eb58e liboctave/numeric/lo-specfun.cc
--- a/liboctave/numeric/lo-specfun.cc	Mon Mar 20 15:32:08 2017 -0700
+++ b/liboctave/numeric/lo-specfun.cc	Sun Feb 05 00:30:18 2017 +0100
@@ -2584,334 +2584,6 @@ namespace octave
       return retval;
     }
 
-    // FIXME: there is still room for improvement here...
-
-    double
-    gammainc (double x, double a, bool& err)
-    {
-      if (a < 0.0 || x < 0.0)
-        (*current_liboctave_error_handler)
-          ("gammainc: A and X must be non-negative");
-
-      err = false;
-
-      double retval;
-
-      F77_XFCN (xgammainc, XGAMMAINC, (a, x, retval));
-
-      return retval;
-    }
-
-    Matrix
-    gammainc (double x, const Matrix& a)
-    {
-      octave_idx_type nr = a.rows ();
-      octave_idx_type nc = a.cols ();
-
-      Matrix retval (nr, nc);
-
-      bool err;
-
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            retval(i,j) = gammainc (x, a(i,j), err);
-
-            if (err)
-              return Matrix ();
-          }
-
-      return retval;
-    }
-
-    Matrix
-    gammainc (const Matrix& x, double a)
-    {
-      octave_idx_type nr = x.rows ();
-      octave_idx_type nc = x.cols ();
-
-      Matrix retval (nr, nc);
-
-      bool err;
-
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            retval(i,j) = gammainc (x(i,j), a, err);
-
-            if (err)
-              return Matrix ();
-          }
-
-      return retval;
-    }
-
-    Matrix
-    gammainc (const Matrix& x, const Matrix& a)
-    {
-      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)
-        (*current_liboctave_error_handler)
-          ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
-           nr, nc, a_nr, a_nc);
-
-      Matrix retval (nr, nc);
-
-      bool err;
-
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            retval(i,j) = gammainc (x(i,j), a(i,j), err);
-
-            if (err)
-              return Matrix ();
-          }
-
-      return retval;
-    }
-
-    NDArray
-    gammainc (double x, const NDArray& a)
-    {
-      dim_vector dv = a.dims ();
-      octave_idx_type nel = dv.numel ();
-
-      NDArray retval (dv);
-
-      bool err;
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        {
-          retval(i) = gammainc (x, a(i), err);
-
-          if (err)
-            return NDArray ();
-        }
-
-      return retval;
-    }
-
-    NDArray
-    gammainc (const NDArray& x, double a)
-    {
-      dim_vector dv = x.dims ();
-      octave_idx_type nel = dv.numel ();
-
-      NDArray retval (dv);
-
-      bool err;
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        {
-          retval(i) = gammainc (x(i), a, err);
-
-          if (err)
-            return NDArray ();
-        }
-
-      return retval;
-    }
-
-    NDArray
-    gammainc (const NDArray& x, const NDArray& a)
-    {
-      dim_vector dv = x.dims ();
-      octave_idx_type nel = dv.numel ();
-
-      if (dv != a.dims ())
-        {
-          std::string x_str = dv.str ();
-          std::string a_str = a.dims ().str ();
-
-          (*current_liboctave_error_handler)
-            ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
-             x_str.c_str (), a_str.c_str ());
-        }
-
-      NDArray retval (dv);
-
-      bool err;
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        {
-          retval(i) = gammainc (x(i), a(i), err);
-
-          if (err)
-            return NDArray ();
-        }
-
-      return retval;
-    }
-
-    float
-    gammainc (float x, float a, bool& err)
-    {
-      if (a < 0.0 || x < 0.0)
-        (*current_liboctave_error_handler)
-          ("gammainc: A and X must be non-negative");
-
-      err = false;
-
-      float retval;
-
-      F77_XFCN (xsgammainc, XSGAMMAINC, (a, x, retval));
-
-      return retval;
-    }
-
-    FloatMatrix
-    gammainc (float x, const FloatMatrix& a)
-    {
-      octave_idx_type nr = a.rows ();
-      octave_idx_type nc = a.cols ();
-
-      FloatMatrix retval (nr,  nc);
-
-      bool err;
-
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            retval(i,j) = gammainc (x, a(i,j), err);
-
-            if (err)
-              return FloatMatrix ();
-          }
-
-      return retval;
-    }
-
-    FloatMatrix
-    gammainc (const FloatMatrix& x, float a)
-    {
-      octave_idx_type nr = x.rows ();
-      octave_idx_type nc = x.cols ();
-
-      FloatMatrix retval (nr, nc);
-
-      bool err;
-
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            retval(i,j) = gammainc (x(i,j), a, err);
-
-            if (err)
-              return FloatMatrix ();
-          }
-
-      return retval;
-    }
-
-    FloatMatrix
-    gammainc (const FloatMatrix& x, const FloatMatrix& a)
-    {
-      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)
-        (*current_liboctave_error_handler)
-          ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
-           nr, nc, a_nr, a_nc);
-
-      FloatMatrix retval (nr, nc);
-
-      bool err;
-
-      for (octave_idx_type j = 0; j < nc; j++)
-        for (octave_idx_type i = 0; i < nr; i++)
-          {
-            retval(i,j) = gammainc (x(i,j), a(i,j), err);
-
-            if (err)
-              return FloatMatrix ();
-          }
-
-      return retval;
-    }
-
-    FloatNDArray
-    gammainc (float x, const FloatNDArray& a)
-    {
-      dim_vector dv = a.dims ();
-      octave_idx_type nel = dv.numel ();
-
-      FloatNDArray retval (dv);
-
-      bool err;
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        {
-          retval(i) = gammainc (x, a(i), err);
-
-          if (err)
-            return FloatNDArray ();
-        }
-
-      return retval;
-    }
-
-    FloatNDArray
-    gammainc (const FloatNDArray& x, float a)
-    {
-      dim_vector dv = x.dims ();
-      octave_idx_type nel = dv.numel ();
-
-      FloatNDArray retval (dv);
-
-      bool err;
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        {
-          retval(i) = gammainc (x(i), a, err);
-
-          if (err)
-            return FloatNDArray ();
-        }
-
-      return retval;
-    }
-
-    FloatNDArray
-    gammainc (const FloatNDArray& x, const FloatNDArray& a)
-    {
-      dim_vector dv = x.dims ();
-      octave_idx_type nel = dv.numel ();
-
-      if (dv != a.dims ())
-        {
-          std::string x_str = dv.str ();
-          std::string a_str = a.dims ().str ();
-
-          (*current_liboctave_error_handler)
-            ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
-             x_str.c_str (), a_str.c_str ());
-        }
-
-      FloatNDArray retval (dv);
-
-      bool err;
-
-      for (octave_idx_type i = 0; i < nel; i++)
-        {
-          retval(i) = gammainc (x(i), a(i), err);
-
-          if (err)
-            return FloatNDArray ();
-        }
-
-      return retval;
-    }
-
     Complex rc_log1p (double x)
     {
       const double pi = 3.14159265358979323846;
@@ -3996,22 +3668,6 @@ Array<float> betainc (const Array<float>
 Array<float> betainc (const Array<float>& x, const Array<float>& a, float b) { return octave::math::betainc (x, a, b); }
 Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b) { return octave::math::betainc (x, a, b); }
 
-Matrix gammainc (double x, const Matrix& a) { return octave::math::gammainc (x, a); }
-Matrix gammainc (const Matrix& x, double a) { return octave::math::gammainc (x, a); }
-Matrix gammainc (const Matrix& x, const Matrix& a) { return octave::math::gammainc (x, a); }
-
-NDArray gammainc (double x, const NDArray& a) { return octave::math::gammainc (x, a); }
-NDArray gammainc (const NDArray& x, double a) { return octave::math::gammainc (x, a); }
-NDArray gammainc (const NDArray& x, const NDArray& a) { return octave::math::gammainc (x, a); }
-
-FloatMatrix gammainc (float x, const FloatMatrix& a) { return octave::math::gammainc (x, a); }
-FloatMatrix gammainc (const FloatMatrix& x, float a) { return octave::math::gammainc (x, a); }
-FloatMatrix gammainc (const FloatMatrix& x, const FloatMatrix& a) { return octave::math::gammainc (x, a); }
-
-FloatNDArray gammainc (float x, const FloatNDArray& a) { return octave::math::gammainc (x, a); }
-FloatNDArray gammainc (const FloatNDArray& x, float a) { return octave::math::gammainc (x, a); }
-FloatNDArray gammainc (const FloatNDArray& x, const FloatNDArray& a) { return octave::math::gammainc (x, a); }
-
 Array<double> betaincinv (double x, double a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); }
 Array<double> betaincinv (double x, const Array<double>& a, double b) { return octave::math::betaincinv (x, a, b); }
 Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); }
diff -r d05f688ae836 -r 4ab8551eb58e scripts/specfun/gammainc.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/specfun/gammainc.m	Sun Feb 05 00:30:18 2017 +0100
@@ -0,0 +1,420 @@
+## Copyright (C) 2017 Marco Caliari
+## Copyright (C) 2017 Nir Krakauer
+##
+## 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/>.
+
+## -*- texinfo -*-
+## @deftypefn  {} {} gammainc (@var{x}, @var{a})
+## @deftypefnx {} {} gammainc (@var{x}, @var{a}, "lower")
+## @deftypefnx {} {} gammainc (@var{x}, @var{a}, "upper")
+## @deftypefnx {} {} gammainc (@var{x}, @var{a}, "scaledlower")
+## @deftypefnx {} {} gammainc (@var{x}, @var{a}, "scaledupper")
+## Compute the normalized incomplete gamma function.
+##
+## This is defined as
+## @tex
+## $$
+##  \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt}
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##                                 x
+##                        1       /
+## gammainc (x, a) = ---------    | exp (-t) t^(a-1) dt
+##                   gamma (a)    /
+##                             t=0
+## @end group
+## @end example
+##
+## @end ifnottex
+## with the limiting value of 1 as @var{x} approaches infinity.
+## The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and
+## @nospell{Stegun} (6.5.1).
+##
+## If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned
+## for each element of @var{x} and vice versa.
+##
+## If neither @var{x} nor @var{a} is scalar, the sizes of @var{x} and
+## @var{a} must agree, and @code{gammainc} is applied element-by-element.
+##
+## By default the incomplete gamma function integrated from 0 to @var{x} is
+## computed.  If @qcode{"upper"} is given then the complementary function
+## integrated from @var{x} to infinity is calculated.
+##
+## Algorithms used include the following.
+## Numerical Recipes in Fortran77, \S 6.2.
+## Regions roughly as in Gautschi, W. (1979).
+## A computational procedure for incomplete gamma functions.
+## ACM Transactions on Mathematical Software (TOMS), 5(4), 466-481.
+##
+## @seealso{gamma, gammaln}
+## @end deftypefn
+
+## P(a,x) = gamma(a,x)/Gamma(a), upper
+## 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a), lower
+
+function y = gammainc (x,a,tail)
+
+  if (nargin >= 4 || nargin <= 1)
+    print_usage ();
+  endif
+  if (nargin == 2)
+    tail = "lower";
+  endif
+
+  if (! isscalar (x) || ! isscalar (a))
+    [err, x, a] = common_size (x, a);
+    if (err > 0)
+      error ("gammainc: x, a must be of common size or scalars");
+    endif
+  endif
+
+  ## Initialize output array.
+  if (isinteger (x))
+    x = double (x);
+  endif
+
+  if (strcmpi (class (a), "single"))
+    x = single (x);
+  endif
+
+  y = zeros (size (x), class (x));
+
+  ## Different x, a combinations are handled by different subfunctions.
+  i_done = false (size (x)); # Track which elements have been calculated.
+  
+  ## Case 1: x == 0, a == 0
+  ii = (x == 0) & (a == 0);
+  if (any (ii(:)))
+    y(ii) = gammainc_00 (tail);
+    i_done(ii) = true;
+  endif
+  
+  ## Case 2: x == 0
+  ii = ~i_done & (x == 0);
+  if (any (ii(:)))
+    y(ii) = gammainc_x0 (tail);
+    i_done(ii) = true;
+  endif
+
+  ## Case 3: a == 0
+  ii = ~i_done & (a == 0);
+  if (any (ii(:)))
+    y(ii) = gammainc_a0 (x(ii), tail);
+    i_done(ii) = true;
+  endif
+
+  ## Case 4: a == 1
+  ii = ~i_done & (a == 1);
+  if (any (ii(:)))
+    y(ii) = gammainc_a1 (x(ii), tail);
+    i_done(ii) = true;
+  endif
+
+  ## Case 5: positive integer a; exp (x) and a! both under 1/eps
+  ii = ~i_done & (a > 1) & (a == fix (a)) & (x <= 36) & (a <= 18);
+  if (any (ii(:)))
+    y(ii) = gammainc_an (x(ii), a(ii), tail);
+    i_done(ii) = true;
+  endif
+
+  ## Case 6: x relatively small
+  ii = ~i_done & (x + 0.25 < a | x < 0);
+  if (any (ii(:)))
+    y(ii) = gammainc_s (x(ii), a(ii), tail);
+    i_done(ii) = true;
+  endif
+
+  ## Case 7: x positive and large relative to a
+  ii = ~i_done;
+  if (any (ii(:)))
+    y(ii) = gammainc_l (x(ii), a(ii), tail);
+    i_done(ii) = true;
+  endif
+
+endfunction
+
+
+
+
+  ##subfunctions to handle each case:
+  
+  ## x == 0, a == 0
+  function y = gammainc_00 (tail)
+    if (strcmpi (tail, "upper") || strcmpi (tail, "scaledupper"))
+      y = 0;
+    else
+      y = 1;
+    endif
+  endfunction
+
+  ## x == 0, a != 0
+  function y = gammainc_x0 (tail)
+    if (strcmpi (tail, "upper") || strcmpi (tail, "scaledlower"))
+      y = 1;
+    elseif (strcmpi (tail, "lower"))
+      y = 0;
+    else
+      y = Inf;
+    endif
+  endfunction
+
+  ## x !=0, a == 0
+  function y = gammainc_a0 (x, tail)
+    if strcmpi (tail, "lower")
+      y = 1;
+    elseif strcmpi (tail, "scaledlower")
+      y = exp (x);
+    else
+      y = 0;
+    endif
+  endfunction
+
+  ## a == 1
+  function y = gammainc_a1 (x, tail)
+    if strcmpi (tail, "lower")
+      y = 1 - exp (-x);
+    elseif strcmpi (tail, "scaledlower")
+      if (abs (x) < 1/2)
+        y = expm1 (x) ./ x;
+      else
+        y = (exp (x) - 1) ./ x;
+      endif
+    elseif strcmpi (tail, "upper")
+      y = exp (-x);
+    else
+      y = 1 ./ x;
+    endif
+  endfunction
+
+  ## positive integer a; exp (x) and a! both under 1/eps
+  ## uses closed-form expressions for nonnegative integer a
+  ## -- http://mathworld.wolfram.com/IncompleteGammaFunction.html
+  function y = gammainc_an (x, a, tail)
+    y = t = ones (size (x), class (x));
+    i = 1;
+    while any (a(:) > i)
+      jj = (a > i);
+      t(jj) .*= (x(jj) / i);
+      y(jj) += t(jj);
+      i++;
+    endwhile
+    if strcmpi (tail, "upper")
+      y .*= exp (-x);
+    elseif strcmpi (tail, "lower")
+      y = 1 - exp (-x) .* y;
+    elseif strcmpi (tail, "scaledupper")
+      y .*= D (x,a) .* exp (-x);
+    elseif strcmpi (tail, "lower")
+      y = D (x,a) .* (1 - exp (-x) .* y);
+    endif
+  endfunction
+
+  ## x + 0.25 < a | x < 0
+  ## Numerical Recipes in Fortran 77 (6.2.5)
+  ## series
+  function y = gammainc_s (x, a, tail)
+    if (strcmpi (tail, "scaledlower") || strcmpi (tail, "scaledupper"))
+      y = ones (size (x), class (x));
+      term = x ./ (a + 1);
+    else
+      ## Of course it is possible to scale at the end, but some tests fail.
+      ## And try gammainc (1,1000), it take 0 iterations if you scale now.
+      y = D (x,a);
+      term = y .* x ./ (a + 1);
+    endif
+    n = 1;
+    while any (abs (term(:)) > abs (y(:)) * eps)
+      ## y can be zero from the beginning (gammainc (1,1000))
+      jj = abs (term) > abs (y) * eps;
+      n += 1;
+      y(jj) += term(jj);
+      term(jj) .*= x(jj) ./ (a(jj) + n);
+    endwhile
+    if strcmpi (tail, "upper")
+      y = 1 - y;
+    elseif strcmpi (tail, "scaledupper")
+      y = 1 ./ D (x,a) - y;
+    endif
+  endfunction
+
+  ## x positive and large relative to a
+  ## NRF77 (6.2.7)
+  ## Gamma (a,x)/Gamma (a)
+  ## Lentz's algorithm
+
+  function y = gammainc_l (x, a, tail)
+    n = numel (x);
+    y = zeros (size (x), class (x));
+    for i = 1:n
+      y(i) = __gammainc_lentz__ (x(i), a(i));
+    endfor
+    if strcmpi (tail, "upper")
+      y .*= D (x,a);
+    elseif strcmpi (tail,  "lower")
+      y = 1 - y .* D (x,a);
+    elseif strcmpi (tail, "scaledlower")
+      y = 1 ./ D (x,a) - y;
+    endif
+  endfunction
+
+function y = D(x, a)
+  ##  Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large
+  ##
+  ## L. Knusel, Computation of the Chi-square and Poisson distribution,
+  ## SIAM J. Sci. Stat. Comput., 7(3), 1986
+  ## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41
+  athresh = 10; ## FIXME: can be better tuned?
+  y = zeros(size(x), class (x));
+  i_done = false (size (x));
+  i_done(x == 0) = true;
+  ii = ~i_done & (x > 0) & (a > athresh) & (a >= x);
+  if any (ii(:))
+    a_tmp = a(ii);
+    lnGa = log (2 * pi * a_tmp) / 2 + 1 ./ (12 * a_tmp) - ...
+           1 ./ (360 * a_tmp .^ 3) + 1 ./ (1260 * a_tmp .^ 5) - ...
+           1 ./ (1680 * a_tmp .^ 7) + 1 ./ (1188 * a_tmp .^ 9)- ...
+           691 ./ (87360 * a_tmp .^ 11) + 1 ./ (156 * a_tmp .^ 13) - ...
+           3617 ./ (122400 * a_tmp .^ 15) + ...
+           43867 ./ (244188 * a_tmp .^ 17) - 174611 ./ (125400 * a_tmp .^ 19);
+    lns = log1p ((a(ii) - x(ii)) ./ x(ii));
+    y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa);
+    i_done(ii) = true;
+  endif
+  ii = ~i_done & (x > 0) & (a > athresh) & (a < x);
+  if any (ii(:))
+    a_tmp = a(ii);
+    lnGa = log (2 * pi * a_tmp) / 2 + 1 ./ (12 * a_tmp) - ...
+           1 ./ (360 * a_tmp .^ 3) + 1 ./ (1260 * a_tmp .^ 5) - ...
+           1 ./ (1680 * a_tmp .^ 7) + 1 ./ (1188 * a_tmp .^ 9)- ...
+           691 ./ (87360 * a_tmp .^ 11) + 1 ./ (156 * a_tmp .^ 13) - ...
+           3617 ./ (122400 * a_tmp .^ 15) + ...
+           43867 ./ (244188 * a_tmp .^ 17) - 174611 ./ (125400 * a_tmp .^ 19);
+    lns = -log1p ((x(ii) - a_tmp) ./ a_tmp);
+    y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa);
+    i_done(ii) = true;
+  endif
+  ii = ~i_done & ((x <= 0) | (a <= athresh));
+  if any (ii(:)) ## standard formula for a not so large
+    y(ii) = exp (a(ii) .* log (x(ii)) - x(ii) - gammaln (a(ii) + 1));
+    i_done(ii) = true;
+  endif
+  ii = (x < 0) & (a == fix (a));
+  if any(ii(:)) ## remove spurious imaginary part
+    y(ii) = real (y(ii));
+  endif
+endfunction
+
+%!test
+%! assert (gammainc(0, [0, 1, 2, 3]), [1, 0, 0, 0]);
+%!test
+%! assert (gammainc(0, [0, 1, 2, 3], "upper"), [0, 1, 1, 1]);
+%!test
+%! assert (gammainc(0, [0, 1, 2, 3], "scaledupper"), [0, Inf, Inf, Inf]);
+%!test
+%! assert (gammainc(0, [0, 1, 2, 3], "scaledlower"), [1, 1, 1, 1]);
+%!test
+%! assert (gammainc([0 0],0),[1 1]);
+%!test
+%! assert (gammainc (0,[0; 0],"upper"),[0; 0]);
+%!test
+%!error (gammainc([0 0],[0; 0]));
+%!test
+%! assert (gammainc ([1 10],[1000 0]),[0 1]);
+%!test
+%! assert (gammainc ([3 2 36],[2 3 18], "upper"),...
+%! [4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)]);
+%!test
+%! assert (gammainc (10,10),1 - (5719087 / 567) * exp (-10),-eps);
+%!test
+%! assert (gammainc (10,10, "upper"),(5719087 / 567) * exp (-10),-eps);
+%!test
+%! assert (gammainc (-10,10),3.112658265341493126871617e7,-2 * eps);
+%!test
+%! % Here Matlab fails
+%! assert (isreal (gammainc (-10,10)),true);
+%!test
+%! assert (gammainc (60,6,"upper"),6.18022358081160257327264261e-20,-10 * eps);
+%!test
+%! assert (gammainc (-10,10.1,"upper"),...
+%!        -2.9582761911890713293e7-1i * 9.612022339061679758e6,-30 * eps);
+%!test
+%! % Here Matlab is better
+%! assert (gammainc (751,750,"upper"),...
+%! 0.4805914320558831327179457887,-12 * eps);
+%!test
+%! % Here Matlab fails
+%! assert (gammainc (-100,1,"upper"),exp (100),-eps);
+%!test
+%! assert (gammainc (200,200,"upper"),0.49059658199276367497217454,-4 * eps);
+%!test
+%! assert (gammainc (200,200),0.509403418007236325027825459574527043,-3 * eps);
+%!test
+%! assert (gammainc (200,200,"scaledupper"),17.3984438553791505135122900,-eps);
+%!test
+%! assert (gammainc (200,200,"scaledlower"),18.065406676779221643065,-6 * eps);
+%!test
+%! assert (gammainc (200,201,"upper"),0.5187943096786844967708766027,-2 * eps);
+%!test
+%! assert (gammainc(201,200,"upper"),0.46249244908276709524913736667,-7 * eps);
+%!test
+%! assert (gammainc (200,201,"scaledupper"),18.4904360746560462660798514,-eps);
+%!test
+%! % here we are very good (no D (x,a)) involved
+%! assert (gammainc(1000,1000.5,"scaledlower"),39.48467539583672271,-2* eps);
+%!test
+%! assert (gammainc (709,1000,"upper"), 0.99999999999999999999999954358, -eps);
+
+## Test on the conservation of the class, five tests for each subroutine
+%!assert (class (gammainc (0, 1)) == "double")
+%!assert (class (gammainc (single (0), 1)) == "single")
+%!assert (class (gammainc (int8 (0), 1)) == "double")
+%!assert (class (gammainc (0, single (1))) == "single")
+%!assert (class (gammainc (0, int8 (1))) == "double")
+%!assert (class (gammainc (1, 0)) == "double")
+%!assert (class (gammainc (single (1), 0)) == "single")
+%!assert (class (gammainc (int8 (1), 0)) == "double")
+%!assert (class (gammainc (1, single (0))) == "single")
+%!assert (class (gammainc (1, int8 (0))) == "double")
+%!assert (class (gammainc (1, 1)) == "double")
+%!assert (class (gammainc (single (1), 1)) == "single")
+%!assert (class (gammainc (int8 (1), 1)) == "double")
+%!assert (class (gammainc (1, single (1))) == "single")
+%!assert (class (gammainc (1, int8(1))) == "double")
+%!assert (class (gammainc (1, 2)) == "double")
+%!assert (class (gammainc (single (1), 2)) == "single")
+%!assert (class (gammainc (int8 (1), 2)) == "double")
+%!assert (class (gammainc (1, single (2))) == "single")
+%!assert (class (gammainc (1, int8 (2))) == "double")
+%!assert (class (gammainc (-1, 0.5)) == "double")
+%!assert (class (gammainc (single (-1), 0.5)) == "single")
+%!assert (class (gammainc (int8 (-1), 0.5)) == "double")
+%!assert (class (gammainc (-1, single (0.5))) == "single")
+%!assert (class (gammainc (-1, int8 (0.5))) == "double")
+%!assert (class (gammainc (1, 0.5)) == "double")
+%!assert (class (gammainc (single (1), 0.5)) == "single")
+%!assert (class (gammainc (int8 (1), 0.5)) == "double")
+%!assert (class (gammainc (1, single (0.5))) == "single")
+%!assert (class (gammainc (1, int8 (0.5))) == "double")
+
+## Test input validation
+%!error gammainc()
+%!error gammainc(1)
+%!error gammainc(1,2,3,4)
diff -r d05f688ae836 -r 4ab8551eb58e scripts/specfun/module.mk
--- a/scripts/specfun/module.mk	Mon Mar 20 15:32:08 2017 -0700
+++ b/scripts/specfun/module.mk	Sun Feb 05 00:30:18 2017 +0100
@@ -8,6 +8,7 @@ scripts_specfun_FCN_FILES = \
   scripts/specfun/expint.m \
   scripts/specfun/factor.m \
   scripts/specfun/factorial.m \
+    scripts/specfun/gammainc.m \
   scripts/specfun/isprime.m \
   scripts/specfun/lcm.m \
   scripts/specfun/legendre.m \
