Mercurial > octave
changeset 24904:740df3fce3fb stable
New strategy to compute the incomplete gamma function (see bug #47800).
--
added libinterp/corefcn/__gammainc_lentz__.cc
added scripts/specfun/gammainc.m
changed libinterp/corefcn/module.mk
changed liboctave/external/slatec-fn/module.mk
changed liboctave/numeric/lo-specfun.cc
changed scripts/specfun/module.mk
removed liboctave/external/slatec-fn/dgami.f
removed liboctave/external/slatec-fn/dgamit.f
removed liboctave/external/slatec-fn/gami.f
removed liboctave/external/slatec-fn/gamit.f
removed liboctave/external/slatec-fn/xdgami.f
removed liboctave/external/slatec-fn/xdgamit.f
removed liboctave/external/slatec-fn/xgmainc.f
removed liboctave/external/slatec-fn/xsgmainc.f
| author | Michele Ginesi <michele.ginesi@gmail.com> |
|---|---|
| date | Thu, 07 Sep 2017 15:47:33 +0200 |
| parents | 8862b95d64ec |
| children | 662faf9de127 |
| files | libinterp/corefcn/__gammainc_lentz__.cc libinterp/corefcn/module.mk liboctave/external/slatec-fn/dgami.f liboctave/external/slatec-fn/dgamit.f liboctave/external/slatec-fn/gami.f liboctave/external/slatec-fn/gamit.f liboctave/external/slatec-fn/module.mk liboctave/external/slatec-fn/xdgami.f liboctave/external/slatec-fn/xdgamit.f liboctave/external/slatec-fn/xgmainc.f liboctave/external/slatec-fn/xsgmainc.f liboctave/numeric/lo-specfun.cc scripts/specfun/gammainc.m scripts/specfun/module.mk |
| diffstat | 14 files changed, 687 insertions(+), 888 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__gammainc_lentz__.cc Thu Sep 07 15:47:33 2017 +0200 @@ -0,0 +1,92 @@ +// Copyright (C) 2017 Nir Krakauer +// 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 (__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; + int maxit = 200; + while((std::abs((Deltaj - 1) / y) > eps) & (j < maxit)) + { + 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; + int maxit = 200; + while((std::abs((Deltaj - 1) / y) > eps) & (j < maxit)) + { + 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 (); +}
--- a/libinterp/corefcn/module.mk Sun Mar 18 18:42:35 2018 -0700 +++ b/libinterp/corefcn/module.mk Thu Sep 07 15:47:33 2017 +0200 @@ -109,6 +109,7 @@ %reldir%/Cell.cc \ %reldir%/__contourc__.cc \ %reldir%/__dsearchn__.cc \ + %reldir%/__gammainc_lentz__.cc \ %reldir%/__ichol__.cc \ %reldir%/__ilu__.cc \ %reldir%/__lin_interpn__.cc \ @@ -155,7 +156,6 @@ %reldir%/filter.cc \ %reldir%/find.cc \ %reldir%/ft-text-renderer.cc \ - %reldir%/gammainc.cc \ %reldir%/gcd.cc \ %reldir%/getgrent.cc \ %reldir%/getpwent.cc \
--- a/liboctave/external/slatec-fn/dgami.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/external/slatec-fn/dgamit.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/external/slatec-fn/gami.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/external/slatec-fn/gamit.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/external/slatec-fn/module.mk Sun Mar 18 18:42:35 2018 -0700 +++ b/liboctave/external/slatec-fn/module.mk Thu Sep 07 15:47:33 2017 +0200 @@ -11,8 +11,6 @@ %reldir%/d9lgmc.f \ %reldir%/dbetai.f \ %reldir%/dcsevl.f \ - %reldir%/dgami.f \ - %reldir%/dgamit.f \ %reldir%/dgamlm.f \ %reldir%/dgamma.f \ %reldir%/dgamr.f \ @@ -23,8 +21,6 @@ %reldir%/dpchim.f \ %reldir%/dpchst.f \ %reldir%/dpsifn.f \ - %reldir%/gami.f \ - %reldir%/gamit.f \ %reldir%/gamlim.f \ %reldir%/gamma.f \ %reldir%/gamr.f \ @@ -38,10 +34,6 @@ %reldir%/r9gmit.f \ %reldir%/r9lgic.f \ %reldir%/xdbetai.f \ - %reldir%/xdgami.f \ - %reldir%/xdgamit.f \ - %reldir%/xgmainc.f \ - %reldir%/xsgmainc.f \ %reldir%/xbetai.f DIRSTAMP_FILES += %reldir%/$(octave_dirstamp)
--- a/liboctave/external/slatec-fn/xdgami.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/external/slatec-fn/xdgamit.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/external/slatec-fn/xgmainc.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/external/slatec-fn/xsgmainc.f Sun Mar 18 18:42:35 2018 -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
--- a/liboctave/numeric/lo-specfun.cc Sun Mar 18 18:42:35 2018 -0700 +++ b/liboctave/numeric/lo-specfun.cc Thu Sep 07 15:47:33 2017 +0200 @@ -2829,334 +2829,6 @@ return result; } - // 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 log1p (const Complex& x) { @@ -3593,22 +3265,6 @@ 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); }
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/gammainc.m Thu Sep 07 15:47:33 2017 +0200 @@ -0,0 +1,593 @@ +## Copyright (C) 2016 Marco Caliari +## Copyright (C) 2016 Nir Krakauer +## 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/>. + +## -*- texinfo -*- +## @deftypefn {} {} gammainc (@var{x}, @var{a}) +## @deftypefnx {} {} gammainc (@var{x}, @var{a}, @var{tail}) +## 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. +## The elements of @var{a} must be nonnegative. +## +## By default or if @var{tail} is @qcode{"lower"} the incomplete gamma +## function integrated from 0 to @var{x} is computed. If @var{tail} +## is @qcode{"upper"} then the complementary function integrated from +## @var{x} to infinity is calculated. +## +## If @var{tail} is @qcode{"scaledlower"}, then the lower incomplete gamma +## function is multiplied by +## @tex +## $\Gamma(a+1)\exp(x)x^{-a}$. +## @end tex +## @ifnottex +## @math{gamma(a+1)*exp(x)/(x^a)}. +## @end ifnottex +## If @var{tail} is @qcode{"scaledupper"}, then the upper incomplete gamma +## function is divided by the same quantity. +## +## References: +## +## @nospell{M. Abramowitz and I. Stegun}, +## @cite{Handbook of mathematical functions} +## @nospell{Dover publications, INC.}, +## 1972. +## +## @nospell{W. Gautschi}, +## @cite{A computational procedure for incomplete gamma functions} +## @nospell{ACM Trans. Math Software}, +## pp. 466--481, Vol 5, No. 4, 2012. +## +## @nospell{W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery} +## @cite{Numerical Recipes in Fortran 77}, +## ch. 6.2, Vol 1, 1992. +## +## @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 = "lower") + + if (nargin >= 4 || nargin <= 1) + print_usage (); + 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 + + if (any (a < 0) || any (imag (a) != 0)) + error ("gammainc: a must be real and non negative"); + endif + + if (any (imag (x) != 0)) + error ("gammainc: x must be real"); + 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 0: x == Inf, a == Inf + + ii = ((x == Inf) & (a == Inf)); + if (any (ii(:))) + y(ii) = NaN; + i_done(ii) = true; + endif + + ## 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: x = Inf + ii = ((! i_done) & (x == Inf)); + if (any (ii(:))) + y(ii) = gammainc_x_inf (tail); + i_done(ii) = true; + endif + + ## Case 4: a = Inf + ii = ((! i_done) & (a == Inf)); + if (any (ii(:))) + y(ii) = gammainc_a_inf (tail); + i_done(ii) = true; + endif + + ## Case 5: a == 0. + ii = ((! i_done) & (a == 0)); + if (any (ii(:))) + y(ii) = gammainc_a0 (x(ii), tail); + i_done(ii) = true; + endif + + ## Case 6: a == 1. + ii = ((! i_done) & (a == 1)); + if (any (ii(:))) + y(ii) = gammainc_a1 (x(ii), tail); + i_done(ii) = true; + endif + + flag_an = ((a > 1) & (a == fix (a)) & (x <= 36) & (abs (x) >= 1e-01) & ... + (a <= 18)); + + ## Case 7: positive integer a; exp (x) and a! both under 1/eps. + ii = ((! i_done) & flag_an); + if (any (ii(:))) + y(ii) = gammainc_an (x(ii), a(ii), tail); + i_done(ii) = true; + endif + + ## For a < 2, x < 0, we increment a by 2 and use a recurrence formula after + ## the computations. + + flag_a_small = ((abs (a) < 2) & (abs(a) > 0) & (! i_done) & (x < 0)); + a(flag_a_small) += 2; + + flag_s = (((x + 0.25 < a | x < 0 | a < 5) & (x > -20)) | (abs (x) < 1)); + + ## Case 8: x, a relatively small. + ii = ((! i_done) & flag_s); + if (any (ii(:))) + y(ii) = gammainc_s (x(ii), a(ii), tail); + i_done(ii) = true; + endif + + ## Case 9: 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 + + if (any (flag_a_small)) + if (strcmpi (tail, "lower")) + y(flag_a_small) += D (x(flag_a_small), a(flag_a_small) - 1) + ... + D (x(flag_a_small), a(flag_a_small) - 2); + elseif (strcmpi (tail, "upper")) + y(flag_a_small) -= D (x(flag_a_small), a(flag_a_small) - 1) + ... + D (x(flag_a_small), a(flag_a_small) - 2); + elseif (strcmpi (tail, "scaledlower")) + y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ... + (a(flag_a_small) .* (a(flag_a_small) - 1)) + (x(flag_a_small) ./ ... + (a(flag_a_small) - 1)) + 1; + elseif (strcmpi (tail, "scaledupper")) + y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ... + (a(flag_a_small) .* (a(flag_a_small) - 1)) - (x(flag_a_small) ./ ... + (a(flag_a_small) - 1)) - 1; + endif + 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. +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 == Inf. +function y = gammainc_x_inf (tail) + if (strcmpi (tail, "lower")) + y = 1; + elseif (strcmpi (tail, "upper") || strcmpi (tail, "scaledupper")) + y = 0; + else + y = Inf; + endif +endfunction + +## a == Inf. +function y = gammainc_a_inf (tail) + if (strcmpi (tail, "lower")) + y = 0; + elseif (strcmpi (tail, "upper") || strcmpi (tail, "scaledlower")) + y = 1; + else + y = Inf; + endif +endfunction + +## 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 .*= exp (-x) ./ D(x, a); + elseif (strcmpi (tail, "scaledlower")) + y = (1 - exp (-x) .* y) ./ D(x, a); + endif +endfunction + +## x + 0.25 < a | x < 0 | x not real | abs(x) < 1 | a < 5. +## 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 +## __gammainc_lentz__ in libinterp/corefcn/__gammainc_lentz__.cc +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(:))) + lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ... + 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ... + 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ... + 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ... + 3617 ./ (122400 * a(ii) .^ 15) + ... + 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 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(:))) + lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ... + 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ... + 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ... + 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ... + 3617 ./ (122400 * a(ii) .^ 15) + ... + 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19); + lns = -log1p ((x(ii) - a(ii)) ./ a(ii)); + 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: case 1,2,5 +%!test +%! assert (gammainc ([0, 0, 1], [0, 1, 0]), [1, 0, 1]); +%!test +%! assert (gammainc ([0, 0, 1], [0, 1, 0], "upper"), [0, 1, 0]); +%!test +%! assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledlower"), [1, 1, exp(1)]); +%!test +%! assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledupper"), [0, Inf, 0]); + +## Test: case 3,4 +%!test +%! assert (gammainc ([2, Inf], [Inf, 2]), [0, 1]); +%!test +%! assert (gammainc ([2, Inf], [Inf, 2], "upper"), [1, 0]); +%!test +%! assert (gammainc ([2, Inf], [Inf, 2], "scaledlower"), [1, Inf]); +%!test +%! assert (gammainc ([2, Inf], [Inf, 2], "scaledupper"), [Inf, 0]); + +## Test: case 5 +%!test +%! % Here Matlab fails +%! assert (gammainc (-100,1,"upper"),exp (100),-eps); + +## Test: case 6 +%!test +%! assert (gammainc ([1, 2, 3], 1), 1 - exp (-[1, 2, 3])); +%!test +%! assert (gammainc ([1, 2, 3], 1, "upper"), exp (- [1, 2, 3])); +%!test +%! assert (gammainc ([1, 2, 3], 1, "scaledlower"), ... +%! (exp ([1, 2, 3]) - 1) ./ [1, 2, 3]); +%! assert (gammainc ([1, 2, 3], 1, "scaledupper"), 1 ./ [1, 2, 3]); + +## Test: case 7 +%!test +%! assert (gammainc (2, 2, "lower"), 0.593994150290162, -2e-15); +%!test +%! assert (gammainc (2, 2, "upper"), 0.406005849709838, -2e-15); +%!test +%! assert (gammainc (2, 2, "scaledlower"), 2.194528049465325, -2e-15); +%!test +%! assert (gammainc (2, 2, "scaledupper"), 1.500000000000000, -2e-15); +%!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: case 8 +%!test +%! assert (gammainc (-10, 10), 3.112658265341493126871617e7, -2 * eps); +%!test +%! % Here Matlab fails +%! assert (isreal (gammainc (-10, 10)), true); +%!test +%! assert (gammainc (-10, 10.1, "upper"),... +%! -2.9582761911890713293e7-1i * 9.612022339061679758e6,-30 * eps); +%!test +%! assert (gammainc (-10, 10, "upper"), -3.112658165341493126871616e7, ... +%! -2 * eps); +%!test +%! assert (gammainc (-10, 10, "scaledlower"), ... +%! 0.5128019364747265,-1e-14); +%!test +%! assert (gammainc (-10, 10, "scaledupper"), ... +%! -0.5128019200000000, -1e-14); +%!test +%! assert (gammainc (200, 201, "upper"),... +%! 0.518794309678684497, -2 * 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: case 9 +%!test <47800> +%! assert (gammainc (60, 6, "upper"), ... +%! 6.18022358081160257327264261e-20, -10 * eps); +%!test +%! % Here Matlab is better +%! assert (gammainc (751, 750, "upper"),... +%! 0.4805914320558831327179457887, -12 * 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 (201, 200, "upper"), 0.46249244908276709524913736667,... +%! -7 * eps); + +## Test small argument +%!test +%! assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.1), ... +%! [0.33239840504050, 0.20972940370977, 0.10511370061022, ... +%! 0.041846517936723], 1e-13); + +%!test +%! assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.2), ... +%! [0.10891226058559, 0.043358823442178, 0.010891244210402, ... +%! 0.0017261458806785], 1e-13); + +%!test +%! assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 0.9), ... +%! [0.016401189184068, 0.0020735998660840, 0.000032879756964708, ... +%! 8.2590606569241e-9, 2.6117443021738e-13], -1e-12); + +%!test +%! assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 2), ... +%! [0.0000496679133402659, 4.99666791633340e-7, 4.99996666679167e-11, ... +%! 4.99999999666667e-19, 4.99999999999997e-29], -1e-12); + +%!xtest +%! assert (gammainc (-20, 1.1, "upper"),... +%! 6.50986687074979e8 + 2.11518396291149e8*i, -1e-13); + +## 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) +%!error gammainc (1, [0, 1i,1]) +%!error gammainc (1, [0, -1, 1]) +%!error gammainc ([0 0],[0; 0]) +%!error gammainc ([1 2 3], [1 2]) +%!error gammainc (1i,1) +%!error gammainc (1,1i)