GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
gammainc.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1997-2015 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "lo-specfun.h"
28 
29 #include "defun.h"
30 #include "error.h"
31 #include "gripes.h"
32 #include "oct-obj.h"
33 #include "utils.h"
34 
35 DEFUN (gammainc, args, ,
36  "-*- texinfo -*-\n\
37 @deftypefn {Mapping Function} {} gammainc (@var{x}, @var{a})\n\
38 @deftypefnx {Mapping Function} {} gammainc (@var{x}, @var{a}, \"lower\")\n\
39 @deftypefnx {Mapping Function} {} gammainc (@var{x}, @var{a}, \"upper\")\n\
40 Compute the normalized incomplete gamma function.\n\
41 \n\
42 This is defined as\n\
43 @tex\n\
44 $$\n\
45  \\gamma (x, a) = {1 \\over {\\Gamma (a)}}\\displaystyle{\\int_0^x t^{a-1} e^{-t} dt}\n\
46 $$\n\
47 @end tex\n\
48 @ifnottex\n\
49 \n\
50 @example\n\
51 @group\n\
52  x\n\
53  1 /\n\
54 gammainc (x, a) = --------- | exp (-t) t^(a-1) dt\n\
55  gamma (a) /\n\
56  t=0\n\
57 @end group\n\
58 @end example\n\
59 \n\
60 @end ifnottex\n\
61 with the limiting value of 1 as @var{x} approaches infinity.\n\
62 The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and\n\
63 @nospell{Stegun} (6.5.1).\n\
64 \n\
65 If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned\n\
66 for each element of @var{x} and vice versa.\n\
67 \n\
68 If neither @var{x} nor @var{a} is scalar, the sizes of @var{x} and\n\
69 @var{a} must agree, and @code{gammainc} is applied element-by-element.\n\
70 \n\
71 By default the incomplete gamma function integrated from 0 to @var{x} is\n\
72 computed. If @qcode{\"upper\"} is given then the complementary function\n\
73 integrated from @var{x} to infinity is calculated. It should be noted that\n\
74 \n\
75 @example\n\
76 gammainc (@var{x}, @var{a}) @equiv{} 1 - gammainc (@var{x}, @var{a}, \"upper\")\n\
77 @end example\n\
78 @seealso{gamma, gammaln}\n\
79 @end deftypefn")
80 {
81  octave_value retval;
82  bool lower = true;
83 
84  int nargin = args.length ();
85 
86  if (nargin == 3)
87  {
88  if (args(2).is_string ())
89  {
90  std::string s = args(2).string_value ();
91  std::transform (s.begin (), s.end (), s.begin (), tolower);
92  if (s == "upper")
93  lower = false;
94  else if (s != "lower")
95  error ("gammainc: third argument must be \"lower\" or \"upper\"");
96  }
97  else
98  error ("gammainc: third argument must be \"lower\" or \"upper\"");
99 
100  }
101 
102  if (!error_state && nargin >= 2 && nargin <= 3)
103  {
104  octave_value x_arg = args(0);
105  octave_value a_arg = args(1);
106 
107  // FIXME: Can we make a template version of the duplicated code below
108  if (x_arg.is_single_type () || a_arg.is_single_type ())
109  {
110  if (x_arg.is_scalar_type ())
111  {
112  float x = x_arg.float_value ();
113 
114  if (! error_state)
115  {
116  if (a_arg.is_scalar_type ())
117  {
118  float a = a_arg.float_value ();
119 
120  if (! error_state)
121  retval = lower ? gammainc (x, a)
122  : 1.0f - gammainc (x, a);
123  }
124  else
125  {
126  FloatNDArray a = a_arg.float_array_value ();
127 
128  if (! error_state)
129  retval = lower ? gammainc (x, a)
130  : 1.0f - gammainc (x, a);
131  }
132  }
133  }
134  else
135  {
136  FloatNDArray x = x_arg.float_array_value ();
137 
138  if (! error_state)
139  {
140  if (a_arg.is_scalar_type ())
141  {
142  float a = a_arg.float_value ();
143 
144  if (! error_state)
145  retval = lower ? gammainc (x, a)
146  : 1.0f - gammainc (x, a);
147  }
148  else
149  {
150  FloatNDArray a = a_arg.float_array_value ();
151 
152  if (! error_state)
153  retval = lower ? gammainc (x, a)
154  : 1.0f - gammainc (x, a);
155  }
156  }
157  }
158  }
159  else
160  {
161  if (x_arg.is_scalar_type ())
162  {
163  double x = x_arg.double_value ();
164 
165  if (! error_state)
166  {
167  if (a_arg.is_scalar_type ())
168  {
169  double a = a_arg.double_value ();
170 
171  if (! error_state)
172  retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
173  }
174  else
175  {
176  NDArray a = a_arg.array_value ();
177 
178  if (! error_state)
179  retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
180  }
181  }
182  }
183  else
184  {
185  NDArray x = x_arg.array_value ();
186 
187  if (! error_state)
188  {
189  if (a_arg.is_scalar_type ())
190  {
191  double a = a_arg.double_value ();
192 
193  if (! error_state)
194  retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
195  }
196  else
197  {
198  NDArray a = a_arg.array_value ();
199 
200  if (! error_state)
201  retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
202  }
203  }
204  }
205  }
206  }
207  else
208  print_usage ();
209 
210  return retval;
211 }
212 
213 /*
214 %!test
215 %! a = [.5 .5 .5 .5 .5];
216 %! x = [0 1 2 3 4];
217 %! v1 = sqrt (pi)*erf (x)./gamma (a);
218 %! v3 = gammainc (x.*x, a);
219 %! assert (v1, v3, sqrt (eps));
220 
221 %!assert (gammainc (0:4,0.5, "upper"), 1-gammainc (0:4,0.5), 1e-10)
222 
223 %!test
224 %! a = single ([.5 .5 .5 .5 .5]);
225 %! x = single ([0 1 2 3 4]);
226 %! v1 = sqrt (pi ("single"))*erf (x)./gamma (a);
227 %! v3 = gammainc (x.*x, a);
228 %! assert (v1, v3, sqrt (eps ("single")));
229 
230 %!assert (gammainc (single (0:4), single (0.5), "upper"),
231 %! single (1)-gammainc (single (0:4), single (0.5)),
232 %! single (1e-7))
233 */
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
bool is_scalar_type(void) const
Definition: ov.h:657
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
float float_value(bool frc_str_conv=false) const
Definition: ov.h:762
double gammainc(double x, double a, bool &err)
Definition: lo-specfun.cc:2575
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:782
int error_state
Definition: error.cc:101
octave_idx_type length(void) const
Definition: ov.cc:1525
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:779
ColumnVector transform(const Matrix &m, double x, double y, double z)
Definition: graphics.cc:5259
bool is_single_type(void) const
Definition: ov.h:611
double double_value(bool frc_str_conv=false) const
Definition: ov.h:759
F77_RET_T const double * x