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