octave: scripts/specfun/gammainc.m annotate

annotate scripts/specfun/gammainc.m @ 31230:6646f2b5a3d1

lcm.m: Emit warnings when results not exact (Bug #32924)

author	Arun Giridhar <arungiridhar@gmail.com>
date	Sat, 17 Sep 2022 04:22:38 -0400
parents	5d3faba0342e
children	c8ad083a5802

rev	line source
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	1 ########################################################################
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	2 ##
30564 796f54d4ddbf update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 29494 diff changeset	3 ## Copyright (C) 2016-2022 The Octave Project Developers
27918 b442ec6dda5c use centralized file for copyright info for individual contributors John W. Eaton <jwe@octave.org> parents: 27800 diff changeset	4 ##
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	6 ## distribution or <https://octave.org/copyright/>.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	7 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	8 ## 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	9 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	10 ## 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	11 ## under the terms of the GNU General Public License as published by
28944 d39b09b4c5db maint: Use identical Copyright text in scripts/. Rik <rik@octave.org> parents: 28912 diff changeset	12 ## the Free Software Foundation, either version 3 of the License, or
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	13 ## (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	14 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	15 ## 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	16 ## 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	17 ## 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	18 ## 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	19 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	20 ## 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	21 ## 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	22 ## <https://www.gnu.org/licenses/>.
27923 bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	23 ##
bd51beb6205e update formatting of copyright notices John W. Eaton <jwe@octave.org> parents: 27919 diff changeset	24 ########################################################################
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	25
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	26 ## -- texinfo --
30875 5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files. Rik <rik@octave.org> parents: 30564 diff changeset	27 ## @deftypefn {} {@var{y} =} gammainc (@var{x}, @var{a})
5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files. Rik <rik@octave.org> parents: 30564 diff changeset	28 ## @deftypefnx {} {@var{y} =} gammainc (@var{x}, @var{a}, @var{tail})
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	29 ## 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	30 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	31 ## 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	32 ## @tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	33 ## $$
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	34 ## \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	35 ## $$
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	36 ## @end tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	37 ## @ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	38 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	39 ## @example
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	40 ## @group
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	41 ## x
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	42 ## 1 /
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	43 ## 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	44 ## gamma (a) /
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	45 ## t=0
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	46 ## @end group
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	47 ## @end example
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	48 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	49 ## @end ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	50 ## 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	51 ## 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	52 ## @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	53 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	54 ## 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	55 ## 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	56 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	57 ## 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	58 ## @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	59 ## 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	60 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	61 ## 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	62 ## 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	63 ## 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	64 ## calculated.
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	65 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	66 ## 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	67 ## 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	68 ## @tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	69 ## $\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	70 ## @end tex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	71 ## @ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	72 ## @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	73 ## @end ifnottex
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	74 ## 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	75 ## 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	76 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	77 ## References:
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	78 ##
24987 a895967fd5a7 doc: grammarcheck manual (.txi) ahead of 4.4 release. Rik <rik@octave.org>* parents: 24939 diff changeset	79 ## @nospell{M. Abramowitz and I.A. Stegun},
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	80 ## @cite{Handbook of mathematical functions},
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	81 ## @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	82 ##
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	83 ## @nospell{W. Gautschi},
24934 1db0b81efafe maint: strip trailing whitespace from source files Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	84 ## @cite{A computational procedure for incomplete gamma functions},
27800 5a6a19a4e3da doc: Use Texinfo non-sentence ending periods in citations. Rik <rik@octave.org> parents: 26946 diff changeset	85 ## @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	86 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	87 ## @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	88 ## @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	89 ##
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	90 ## @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	91 ## @end deftypefn
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	92
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	93 ## 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	94 ## 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	95
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	96 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	97
28789 28de41192f3c Eliminate unneeded verification of nargin, nargout in m-files. Rik <rik@octave.org> parents: 27923 diff changeset	98 if (nargin < 2)
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	99 print_usage ();
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	100 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	101
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	102 [err, x, a] = common_size (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	103 if (err > 0)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	104 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	105 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	106
24934 1db0b81efafe maint: strip trailing whitespace from source files Mike Miller <mtmiller@octave.org> parents: 24927 diff changeset	107 if (iscomplex (x) \|\| iscomplex (a))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	108 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	109 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	110
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	111 ## 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	112 x_sz = size (x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	113 x = x(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	114 a = a(:);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	115
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	116 if (any (a < 0))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	117 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	118 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	119
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	120 if (nargin == 3
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	121 && ! any (strcmpi (tail, {"lower","upper","scaledlower","scaledupper"})))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	122 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	123 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	124 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	125
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	126 ## 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	127 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	128 x = single (x);
24918 ea3edda05b66 Lentz, gammainc: added single precision Michele Ginesi <michele.ginesi@gmail.com> parents: 24917 diff changeset	129 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	130 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	131
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	132 ## Convert to floating point if necessary
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	133 if (isinteger (x))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	134 x = double (x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	135 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	136 if (isinteger (a))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	137 a = double (a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	138 endif
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	139
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	140 ## Initialize output array
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	141 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	142
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	143 ## 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	144 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	145
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	146 ## Case 0: x == Inf, a == Inf
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	147 idx = (x == Inf) & (a == Inf);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	148 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	149 y(idx) = NaN;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	150 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	151 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	152
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	153 ## Case 1: x == 0, a == 0.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	154 idx = (x == 0) & (a == 0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	155 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	156 y(idx) = gammainc_00 (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	157 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	158 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	159
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	160 ## Case 2: x == 0.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	161 idx = todo & (x == 0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	162 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	163 y(idx) = gammainc_x0 (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	164 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	165 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	166
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	167 ## Case 3: x = Inf
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	168 idx = todo & (x == Inf);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	169 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	170 y(idx) = gammainc_x_inf (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	171 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	172 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	173
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	174 ## Case 4: a = Inf
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	175 idx = todo & (a == Inf);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	176 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	177 y(idx) = gammainc_a_inf (tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	178 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	179 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	180
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	181 ## Case 5: a == 0.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	182 idx = todo & (a == 0);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	183 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	184 y(idx) = gammainc_a0 (x(idx), tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	185 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	186 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	187
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	188 ## Case 6: a == 1.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	189 idx = todo & (a == 1);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	190 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	191 y(idx) = gammainc_a1 (x(idx), tail);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	192 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	193 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	194
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	195 ## 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	196 idx = (todo
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	197 & (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	198 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	199 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	200 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	201 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	202
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	203 ## 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	204 ## the computations.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	205
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	206 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	207 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	208
25825 b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	209 flag_s = (((x + 0.25 < a) \| (x < 0)) & (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	210
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	211 ## Case 8: x, a relatively small.
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	212 idx = todo & flag_s;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	213 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	214 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	215 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	216 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	217
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	218 ## 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	219 idx = todo;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	220 if (any (idx))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	221 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	222 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	223 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	224
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	225 if (any (flag_a_small))
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	226 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	227 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	228 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	229 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	230 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	231 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	232 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	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;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	236 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	237 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	238 (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	239 (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	240 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	241 endif
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 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	244
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	245 ## 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	246
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	247 ## 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	248 function y = gammainc_00 (tail)
24939 5f4550b5d31b gammainc.m: Fix incorrect documentation. Rik <rik@octave.org> parents: 24934 diff changeset	249 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	250 y = 0;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	251 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	252 y = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	253 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	254 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	255
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	256 ## x == 0.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	257 function y = gammainc_x0 (tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	258 if (strcmp (tail, "lower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	259 y = 0;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	260 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	261 y = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	262 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	263 y = Inf;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	264 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	265 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	266
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	267 ## x == Inf.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	268 function y = gammainc_x_inf (tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	269 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	270 y = 1;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	271 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	272 y = 0;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	273 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	274 y = Inf;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	275 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	276 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	277
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	278 ## a == Inf.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	279 function y = gammainc_a_inf (tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	280 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	281 y = 0;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	282 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	283 y = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	284 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	285 y = Inf;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	286 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	287 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	288
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	289 ## a == 0.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	290 function y = gammainc_a0 (x, tail)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	291 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	292 y = 1;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	293 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	294 y = exp (x);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	295 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	296 y = 0;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	297 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	298 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	299
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	300 ## a == 1.
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	301 function y = gammainc_a1 (x, tail)
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	302
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	303 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	304 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	305 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	306 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	307 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	308 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	309 elseif (strcmp (tail, "upper"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	310 y = exp (-x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	311 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	312 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	313 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	314 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	315 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	316 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	317 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	318 y = 1 ./ x;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	319 endif
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	320
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	321 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	322
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	323 ## 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	324 ## 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	325 ## -- 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	326 function y = gammainc_an (x, a, tail)
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	327
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	328 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	329 i = 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	330 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	331 jj = (a > i);
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	332 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	333 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	334 i++;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	335 endwhile
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	336 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	337 y = 1 - exp (-x) .* y;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	338 elseif (strcmp (tail, "upper"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	339 y .*= exp (-x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	340 elseif (strcmp (tail, "scaledlower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	341 y = (1 - exp (-x) .* y) ./ D(x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	342 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	343 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	344 endif
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	345
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	346 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	347
25825 b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	348 ## x + 0.25 < a \| x < 0 \| 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	349 ## 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	350 ## series
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	351 function y = gammainc_s (x, a, tail)
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	352
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	353 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	354 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	355 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	356 else
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	357 ## 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	358 ## 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	359 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	360 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	361 endif
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	362 n = 1;
24925 24ae3461fb85 Fixed a typo introduced in gammainc.m Michele Ginesi <michele.ginesi@gmail.com> parents: 24923 diff changeset	363 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	364 ## 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	365 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	366 n += 1;
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	367 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	368 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	369 endwhile
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	370 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	371 y = 1 - y;
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	372 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	373 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	374 endif
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	375
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	376 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	377
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	378 ## 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	379 ## 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	380 ## 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	381 ## Lentz's algorithm
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	382 ## __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	383 function y = gammainc_l (x, a, tail)
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	384
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	385 y = __gammainc__ (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	386 if (strcmp (tail, "lower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	387 y = 1 - y .* D (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	388 elseif (strcmp (tail, "upper"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	389 y .*= D (x, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	390 elseif (strcmp (tail, "scaledlower"))
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	391 y = 1 ./ D (x, a) - y;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	392 endif
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	393
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	394 endfunction
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	395
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	396 ## 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	397 ##
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	398 ## 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	399 ## SIAM J. Sci. Stat. Comput., 7(3), 1986
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	400 ## 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	401 function y = D (x, a)
28945 6e460773bdda maint: Use newlines after "function" and before "endfunction" for clarity. Rik <rik@octave.org> parents: 28944 diff changeset	402
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	403 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	404 y = zeros (size (x), class (x));
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	405
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	406 todo = true (size (x));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	407 todo(x == 0) = false;
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 ((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	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) & (a < x);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	423 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	424 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	425 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	426 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	427 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	428 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	429 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	430 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	431 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	432 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	433 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	434
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	435 ii = todo & ((x <= 0) \| (a <= athresh));
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	436 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	437 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	438 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	439 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	440
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	441 ii = (x < 0) & (a == fix (a));
28912 0de38a6ef693 maint: Use Octave convention of space after function name in scripts dir. Rik <rik@octave.org> parents: 28907 diff changeset	442 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	443 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	444 endif
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	445
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	446 endfunction
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	447
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	448
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	449 ## Test: case 1,2,5
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	450 %!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	451 %!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	452 %!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	453 %!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	454
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	455 ## Test: case 3,4
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	456 %!assert (gammainc ([2, Inf], [Inf, 2]), [0, 1])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	457 %!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	458 %!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	459 %!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	460
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	461 ## Test: case 5
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	462 ## Matlab fails for this test
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	463 %!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	464
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	465 ## Test: case 6
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	466 %!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	467 %!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	468 %!assert (gammainc ([1, 2, 3], 1, "scaledlower"), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	469 %! (exp ([1, 2, 3]) - 1) ./ [1, 2, 3])
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	470 %!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	471
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	472 ## Test: case 7
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	473 %!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	474 %!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	475 %!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	476 %!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	477 %!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	478 %! [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	479 %!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	480 %!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	481
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	482 ## Test: case 8
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	483 %!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	484 ## Matlab fails this next one%! %!
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	485 %!assert (isreal (gammainc (-10, 10)), true)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	486 %!assert (gammainc (-10, 10.1, "upper"), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	487 %! -2.9582761911890713293e7-1i * 9.612022339061679758e6, -30*eps)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	488 %!assert (gammainc (-10, 10, "upper"), -3.112658165341493126871616e7, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	489 %! -2*eps)
28907 11f1207111c5 maint: Don't use semicolon at end of single-line BIST tests. Rik <rik@octave.org> parents: 28896 diff changeset	490 %!assert (gammainc (-10, 10, "scaledlower"), 0.5128019364747265, -1e-14)
11f1207111c5 maint: Don't use semicolon at end of single-line BIST tests. Rik <rik@octave.org> parents: 28896 diff changeset	491 %!assert (gammainc (-10, 10, "scaledupper"), -0.5128019200000000, -1e-14)
11f1207111c5 maint: Don't use semicolon at end of single-line BIST tests. Rik <rik@octave.org> parents: 28896 diff changeset	492 %!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	493 %!assert (gammainc (200, 201, "scaledupper"),
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	494 %! 18.4904360746560462660798514, -eps)
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	495 ## Here we are very good (no D (x,a)) involved
26680 eb1d4a619260 gammainc.m: relax tolerances of BIST tests Mike Miller <mtmiller@octave.org> parents: 26376 diff changeset	496 %!assert (gammainc (1000, 1000.5, "scaledlower"), 39.48467539583672271, -2*eps)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	497 %!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	498
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	499 ## Test: case 9
26946 04e5cb5e2cb3 update bug status in tests John W. Eaton <jwe@octave.org> parents: 26680 diff changeset	500 %!test <*47800>
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	501 %! assert (gammainc (60, 6, "upper"), 6.18022358081160257327264261e-20,
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	502 %! -10*eps);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	503 ## Matlab is better here than Octave
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	504 %!assert (gammainc (751, 750, "upper"), 0.4805914320558831327179457887, -12*eps)
26680 eb1d4a619260 gammainc.m: relax tolerances of BIST tests Mike Miller <mtmiller@octave.org> parents: 26376 diff changeset	505 %!assert (gammainc (200, 200, "upper"), 0.49059658199276367497217454, -6*eps)
24996 f80e68529bbf test: relax tolerances on specfun tests to pass on i386 (bug #53437) Mike Miller <mtmiller@octave.org> parents: 24987 diff changeset	506 %!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	507 %!assert (gammainc (200, 200, "scaledupper"), 17.3984438553791505135122900,
26680 eb1d4a619260 gammainc.m: relax tolerances of BIST tests Mike Miller <mtmiller@octave.org> parents: 26376 diff changeset	508 %! -3*eps)
eb1d4a619260 gammainc.m: relax tolerances of BIST tests Mike Miller <mtmiller@octave.org> parents: 26376 diff changeset	509 %!assert (gammainc (200, 200, "scaledlower"), 18.065406676779221643065, -8*eps)
25825 b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	510 %!assert (gammainc (201, 200, "upper"), 0.46249244908276709524913736667,
b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	511 %! -7*eps)
b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	512 %!assert <*54550> (gammainc (77, 2), 1)
b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	513
b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	514 %!assert (gammainc (77, 2, "upper"), 0, -eps)
b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	515 %!assert (gammainc (1000, 3.1), 1)
b7e0ed1b2b32 Fix inaccurate result of gammainc for x large (bug #54550). Marco Caliari <marco.caliari@univr.it> parents: 25216 diff changeset	516 %!assert (gammainc (1000, 3.1, "upper"), 0)
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	517
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	518 ## Test small argument
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	519 %!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	520 %! [0.33239840504050, 0.20972940370977, 0.10511370061022, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	521 %! 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	522
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	523 %!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	524 %! [0.10891226058559, 0.043358823442178, 0.010891244210402, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	525 %! 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	526
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	527 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	528 %!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	529 %! [0.016401189184068, 0.0020735998660840, 0.000032879756964708, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	530 %! 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	531
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	532 %!test
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	533 %!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	534 %! [0.0000496679133402659, 4.99666791633340e-7, 4.99996666679167e-11, ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	535 %! 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	536
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	537 %!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	538 %! 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	539 %! 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	540
29494 76269aa97398 if bug number is supplied, prefer '%!test' over '%!xtest' John W. Eaton <jwe@octave.org> parents: 29359 diff changeset	541 %!test <53612>
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	542 %! assert (gammainc (-20, 1.1, "upper"), ...
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	543 %! 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	544
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	545 ## 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	546 %!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	547 %!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	548 %!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	549 %!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	550 %!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	551 %!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	552 %!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	553 %!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	554 %!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	555 %!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	556 %!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	557 %!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	558 %!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	559 %!assert (class (gammainc (1, single (1))) == "single")
26680 eb1d4a619260 gammainc.m: relax tolerances of BIST tests Mike Miller <mtmiller@octave.org> parents: 26376 diff changeset	560 %!assert (class (gammainc (1, int8 (1))) == "double")
24904 740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	561 %!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	562 %!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	563 %!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	564 %!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	565 %!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	566 %!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	567 %!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	568 %!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	569 %!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	570 %!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	571 %!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	572 %!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	573 %!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	574 %!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	575 %!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	576
740df3fce3fb New strategy to compute the incomplete gamma function (see bug #47800). Michele Ginesi <michele.ginesi@gmail.com> parents: diff changeset	577 ## Test input validation
28896 90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin. Rik <rik@octave.org> parents: 28789 diff changeset	578 %!error <Invalid call> gammainc ()
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin. Rik <rik@octave.org> parents: 28789 diff changeset	579 %!error <Invalid call> gammainc (1)
24927 c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	580 %!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	581 %!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	582 %!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	583 %!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	584 %!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	585 %!error <A must be non-negative>
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	586 %! a = ones (2,2,2);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	587 %! a(1,1,2) = -1;
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	588 %! gammainc (1, a);
c280560d9c96 Overhaul special functions modified by GSOC2018 project. Rik <rik@octave.org> parents: 24925 diff changeset	589 %!error <invalid value for TAIL> gammainc (1,2, "foobar")

Mercurial > octave

annotate scripts/specfun/gammainc.m @ 31230:6646f2b5a3d1