scilab/modules/special_functions/macros/%s_gamma.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) 2018 - 2020 - Samuel GOUGEON
   3 //
   4 // This file is hereby licensed under the terms of the GNU GPL v2.0,
   5 // pursuant to article 5.3.4 of the CeCILL v.2.1.
   6 // This file was originally licensed under the terms of the CeCILL v2.1,
   7 // and continues to be available under such terms.
   8 // For more information, see the COPYING file which you should have received
   9 // along with this program.
  10
  11 function y = %s_gamma(varargin)
  12     // call for complex numbers or/and for hypermatrix of reals
  13     //      or/and for incomplete integral
  14     // gamma(a)          // with a complex or real hypermat
  15     // gamma(x, a)       // with real x >= 0
  16     // gamma(x, a, b)
  17     // gamma(x,.., option)  // "upper": complementary integral
  18
  19     msg = _("%s: Function not defined for the given argument type.\n  Check arguments or define function %s for overloading.\n")
  20     rhs = argn(2)
  21     a = varargin(min(2,rhs))
  22     sa = size(a)
  23     y = []
  24
  25     // gamma(..)  with complex numbers (all cases: complete, uncomplete, hypermat)
  26     // -------------------------------
  27     if or(type(a)==[1 5]) & ~isreal(a,0)
  28         if isdef("%s_gamma_user") & type(%s_gamma_user)==13
  29             y = %s_gamma_user(varargin(:));
  30             // Note: the overload must be able to process hypermatrices
  31             return
  32         else
  33             error(msprintf(msg, "%s_gamma", "%s_gamma_user()"))
  34         end
  35     end
  36     if a==[] then
  37         return
  38     end
  39     a = real(a)
  40
  41     // gamma(a) with hypermatrix of real numbers
  42     // -----------------------------------------
  43     if rhs==1 then
  44         y = matrix(gamma(a(:)), sa)
  45         return
  46     end
  47
  48     // gamma(x,a,..) : uncomplete gamma (all cases with real numbers)
  49     // -------------
  50     x = varargin(1)
  51     if type(x)<>1 | ~isreal(x,0) then
  52         msg = _("%s: Argument #%d: Decimal numbers expected.\n")
  53         error(msprintf(msg, "gamma", 1))
  54     end
  55     if x==[] then
  56         return
  57     end
  58     x = real(x)
  59     if min(x) <= 0 then
  60         msg = _("%s: Argument #%d: Must be > %d.\n")
  61         error(msprintf(msg, "gamma", 1, 0))
  62     end
  63     if ~isscalar(x) & ~isscalar(a) & or(size(x)<>size(a)) then
  64         msg = _("%s: Arguments #%d and #%d: Same sizes expected.\n")
  65         error(msprintf(msg, "gamma", 1, 2))
  66     end
  67     if isscalar(a) then
  68         a = a * ones(x)
  69     elseif isscalar(x)
  70         x = x * ones(a)
  71     end
  72     // "upper"
  73     upper = %f
  74     u = varargin($)
  75     if type(u)==10
  76         if convstr(u(1))<>"upper" then
  77             msg = _("%s: Argument #%d: ''%s'' expected .\n")
  78             error(msprintf(msg, "gamma", rhs, "upper"))
  79         end
  80         upper = %t
  81         varargin($) = null()
  82     end
  83     // b
  84     if length(varargin) > 2 then
  85         b = varargin(3)
  86         if type(b)<>1 | ~isreal(b,0) then
  87             msg = _("%s: Argument #%d: Decimal numbers expected.\n")
  88             error(msprintf(msg, "gamma", 3))
  89         end
  90         b = real(b)
  91         if min(b) <= 0 then
  92             msg = _("%s: Argument #%d: Must be > %d.\n")
  93             error(msprintf(msg, "gamma", 3, 0))
  94         end
  95         if ~isscalar(a) & ~isscalar(b) & or(size(a)<>size(b)) then
  96             msg = _("%s: Arguments #%d and #%d: Same sizes expected.\n")
  97             error(msprintf(msg, "gamma", 2, 3))
  98         end
  99     else
 100         b = 1
 101     end
 102     if b==[] then
 103         return
 104     end
 105     if isscalar(b) then
 106         b = b * ones(a)
 107     end
 108
 109     // PROCESSING
 110     // ==========
 111     n = ndims(x)
 112     sa = size(a)
 113     if n > 2 then
 114         [x, a, b] = (x(:), a(:), b(:))
 115     end
 116     z = find(a==0)
 117     a(z) = 1
 118
 119     if upper then
 120         [?,y] = cdfgam("PQ", x, a, b)
 121     else
 122         y = cdfgam("PQ", x, a, b)
 123     end
 124
 125     if z<>[] then
 126         y(z) = 1
 127     end
 128     if n > 2 then
 129         y = matrix(y, sa)
 130     end
 131 endfunction