gammainc.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1997-2012 John W. Eaton
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026 
00027 #include "lo-specfun.h"
00028 
00029 #include "defun-dld.h"
00030 #include "error.h"
00031 #include "gripes.h"
00032 #include "oct-obj.h"
00033 #include "utils.h"
00034 
00035 DEFUN_DLD (gammainc, args, ,
00036   "-*- texinfo -*-\n\
00037 @deftypefn  {Mapping Function} {} gammainc (@var{x}, @var{a})\n\
00038 @deftypefnx {Mapping Function} {} gammainc (@var{x}, @var{a}, \"lower\")\n\
00039 @deftypefnx {Mapping Function} {} gammainc (@var{x}, @var{a}, \"upper\")\n\
00040 Compute the normalized incomplete gamma function,\n\
00041 @tex\n\
00042 $$\n\
00043  \\gamma (x, a) = {1 \\over {\\Gamma (a)}}\\displaystyle{\\int_0^x t^{a-1} e^{-t} dt}\n\
00044 $$\n\
00045 @end tex\n\
00046 @ifnottex\n\
00047 \n\
00048 @example\n\
00049 @group\n\
00050                                  x\n\
00051                        1        /\n\
00052 gammainc (x, a) = ---------    | exp (-t) t^(a-1) dt\n\
00053                    gamma (a)    /\n\
00054                              t=0\n\
00055 @end group\n\
00056 @end example\n\
00057 \n\
00058 @end ifnottex\n\
00059 with the limiting value of 1 as @var{x} approaches infinity.\n\
00060 The standard notation is @math{P(a,x)}, e.g., Abramowitz and Stegun (6.5.1).\n\
00061 \n\
00062 If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned\n\
00063 for each element of @var{x} and vice versa.\n\
00064 \n\
00065 If neither @var{x} nor @var{a} is scalar, the sizes of @var{x} and\n\
00066 @var{a} must agree, and @code{gammainc} is applied element-by-element.\n\
00067 \n\
00068 By default the incomplete gamma function integrated from 0 to @var{x} is\n\
00069 computed.  If \"upper\" is given then the complementary function integrated\n\
00070 from @var{x} to infinity is calculated.  It should be noted that\n\
00071 \n\
00072 @example\n\
00073 gammainc (@var{x}, @var{a}) @equiv{} 1 - gammainc (@var{x}, @var{a}, \"upper\")\n\
00074 @end example\n\
00075 @seealso{gamma, lgamma}\n\
00076 @end deftypefn")
00077 {
00078   octave_value retval;
00079   bool lower = true;
00080 
00081   int nargin = args.length ();
00082 
00083   if (nargin == 3)
00084     {
00085       if (args(2).is_string ())
00086         {
00087           std::string s = args(2).string_value ();
00088           std::transform (s.begin (), s.end (), s.begin (), tolower);
00089           if (s == "upper")
00090             lower = false;
00091           else if (s != "lower")
00092             error ("gammainc: third argument must be \"lower\" or \"upper\"");
00093         }
00094       else
00095         error ("gammainc: third argument must be \"lower\" or \"upper\"");
00096 
00097     }
00098 
00099   if (!error_state && nargin >= 2  && nargin <= 3)
00100     {
00101       octave_value x_arg = args(0);
00102       octave_value a_arg = args(1);
00103 
00104       // FIXME Can we make a template version of the duplicated code below
00105       if (x_arg.is_single_type () || a_arg.is_single_type ())
00106         {
00107           if (x_arg.is_scalar_type ())
00108             {
00109               float x = x_arg.float_value ();
00110 
00111               if (! error_state)
00112                 {
00113                   if (a_arg.is_scalar_type ())
00114                     {
00115                       float a = a_arg.float_value ();
00116 
00117                       if (! error_state)
00118                         retval = lower ? gammainc (x, a)
00119                           : static_cast<float>(1) - gammainc (x, a);
00120                     }
00121                   else
00122                     {
00123                       FloatNDArray a = a_arg.float_array_value ();
00124 
00125                       if (! error_state)
00126                         retval = lower ? gammainc (x, a)
00127                           : static_cast<float>(1) - gammainc (x, a);
00128                     }
00129                 }
00130             }
00131           else
00132             {
00133               FloatNDArray x = x_arg.float_array_value ();
00134 
00135               if (! error_state)
00136                 {
00137                   if (a_arg.is_scalar_type ())
00138                     {
00139                       float a = a_arg.float_value ();
00140 
00141                       if (! error_state)
00142                         retval = lower ? gammainc (x, a)
00143                           : static_cast<float>(1) - gammainc (x, a);
00144                     }
00145                   else
00146                     {
00147                       FloatNDArray a = a_arg.float_array_value ();
00148 
00149                       if (! error_state)
00150                         retval = lower ? gammainc (x, a)
00151                           : static_cast<float>(1) - gammainc (x, a);
00152                     }
00153                 }
00154             }
00155         }
00156       else
00157         {
00158           if (x_arg.is_scalar_type ())
00159             {
00160               double x = x_arg.double_value ();
00161 
00162               if (! error_state)
00163                 {
00164                   if (a_arg.is_scalar_type ())
00165                     {
00166                       double a = a_arg.double_value ();
00167 
00168                       if (! error_state)
00169                         retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
00170                     }
00171                   else
00172                     {
00173                       NDArray a = a_arg.array_value ();
00174 
00175                       if (! error_state)
00176                         retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
00177                     }
00178                 }
00179             }
00180           else
00181             {
00182               NDArray x = x_arg.array_value ();
00183 
00184               if (! error_state)
00185                 {
00186                   if (a_arg.is_scalar_type ())
00187                     {
00188                       double a = a_arg.double_value ();
00189 
00190                       if (! error_state)
00191                         retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
00192                     }
00193                   else
00194                     {
00195                       NDArray a = a_arg.array_value ();
00196 
00197                       if (! error_state)
00198                         retval = lower ? gammainc (x, a) : 1. - gammainc (x, a);
00199                     }
00200                 }
00201             }
00202         }
00203     }
00204   else
00205     print_usage ();
00206 
00207   return retval;
00208 }
00209 
00210 /*
00211 
00212 %!test
00213 %! a = [.5 .5 .5 .5 .5];
00214 %! x = [0 1 2 3 4];
00215 %! v1 = sqrt(pi)*erf(x)./gamma(a);
00216 %! v3 = gammainc(x.*x,a);
00217 %! assert(v1, v3, sqrt(eps));
00218 
00219 %!assert (gammainc(0:4,0.5,"upper"), 1-gammainc(0:4,0.5),1e-10)
00220 
00221 %!test
00222 %! a = single ([.5 .5 .5 .5 .5]);
00223 %! x = single([0 1 2 3 4]);
00224 %! v1 = sqrt(pi('single'))*erf(x)./gamma(a);
00225 %! v3 = gammainc(x.*x,a);
00226 %! assert(v1, v3, sqrt(eps('single')));
00227 
00228 %!assert (gammainc(single(0:4),single(0.5),"upper"), single(1)-gammainc(single(0:4),single(0.5)),single(1e-7))
00229 
00230 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines