* Bug 7589 fixed: nchoosek() introduced
[scilab.git] / scilab / modules / elementary_functions / macros / nchoosek.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 //
3 // Copyright (C) 2009 - John Burkardt
4 // Copyright (C) 2009 - 2010 - DIGITEO - Michael Baudin
5 // Copyright (C) 2012 - Michael Baudin
6 // Copyright (C) 2019 - Samuel GOUGEON - Le Mans Universit√©
7 //
8 // This file is hereby licensed under the terms of the GNU GPL v2.0,
9 // pursuant to article 5.3.4 of the CeCILL v.2.1.
10 // This file was originally licensed under the terms of the CeCILL v2.1,
11 // and continues to be available under such terms.
12 // For more information, see the COPYING file which you should have received
13 // along with this program.
14
15 function varargout = nchoosek(n, k, logFormat)
16
17     //   Returns the binomial number (n,k) and/or its log10.
18     //
19     // Parameters
20     //   n : a matrix of floating point integers, must be positive
21     //   k : a matrix of floating point integers, must be positive
22     //
23     // Syntax
24     //   b = nchoosek(n, k)
25     //   [logb, b] = nchoosek(n, k)
26     //   [logb, b] = nchoosek(n, k, logFormat)
27     //   logb = nchoosek(n, k, logFormat)
28
29     fname = "nchoosek"
30
31     // CHECKING ARGUMENTS
32     // ==================
33     // Number of input argument
34     [lhs, rhs] = argn()
35     if rhs < 2
36         msg = _("%s: Wrong number of input arguments: %d or %d expected.\n")
37         error(msprintf(msg, fname, 2, 3))
38     end
39     // Type and value of n
40     msg1 = _("%s: Argument #%d: Decimal integer expected.\n")
41     msg2 = _("%s: Argument #%d: Non-negative integers expected.\n")
42     if type(n) <> 1 | or(round(n)<>n) then
43         error(msprintf(msg1, fname, 1))
44     end
45     if or(n<0)
46         error(msprintf(msg2, fname, 1))
47     end
48     // Type and value of k
49     if type(k) <> 1 | or(round(k)<>k) then
50         error(msprintf(msg1, fname, 2))
51     end
52     if or(k<0)
53         error(msprintf(msg2, fname, 2))
54     end
55     if length(n)>1 & length(k)>1 & or(size(n)<>size(k)) then
56         msg = _("%s: Arguments #%d and #%d: Incompatible sizes.\n")
57         error(msprintf(msg, fname, 1, 2))
58     end
59
60     i = find(n<k, 1)
61     if i <> [] then
62         msg = _("%s: n(%d) < k(%d) is forbidden.\n")
63         error(msprintf(msg, fname, i, i))
64     end
65
66     if length(k)==1 then
67         k = ones(n) * k
68     elseif length(n)==1
69         n = ones(k) * n
70     end
71     // logFormat
72     if isdef("logFormat","l") then
73         if type(logFormat) <> 10
74             msg = _("%s: Argument #%d: Text expected.\n")
75             error(msprintf(msg, fname, 3))
76         end
77         logFormat = convstr(logFormat(1))
78         if and(logFormat<>["log10" "10.mant"])
79             msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
80             error(msprintf(msg, fname, 3, "''log10'',''10.mant''"))
81         end
82     end
83
84     // PROCESSING
85     // ==========
86     b = ones(n)
87     i = find(n>k)
88     if (i<>[]) then
89         b(i) = 1 ./ ( (n(i)-k(i)) .* beta(k(i)+1, n(i)-k(i)) )
90     end
91     b(k==1) = n(k==1)  // C(n,1) = n
92     b(n==k) = 1        // C(n,n) = 1
93     b = round ( b )
94
95     if argn(1)>1 | isdef("logFormat","l")then
96         p = isinf(b)
97         logb = b;
98         logb(~p) = log10(b(~p));
99         if or(p)
100             logb(p) = (gammaln(n(p)+1) - gammaln(k(p)+1) - gammaln(n(p)-k(p))  - log(n(p)-k(p)))/log(10)
101         end
102         if isdef("logFormat","l") & logFormat=="10.mant" then
103             logb = int(logb) + 10^(logb-int(logb)-1)
104         end
105         logb(logb<0) = %nan // Too big n. Example: [c, l] = nchoosek(1e110,2)
106     end
107
108     // OUTPUT
109     // ======
110     //   b = nchoosek(n, k)
111     //   [logb, b] = nchoosek(n, k)
112     //   [logb, b] = nchoosek(n, k, logFormat)
113     //   logb = nchoosek(n, k, logFormat)
114     varargout = list()
115     if lhs==2 | isdef("logFormat","l") then
116         varargout(1) = logb
117     end
118     if lhs==1 & ~isdef("logFormat","l") | lhs==2 then
119         varargout($+1) = b
120     end
121 endfunction