Mercurial > hg > octave-jordi
changeset 10698:fa00ccf7b1f9
Implementary complementary incomplete gamma function (bug #30088)
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Sat, 12 Jun 2010 01:01:39 +0200 |
| parents | 1215ab6f3491 |
| children | da51afafca80 |
| files | src/ChangeLog src/DLD-FUNCTIONS/gammainc.cc |
| diffstat | 2 files changed, 49 insertions(+), 9 deletions(-) [+] |
line wrap: on
line diff
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2010-06-11 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/gammainc.cc: Implement complementary incomplete + gamma function. + 2010-06-10 Ben Abbott <bpabbott@mac.com> * data.cc: Fix test for concatentating empty nd-arrays.
--- a/src/DLD-FUNCTIONS/gammainc.cc +++ b/src/DLD-FUNCTIONS/gammainc.cc @@ -36,6 +36,8 @@ DEFUN_DLD (gammainc, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} gammainc (@var{x}, @var{a})\n\ +@deftypefnx {Mapping Function} {} gammainc (@var{x}, @var{a}, \"lower\")\n\ +@deftypefnx {Mapping Function} {} gammainc (@var{x}, @var{a}, \"upper\")\n\ Compute the normalized incomplete gamma function,\n\ @tex\n\ $$\n\ @@ -61,14 +63,39 @@ \n\ If neither @var{x} nor @var{a} is scalar, the sizes of @var{x} and\n\ @var{a} must agree, and @var{gammainc} is applied element-by-element.\n\ +\n\ +By default the incomplete gamma function integrated from 0 to @var{x} is\n\ +computed. If \"upper\" is given then the complementary function integrated\n\ +for @var{x} to infinity is calculated. It should be noted that\n\ +\n\ +@example\n\ +gammainc (@var{x}, @var{a}) = 1 - gammainc (@var{x}, @vaar{a}, \"upper\")\n\ +@end example\n\ @seealso{gamma, lgamma}\n\ @end deftypefn") { octave_value retval; + bool lower = true; int nargin = args.length (); - if (nargin == 2) + if (nargin == 3) + { + if (args(2).is_string ()) + { + std::string s = args(2).string_value (); + std::transform (s.begin (), s.end (), s.begin (), tolower); + if (s == "upper") + lower = false; + else if (s != "lower") + error ("expecting third argument to be \"lower\" or \"upper\""); + } + else + error ("expecting third argument to be \"lower\" or \"upper\""); + + } + + if (!error_state && nargin >= 2 && nargin <= 3) { octave_value x_arg = args(0); octave_value a_arg = args(1); @@ -87,14 +114,16 @@ float a = a_arg.float_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) + : static_cast<float>(1) - gammainc (x, a); } else { FloatNDArray a = a_arg.float_array_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) + : static_cast<float>(1) - gammainc (x, a); } } } @@ -109,14 +138,16 @@ float a = a_arg.float_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) + : static_cast<float>(1) - gammainc (x, a); } else { FloatNDArray a = a_arg.float_array_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) + : static_cast<float>(1) - gammainc (x, a); } } } @@ -134,14 +165,14 @@ double a = a_arg.double_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) : 1. - gammainc (x, a); } else { NDArray a = a_arg.array_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) : 1. - gammainc (x, a); } } } @@ -156,14 +187,14 @@ double a = a_arg.double_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) : 1. - gammainc (x, a); } else { NDArray a = a_arg.array_value (); if (! error_state) - retval = gammainc (x, a); + retval = lower ? gammainc (x, a) : 1. - gammainc (x, a); } } } @@ -184,6 +215,8 @@ %! v3 = gammainc(x.*x,a); %! assert(v1, v3, sqrt(eps)); +%!assert (gammainc(0:4,0.5,"upper"), 1-gammainc(0:4,0.5),1e-10) + %!test %! a = single ([.5 .5 .5 .5 .5]); %! x = single([0 1 2 3 4]); @@ -191,4 +224,6 @@ %! v3 = gammainc(x.*x,a); %! assert(v1, v3, sqrt(eps('single'))); +%!assert (gammainc(single(0:4),single(0.5),"upper"), single(1)-gammainc(single(0:4),single(0.5)),single(1e-7)) + */