octave: scripts/specfun/gammainc.m annotate

annotate scripts/specfun/gammainc.m @ 25216:b72972ab83f1 stable

test: Add bug IDs to failing %!xtest blocks. * gammainc.m, test.m, classdef.tst: Add bug ID to failing %!xtest blocks. * jit.tst: Add bug ID to failing %!xtest blocks. Add code to immediately exit from test if jit is not enabled. Need to do this because there is no %!xtestif keyword.

author	Rik <rik@octave.org>
date	Wed, 11 Apr 2018 11:52:49 -0700
parents	f98ef2b55641
children	b7e0ed1b2b32

rev	line source
25054 6652d3823428 maint: Update copyright dates in all source files. John W. Eaton <jwe@octave.org> parents: 24996 diff changeset	1 ## Copyright (C) 2016-2018 Marco Caliari
6652d3823428 maint: Update copyright dates in all source files. John W. Eaton <jwe@octave.org> parents: 24996 diff changeset	2 ## Copyright (C) 2016-2018 Nir Krakauer
24917 509e78c26e82 improved performance in gammainc Michele Ginesi <michele.ginesi@gmail.com> parents: 24905 diff changeset	3 ## Copyright (C) 2018 Stefan Schlögl
24923 40ab8be7d7ec Fixed style in specfun scripts Michele Ginesi <michele.ginesi@gmail.com> parents: 24918 diff changeset	4 ## Copyright (C) 2018 Michele Ginesi
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	5 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	6 ## This file is part of Octave.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	7 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	8 ## Octave is free software: you can redistribute it and/or modify it
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	9 ## under the terms of the GNU General Public License as published by
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	10 ## the Free Software Foundation; either version 3 of the License, or
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	11 ## (at your option) any later version.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	12 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	13 ## Octave is distributed in the hope that it will be useful, but
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	14 ## WITHOUT ANY WARRANTY; without even the implied warranty of
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	15 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	16 ## GNU General Public License for more details.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	17 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	18 ## You should have received a copy of the GNU General Public License
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	19 ## along with Octave; see the file COPYING. If not, see
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	20 ## <https://www.gnu.org/licenses/>.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	21
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	22 ## -- texinfo --
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	23 ## @deftypefn {} {} gammainc (@var{x}, @var{a})
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	24 ## @deftypefnx {} {} gammainc (@var{x}, @var{a}, @var{tail})
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	25 ## Compute the normalized incomplete gamma function.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	26 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	27 ## This is defined as
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	28 ## @tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	29 ## $$
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	30 ## \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt}
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	31 ## $$
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	32 ## @end tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	33 ## @ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	34 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	35 ## @example
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	36 ## @group
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	37 ## x
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	38 ## 1 /
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	39 ## gammainc (x, a) = --------- \| exp (-t) t^(a-1) dt
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	40 ## gamma (a) /
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	41 ## t=0
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	42 ## @end group
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	43 ## @end example
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	44 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	45 ## @end ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	46 ## with the limiting value of 1 as @var{x} approaches infinity.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	47 ## The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	48 ## @nospell{Stegun} (6.5.1).
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	49 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	50 ## If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	51 ## for each element of @var{x} and vice versa.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	52 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	53 ## If neither @var{x} nor @var{a} is scalar then the sizes of @var{x} and
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	54 ## @var{a} must agree, and @code{gammainc} is applied element-by-element.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	55 ## The elements of @var{a} must be non-negative.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	56 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	57 ## By default, @var{tail} is @qcode{"lower"} and the incomplete gamma function
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	58 ## integrated from 0 to @var{x} is computed. If @var{tail} is @qcode{"upper"}
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	59 ## then the complementary function integrated from @var{x} to infinity is
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	60 ## calculated.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	61 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	62 ## If @var{tail} is @qcode{"scaledlower"}, then the lower incomplete gamma
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	63 ## function is multiplied by
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	64 ## @tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	65 ## $\Gamma(a+1)\exp(x)x^{-a}$.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	66 ## @end tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	67 ## @ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	68 ## @math{gamma(a+1)*exp(x)/(x^a)}.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	69 ## @end ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	70 ## If @var{tail} is @qcode{"scaledupper"}, then the upper incomplete gamma
24939 5f4550b5d31b gammainc.m: Fix incorrect documentation. Rik <rik@octave.org> parents: 24934 diff changeset	71 ## function is multiplied by the same quantity.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	72 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	73 ## References:
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	74 ##
24987 a895967fd5a7 doc: grammarcheck manual (.txi) ahead of 4.4 release. Rik <rik@octave.org>* parents: 24939 diff changeset	75 ## @nospell{M. Abramowitz and I.A. Stegun},
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	76 ## @cite{Handbook of mathematical functions},
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	77 ## @nospell{Dover publications, Inc.}, 1972.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	78 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	79 ## @nospell{W. Gautschi},
24934 1db0b81efafe maint: strip trailing whitespace from source files Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	80 ## @cite{A computational procedure for incomplete gamma functions},
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	81 ## @nospell{ACM Trans. Math Software}, pp. 466--481, Vol 5, No. 4, 2012.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	82 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	83 ## @nospell{W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery},
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	84 ## @cite{Numerical Recipes in Fortran 77}, ch.@: 6.2, Vol 1, 1992.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	85 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	86 ## @seealso{gamma, gammaincinv, gammaln}
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	87 ## @end deftypefn
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	88
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	89 ## P(a,x) = gamma(a,x)/Gamma(a), upper
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	90 ## 1-P(a,x)=Q(a,x)=Gamma(a,x)/Gamma(a), lower
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	91
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	92 function y = gammainc (x, a, tail = "lower")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	93
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	94 if (nargin < 2 \|\| nargin > 3)
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	95 print_usage ();
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	96 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	97
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	98 [err, x, a] = common_size (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	99 if (err > 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	100 error ("gammainc: X and A must be of common size or scalars");
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	101 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	102
24934 1db0b81efafe maint: strip trailing whitespace from source files Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	103 if (iscomplex (x) \|\| iscomplex (a))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	104 error ("gammainc: all inputs must be real");
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	105 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	106
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	107 ## Remember original shape of data, but convert to column vector for calcs.
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	108 x_sz = size (x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	109 x = x(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	110 a = a(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	111
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	112 if (any (a < 0))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	113 error ("gammainc: A must be non-negative");
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	114 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	115
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	116 if (nargin == 3
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	117 && ! any (strcmpi (tail, {"lower","upper","scaledlower","scaledupper"})))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	118 error ("gammainc: invalid value for TAIL");
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	119 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	120 tail = tolower (tail);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	121
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	122 ## If any of the arguments is single then the output should be as well.
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	123 if (strcmp (class (x), "single") \|\| strcmp (class (a), "single"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	124 x = single (x);
24918 ea3edda05b66 Lentz, gammainc: added single precision Michele Ginesi <michele.ginesi@gmail.com> parents: 24917 diff changeset	125 a = single (a);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	126 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	127
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	128 ## Convert to floating point if necessary
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	129 if (isinteger (x))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	130 x = double (x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	131 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	132 if (isinteger (a))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	133 a = double (a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	134 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	135
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	136 ## Initialize output array
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	137 y = zeros (x_sz, class (x));
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	138
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	139 ## Different x, a combinations are handled by different subfunctions.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	140 todo = true (size (x)); # Track which elements need to be calculated.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	141
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	142 ## Case 0: x == Inf, a == Inf
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	143 idx = (x == Inf) & (a == Inf);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	144 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	145 y(idx) = NaN;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	146 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	147 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	148
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	149 ## Case 1: x == 0, a == 0.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	150 idx = (x == 0) & (a == 0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	151 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	152 y(idx) = gammainc_00 (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	153 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	154 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	155
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	156 ## Case 2: x == 0.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	157 idx = todo & (x == 0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	158 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	159 y(idx) = gammainc_x0 (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	160 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	161 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	162
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	163 ## Case 3: x = Inf
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	164 idx = todo & (x == Inf);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	165 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	166 y(idx) = gammainc_x_inf (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	167 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	168 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	169
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	170 ## Case 4: a = Inf
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	171 idx = todo & (a == Inf);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	172 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	173 y(idx) = gammainc_a_inf (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	174 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	175 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	176
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	177 ## Case 5: a == 0.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	178 idx = todo & (a == 0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	179 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	180 y(idx) = gammainc_a0 (x(idx), tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	181 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	182 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	183
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	184 ## Case 6: a == 1.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	185 idx = todo & (a == 1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	186 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	187 y(idx) = gammainc_a1 (x(idx), tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	188 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	189 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	190
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	191 ## Case 7: positive integer a; exp (x) and a! both under 1/eps.
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	192 idx = (todo
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	193 & (a == fix (a)) & (a > 1) & (a <= 18) & (x <= 36) & (abs (x) >= .1));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	194 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	195 y(idx) = gammainc_an (x(idx), a(idx), tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	196 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	197 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	198
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	199 ## For a < 2, x < 0, we increment a by 2 and use a recurrence formula after
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	200 ## the computations.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	201
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	202 flag_a_small = todo & (abs (a) > 0) & (abs (a) < 2) & (x < 0);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	203 a(flag_a_small) += 2;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	204
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	205 flag_s = (((x + 0.25 < a) \| (x < 0) \| (a < 5)) & (x > -20)) \| (abs (x) < 1);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	206
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	207 ## Case 8: x, a relatively small.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	208 idx = todo & flag_s;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	209 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	210 y(idx) = gammainc_s (x(idx), a(idx), tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	211 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	212 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	213
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	214 ## Case 9: x positive and large relative to a.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	215 idx = todo;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	216 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	217 y(idx) = gammainc_l (x(idx), a(idx), tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	218 todo(idx) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	219 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	220
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	221 if (any (flag_a_small))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	222 if (strcmp (tail, "lower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	223 y(flag_a_small) += D (x(flag_a_small), a(flag_a_small) - 1) + ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	224 D (x(flag_a_small), a(flag_a_small) - 2);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	225 elseif (strcmp (tail, "upper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	226 y(flag_a_small) -= D (x(flag_a_small), a(flag_a_small) - 1) + ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	227 D (x(flag_a_small), a(flag_a_small) - 2);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	228 elseif (strcmp (tail, "scaledlower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	229 y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	230 (a(flag_a_small) .* (a(flag_a_small) - 1)) + (x(flag_a_small) ./ ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	231 (a(flag_a_small) - 1)) + 1;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	232 elseif (strcmp (tail, "scaledupper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	233 y(flag_a_small) = y(flag_a_small) .* (x(flag_a_small) .^ 2) ./ ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	234 (a(flag_a_small) .* (a(flag_a_small) - 1)) - (x(flag_a_small) ./ ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	235 (a(flag_a_small) - 1)) - 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	236 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	237 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	238
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	239 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	240
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	241 ## Subfunctions to handle each case:
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	242
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	243 ## x == 0, a == 0.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	244 function y = gammainc_00 (tail)
24939 5f4550b5d31b gammainc.m: Fix incorrect documentation. Rik <rik@octave.org> parents: 24934 diff changeset	245 if (strcmp (tail, "upper") \|\| strcmp (tail, "scaledupper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	246 y = 0;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	247 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	248 y = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	249 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	250 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	251
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	252 ## x == 0.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	253 function y = gammainc_x0 (tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	254 if (strcmp (tail, "lower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	255 y = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	256 elseif (strcmp (tail, "upper") \|\| strcmp (tail, "scaledlower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	257 y = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	258 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	259 y = Inf;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	260 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	261 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	262
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	263 ## x == Inf.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	264 function y = gammainc_x_inf (tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	265 if (strcmp (tail, "lower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	266 y = 1;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	267 elseif (strcmp (tail, "upper") \|\| strcmp (tail, "scaledupper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	268 y = 0;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	269 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	270 y = Inf;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	271 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	272 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	273
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	274 ## a == Inf.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	275 function y = gammainc_a_inf (tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	276 if (strcmp (tail, "lower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	277 y = 0;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	278 elseif (strcmp (tail, "upper") \|\| strcmp (tail, "scaledlower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	279 y = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	280 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	281 y = Inf;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	282 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	283 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	284
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	285 ## a == 0.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	286 function y = gammainc_a0 (x, tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	287 if (strcmp (tail, "lower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	288 y = 1;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	289 elseif (strcmp (tail, "scaledlower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	290 y = exp (x);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	291 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	292 y = 0;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	293 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	294 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	295
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	296 ## a == 1.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	297 function y = gammainc_a1 (x, tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	298 if (strcmp (tail, "lower"))
25077 f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	299 if (abs (x) < 1/2)
f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	300 y = - expm1 (-x);
f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	301 else
f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	302 y = 1 - exp (-x);
f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	303 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	304 elseif (strcmp (tail, "upper"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	305 y = exp (-x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	306 elseif (strcmp (tail, "scaledlower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	307 if (abs (x) < 1/2)
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	308 y = expm1 (x) ./ x;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	309 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	310 y = (exp (x) - 1) ./ x;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	311 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	312 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	313 y = 1 ./ x;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	314 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	315 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	316
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	317 ## positive integer a; exp (x) and a! both under 1/eps
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	318 ## uses closed-form expressions for nonnegative integer a
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	319 ## -- http://mathworld.wolfram.com/IncompleteGammaFunction.html.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	320 function y = gammainc_an (x, a, tail)
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	321 y = t = ones (size (x), class (x));
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	322 i = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	323 while (any (a(:) > i))
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	324 jj = (a > i);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	325 t(jj) .*= (x(jj) / i);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	326 y(jj) += t(jj);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	327 i++;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	328 endwhile
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	329 if (strcmp (tail, "lower"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	330 y = 1 - exp (-x) .* y;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	331 elseif (strcmp (tail, "upper"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	332 y .*= exp (-x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	333 elseif (strcmp (tail, "scaledlower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	334 y = (1 - exp (-x) .* y) ./ D(x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	335 elseif (strcmp (tail, "scaledupper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	336 y .*= exp (-x) ./ D(x, a);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	337 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	338 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	339
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	340 ## x + 0.25 < a \| x < 0 \| x not real \| abs(x) < 1 \| a < 5.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	341 ## Numerical Recipes in Fortran 77 (6.2.5)
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	342 ## series
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	343 function y = gammainc_s (x, a, tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	344 if (strcmp (tail, "scaledlower") \|\| strcmp (tail, "scaledupper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	345 y = ones (size (x), class (x));
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	346 term = x ./ (a + 1);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	347 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	348 ## Of course it is possible to scale at the end, but some tests fail.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	349 ## And try gammainc (1,1000), it take 0 iterations if you scale now.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	350 y = D (x,a);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	351 term = y .* x ./ (a + 1);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	352 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	353 n = 1;
24925 24ae3461fb85 Fixed a typo introduced in gammainc.m Michele Ginesi <michele.ginesi@gmail.com> parents: 24923 diff changeset	354 while (any (abs (term(:)) > (abs (y(:)) * eps)))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	355 ## y can be zero from the beginning (gammainc (1,1000))
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	356 jj = abs (term) > abs (y) * eps;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	357 n += 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	358 y(jj) += term(jj);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	359 term(jj) .*= x(jj) ./ (a(jj) + n);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	360 endwhile
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	361 if (strcmp (tail, "upper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	362 y = 1 - y;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	363 elseif (strcmp (tail, "scaledupper"))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	364 y = 1 ./ D (x,a) - y;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	365 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	366 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	367
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	368 ## x positive and large relative to a
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	369 ## NRF77 (6.2.7)
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	370 ## Gamma (a,x)/Gamma (a)
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	371 ## Lentz's algorithm
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	372 ## __gammainc__ in libinterp/corefcn/__gammainc__.cc
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	373 function y = gammainc_l (x, a, tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	374 y = __gammainc__ (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	375 if (strcmp (tail, "lower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	376 y = 1 - y .* D (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	377 elseif (strcmp (tail, "upper"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	378 y .*= D (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	379 elseif (strcmp (tail, "scaledlower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	380 y = 1 ./ D (x, a) - y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	381 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	382 endfunction
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	383
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	384 ## Compute exp(-x)*x^a/Gamma(a+1) in a stable way for x and a large.
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	385 ##
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	386 ## L. Knusel, Computation of the Chi-square and Poisson distribution,
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	387 ## SIAM J. Sci. Stat. Comput., 7(3), 1986
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	388 ## which quotes Section 5, Abramowitz&Stegun 6.1.40, 6.1.41.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	389 function y = D (x, a)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	390 athresh = 10; # FIXME: can this be better tuned?
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	391 y = zeros (size (x), class (x));
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	392
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	393 todo = true (size (x));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	394 todo(x == 0) = false;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	395
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	396 ii = todo & (x > 0) & (a > athresh) & (a >= x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	397 if (any (ii))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	398 lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	399 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	400 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	401 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	402 3617 ./ (122400 * a(ii) .^ 15) + ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	403 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	404 lns = log1p ((a(ii) - x(ii)) ./ x(ii));
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	405 y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	406 todo(ii) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	407 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	408
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	409 ii = todo & (x > 0) & (a > athresh) & (a < x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	410 if (any (ii))
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	411 lnGa = log (2 * pi * a(ii)) / 2 + 1 ./ (12 * a(ii)) - ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	412 1 ./ (360 * a(ii) .^ 3) + 1 ./ (1260 * a(ii) .^ 5) - ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	413 1 ./ (1680 * a(ii) .^ 7) + 1 ./ (1188 * a(ii) .^ 9)- ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	414 691 ./ (87360 * a(ii) .^ 11) + 1 ./ (156 * a(ii) .^ 13) - ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	415 3617 ./ (122400 * a(ii) .^ 15) + ...
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	416 43867 ./ (244188 * a(ii) .^ 17) - 174611 ./ (125400 * a(ii) .^ 19);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	417 lns = -log1p ((x(ii) - a(ii)) ./ a(ii));
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	418 y(ii) = exp ((a(ii) - x(ii)) - a(ii) .* lns - lnGa);
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	419 todo(ii) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	420 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	421
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	422 ii = todo & ((x <= 0) \| (a <= athresh));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	423 if (any (ii)) # standard formula for a not so large.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	424 y(ii) = exp (a(ii) .* log (x(ii)) - x(ii) - gammaln (a(ii) + 1));
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	425 todo(ii) = false;
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	426 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	427
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	428 ii = (x < 0) & (a == fix (a));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	429 if (any(ii)) # remove spurious imaginary part.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	430 y(ii) = real (y(ii));
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	431 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	432
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	433 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	434
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	435
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	436 ## Test: case 1,2,5
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	437 %!assert (gammainc ([0, 0, 1], [0, 1, 0]), [1, 0, 1])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	438 %!assert (gammainc ([0, 0, 1], [0, 1, 0], "upper"), [0, 1, 0])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	439 %!assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledlower"), [1, 1, exp(1)])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	440 %!assert (gammainc ([0, 0, 1], [0, 1, 0], "scaledupper"), [0, Inf, 0])
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	441
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	442 ## Test: case 3,4
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	443 %!assert (gammainc ([2, Inf], [Inf, 2]), [0, 1])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	444 %!assert (gammainc ([2, Inf], [Inf, 2], "upper"), [1, 0])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	445 %!assert (gammainc ([2, Inf], [Inf, 2], "scaledlower"), [1, Inf])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	446 %!assert (gammainc ([2, Inf], [Inf, 2], "scaledupper"), [Inf, 0])
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	447
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	448 ## Test: case 5
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	449 ## Matlab fails for this test
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	450 %!assert (gammainc (-100,1,"upper"), exp (100), -eps)
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	451
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	452 ## Test: case 6
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	453 %!assert (gammainc ([1, 2, 3], 1), 1 - exp (-[1, 2, 3]))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	454 %!assert (gammainc ([1, 2, 3], 1, "upper"), exp (- [1, 2, 3]))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	455 %!assert (gammainc ([1, 2, 3], 1, "scaledlower"), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	456 %! (exp ([1, 2, 3]) - 1) ./ [1, 2, 3])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	457 %!assert (gammainc ([1, 2, 3], 1, "scaledupper"), 1 ./ [1, 2, 3])
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	458
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	459 ## Test: case 7
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	460 %!assert (gammainc (2, 2, "lower"), 0.593994150290162, -2e-15)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	461 %!assert (gammainc (2, 2, "upper"), 0.406005849709838, -2e-15)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	462 %!assert (gammainc (2, 2, "scaledlower"), 2.194528049465325, -2e-15)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	463 %!assert (gammainc (2, 2, "scaledupper"), 1.500000000000000, -2e-15)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	464 %!assert (gammainc ([3 2 36],[2 3 18], "upper"), ...
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24987 diff changeset	465 %! [4/exp(3) 5*exp(-2) (4369755579265807723 / 2977975)/exp(36)], -eps)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	466 %!assert (gammainc (10, 10), 1 - (5719087 / 567) * exp (-10), -eps)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	467 %!assert (gammainc (10, 10, "upper"), (5719087 / 567) * exp (-10), -eps)
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	468
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	469 ## Test: case 8
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	470 %!assert (gammainc (-10, 10), 3.112658265341493126871617e7, -2*eps)
24934 1db0b81efafe maint: strip trailing whitespace from source files Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	471 ## Matlab fails this next one%! %!
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	472 %!assert (isreal (gammainc (-10, 10)), true)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	473 %!assert (gammainc (-10, 10.1, "upper"), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	474 %! -2.9582761911890713293e7-1i * 9.612022339061679758e6, -30*eps)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	475 %!assert (gammainc (-10, 10, "upper"), -3.112658165341493126871616e7, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	476 %! -2*eps)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	477 %!assert (gammainc (-10, 10, "scaledlower"), 0.5128019364747265, -1e-14);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	478 %!assert (gammainc (-10, 10, "scaledupper"), -0.5128019200000000, -1e-14);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	479 %!assert (gammainc (200, 201, "upper"), 0.518794309678684497, -2 * eps);
24934 1db0b81efafe maint: strip trailing whitespace from source files Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	480 %!assert (gammainc (200, 201, "scaledupper"),
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	481 %! 18.4904360746560462660798514, -eps)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	482 ## Here we are very good (no D (x,a)) involved
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	483 %!assert (gammainc(1000, 1000.5, "scaledlower"), 39.48467539583672271, -2*eps)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	484 %!assert (gammainc (709, 1000, "upper"), 0.99999999999999999999999954358, -eps)
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	485
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	486 ## Test: case 9
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	487 %!test <47800>
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	488 %! assert (gammainc (60, 6, "upper"), 6.18022358081160257327264261e-20,
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	489 %! -10*eps);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	490 ## Matlab is better here than Octave
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	491 %!assert (gammainc (751, 750, "upper"), 0.4805914320558831327179457887, -12*eps)
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24987 diff changeset	492 %!assert (gammainc (200, 200, "upper"), 0.49059658199276367497217454, -5*eps)
f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24987 diff changeset	493 %!assert (gammainc (200, 200), 0.509403418007236325027825459574527043, -5*eps)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	494 %!assert (gammainc (200, 200, "scaledupper"), 17.3984438553791505135122900,
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24987 diff changeset	495 %! -2*eps)
f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24987 diff changeset	496 %!assert (gammainc (200, 200, "scaledlower"), 18.065406676779221643065, -7*eps)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	497 %!assert (gammainc (201, 200, "upper"), 0.46249244908276709524913736667, -7*eps)
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	498
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	499 ## Test small argument
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	500 %!assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.1), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	501 %! [0.33239840504050, 0.20972940370977, 0.10511370061022, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	502 %! 0.041846517936723], 1e-13);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	503
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	504 %!assert (gammainc ([1e-05, 1e-07,1e-10,1e-14], 0.2), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	505 %! [0.10891226058559, 0.043358823442178, 0.010891244210402, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	506 %! 0.0017261458806785], 1e-13);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	507
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	508 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	509 %!assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 0.9), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	510 %! [0.016401189184068, 0.0020735998660840, 0.000032879756964708, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	511 %! 8.2590606569241e-9, 2.6117443021738e-13], -1e-12);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	512
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	513 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	514 %!assert (gammainc ([1e-02, 1e-03, 1e-5, 1e-9, 1e-14], 2), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	515 %! [0.0000496679133402659, 4.99666791633340e-7, 4.99996666679167e-11, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	516 %! 4.99999999666667e-19, 4.99999999999997e-29], -1e-12);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	517
25077 f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	518 %!test <*53543>
f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	519 %! y_exp = 9.995001666250085e-04;
f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	520 %! assert (gammainc (1/1000, 1), y_exp, -eps);
f98ef2b55641 gammainc.m: Fixed issue with a=1, and x is small (bug #53543). Michele Ginesi <michele.ginesi@gmail.com> parents: 25054 diff changeset	521
25216 b72972ab83f1 test: Add bug IDs to failing %!xtest blocks. Rik <rik@octave.org> parents: 25077 diff changeset	522 %!xtest <53612>
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	523 %! assert (gammainc (-20, 1.1, "upper"), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	524 %! 6.50986687074979e8 + 2.11518396291149e8*i, -1e-13);
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	525
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	526 ## Test conservation of the class (five tests for each subroutine).
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	527 %!assert (class (gammainc (0, 1)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	528 %!assert (class (gammainc (single (0), 1)) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	529 %!assert (class (gammainc (int8 (0), 1)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	530 %!assert (class (gammainc (0, single (1))) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	531 %!assert (class (gammainc (0, int8 (1))) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	532 %!assert (class (gammainc (1, 0)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	533 %!assert (class (gammainc (single (1), 0)) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	534 %!assert (class (gammainc (int8 (1), 0)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	535 %!assert (class (gammainc (1, single (0))) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	536 %!assert (class (gammainc (1, int8 (0))) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	537 %!assert (class (gammainc (1, 1)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	538 %!assert (class (gammainc (single (1), 1)) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	539 %!assert (class (gammainc (int8 (1), 1)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	540 %!assert (class (gammainc (1, single (1))) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	541 %!assert (class (gammainc (1, int8(1))) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	542 %!assert (class (gammainc (1, 2)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	543 %!assert (class (gammainc (single (1), 2)) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	544 %!assert (class (gammainc (int8 (1), 2)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	545 %!assert (class (gammainc (1, single (2))) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	546 %!assert (class (gammainc (1, int8 (2))) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	547 %!assert (class (gammainc (-1, 0.5)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	548 %!assert (class (gammainc (single (-1), 0.5)) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	549 %!assert (class (gammainc (int8 (-1), 0.5)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	550 %!assert (class (gammainc (-1, single (0.5))) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	551 %!assert (class (gammainc (-1, int8 (0.5))) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	552 %!assert (class (gammainc (1, 0.5)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	553 %!assert (class (gammainc (single (1), 0.5)) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	554 %!assert (class (gammainc (int8 (1), 0.5)) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	555 %!assert (class (gammainc (1, single (0.5))) == "single")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	556 %!assert (class (gammainc (1, int8 (0.5))) == "double")
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	557
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	558 ## Test input validation
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	559 %!error gammainc ()
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	560 %!error gammainc (1)
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	561 %!error gammainc (1,2,3,4)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	562 %!error <must be of common size or scalars> gammainc ([0, 0],[0; 0])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	563 %!error <must be of common size or scalars> gammainc ([1 2 3], [1 2])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	564 %!error <all inputs must be real> gammainc (2+i, 1)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	565 %!error <all inputs must be real> gammainc (1, 2+i)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	566 %!error <A must be non-negative> gammainc (1, [0, -1, 1])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	567 %!error <A must be non-negative>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	568 %! a = ones (2,2,2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	569 %! a(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	570 %! gammainc (1, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	571 %!error <invalid value for TAIL> gammainc (1,2, "foobar")

Mercurial > octave

annotate scripts/specfun/gammainc.m @ 25216:b72972ab83f1 stable