* Bug 16617 fixed: gamma() extended to incomplete cases
[scilab.git] / 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