* Bug 7589 fixed: nchoosek() introduced 28/21128/10
Samuel Gougeon [Wed, 13 Nov 2019 23:24:25 +0000 (00:24 +0100)]
  http://bugzilla.scilab.org/7589

  help page (PDF): http://bugzilla.scilab.org/attachment.cgi?id=5025

Change-Id: If28b816f3481c6fd48cd8c7a0c2d499bb871b322

scilab/CHANGES.md
scilab/modules/elementary_functions/help/en_US/discrete/nchoosek.xml [new file with mode: 0644]
scilab/modules/elementary_functions/macros/nchoosek.sci [new file with mode: 0644]
scilab/modules/elementary_functions/tests/unit_tests/nchoosek.tst [new file with mode: 0644]
scilab/modules/helptools/data/pages/homepage-en_US.html
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/_LaTeX_nchoosek.xml_1.png [new file with mode: 0644]
scilab/modules/helptools/images/_LaTeX_nchoosek.xml_2.png [new file with mode: 0644]
scilab/modules/helptools/images/nchoosek_1.png [new file with mode: 0644]
scilab/modules/helptools/images/nchoosek_2.png [new file with mode: 0644]
scilab/modules/helptools/images/nchoosek_3.png [new file with mode: 0644]

index 795ec3c..460302c 100644 (file)
@@ -172,12 +172,12 @@ Feature changes and additions
   - Returned indices can now be formatted with the new option `indType`.
   - There were no unit tests. More than 100 tests are added.
 * `datafit` is now able to fit weighted data. It now supports any gap function vectorized for Data points, and so is much faster. It now accepts any matrix of parameters (not necessarily a colum). It now returns the average Mode-to-data distance, and the termination status for the quasi-Newton algo.
-* `tree_show()` is upgraded:
+* `tree_show` is upgraded:
   - New `rootTitle` and `styles` optional inputs.
   - New `arrayByFields` option, to allow displaying each array component as an object in its whole.
   - Improved layout: detailled indices for 2D arrays, simplified symbols, etc.
   - The content of implicitlist objects, and information for Scilab functions and libraries of functions are now displayed.
-
+* `nchoosek` is introduced, to compute the binomial coefficients.
 
 Help pages:
 -----------
@@ -254,6 +254,7 @@ Bug Fixes
 * [#5824](https://bugzilla.scilab.org/5824): The `datafit` algorithm was not documented.
 * [#6070](https://bugzilla.scilab.org/6070): How to make multiscaled plots was not documented.
 * [#7562](https://bugzilla.scilab.org/7562): `factorial` could use a huge memory amount even for a scalar argument.
+* [#7589](https://bugzilla.scilab.org/7589): There was no function computing the binomial coefficients.
 * [#7657](https://bugzilla.scilab.org/7657): `lstsize` was a duplicate of `size` and should be removed.
 * [#7724](https://bugzilla.scilab.org/7724): When a figure is created in .auto_resize="on" mode, its .axes_size sets its .figure_size accordingly, not the reverse. But this was not documented.
 * [#7732](https://bugzilla.scilab.org/7732): The `datafit` help page needed to be fixed and overhauled.
diff --git a/scilab/modules/elementary_functions/help/en_US/discrete/nchoosek.xml b/scilab/modules/elementary_functions/help/en_US/discrete/nchoosek.xml
new file mode 100644 (file)
index 0000000..d7fdb73
--- /dev/null
@@ -0,0 +1,534 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2009 - 2010 - DIGITEO - Michael Baudin
+ * Copyright (C) 2012 - Michael Baudin
+ * Copyright (C) 2019 - Samuel GOUGEON - Le Mans Université
+ *
+ * This file is hereby licensed under the terms of the GNU GPL v2.0,
+ * pursuant to article 5.3.4 of the CeCILL v.2.1.
+ * This file was originally licensed under the terms of the CeCILL v2.1,
+ * and continues to be available under such terms.
+ * For more information, see the COPYING file which you should have received
+ * along with this program.
+ *
+ -->
+ <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml"
+          xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
+          xmlns:scilab="http://www.scilab.org" xml:id="nchoosek" xml:lang="en">
+  <refnamediv>
+    <refname>nchoosek</refname>
+    <refpurpose>
+        Computes binomial numbers (n,k) = numbers of combinations
+    </refpurpose>
+  </refnamediv>
+
+<refsynopsisdiv>
+   <title>Syntax</title>
+   <synopsis>
+    b = nchoosek(n, k)
+    logb = nchoosek(n, k, logFormat)
+    [logb, b] = nchoosek(n, k, logFormat)
+    [logb, b] = nchoosek(n, k)
+   </synopsis>
+</refsynopsisdiv>
+
+<refsection>
+    <title>Arguments</title>
+    <variablelist>
+        <varlistentry>
+            <term>n, k</term>
+            <listitem>
+                arrays of positive decimal integers. If both <varname>n</varname> and
+                <varname>k</varname> are not scalars, they must have the same size.
+            </listitem>
+        </varlistentry>
+        <varlistentry>
+            <term>b</term>
+            <listitem>
+                array of positive decimal integers, with the size of the biggest array
+                <varname>n</varname> or <varname>k</varname> :
+                C<subscript>n</subscript><superscript>k</superscript>
+            </listitem>
+        </varlistentry>
+        <varlistentry>
+            <term>logb</term>
+            <listitem>
+                array of positive decimal numbers, with the size of <varname>b</varname>:
+                <literal>log10(C<subscript>n</subscript><superscript>k</superscript>)</literal>
+            </listitem>
+        </varlistentry>
+        <varlistentry>
+            <term>logFormat</term>
+            <listitem>
+                single word "log10" | "10.mant". "log10" by default.
+                <itemizedlist>
+                    <listitem>
+                        If "log10" is used, <varname>logb</varname> returns <literal>log10(b)</literal>.
+                        <para/>
+                    </listitem>
+                    <listitem>
+                        If "10.mant" is used, then <varname>logb</varname> returns
+                        <literal>int(log10(b)) + 10^(log10(b)-int(log10(b)))/10</literal>:
+                        The 10-exponent of b is the logb integer part, while its fractional part
+                        directly shows the mantissa/10 of b, in [1.0, 10)/10.
+                    </listitem>
+                </itemizedlist>
+            </listitem>
+        </varlistentry>
+    </variablelist>
+</refsection>
+
+<refsection>
+    <title>Description</title>
+    <para>
+        For every (n(i),k(i)) element-wise pair, nchoosek() computes the number of
+        k<subscript>i</subscript>-element subsets of an n<subscript>i</subscript>-element set.
+        It is mathematically defined by
+        <latex alt="C(n,k) = n! / [ k! (n-k)! ]" fontsize="18">C_n^k = {{n!} \over {k! \,(n-k)!}}
+        = {{(n-k+1)\cdot ..(n-1)\cdot n} \over {1\cdot 2\cdot ...(k-1)\cdot k}}</latex>.
+        The implementation, though, uses the βeta function.
+    </para>
+    <para>
+        If n &lt; 0, or k &lt; 0, or k &gt; n, then an error is generated.
+    </para>
+    <refsect3>
+        <title>Note about floating point accuracy</title>
+        <para>
+            Despite n! overflows for n>170, nchoosek(n,k) is computed within almost the %eps relative
+            accuracy for n much beyond 170.
+        </para>
+        <para>
+            For any given n value, we know that C<subscript>n</subscript><superscript>k</superscript>
+            is maximal for k=n/2. Scilab computes b = nchoosek(n,n/2) close to the %eps relative
+            accuracy for n up to 1029. Beyond this value, b=%inf is returned.
+        </para>
+        <para></para>
+        <para>
+            This narrow n=1029 limit can be overcome by computing
+            log10(C<subscript>n</subscript><superscript>k</superscript>) and returning
+            it through the second output argument logb. This allows to use still bigger n values,
+            and to get some valid information about the exponent and the mantissa of the true result.
+        </para>
+        <para>
+            However, we must keep in mind that beyond n = 1/%eps ~ 4.10<superscript>15</superscript>,
+            the numerical step between consecutive integers m and m+1 stored as floating point
+            numbers become > 1. Then <literal>n*(n-1)</literal> is numerically equal to
+            <literal>n*n</literal>, and getting reliable results becomes harder.
+       </para>
+        <para>
+            The integer part of logb represents the exponent (of 10) of
+            C<subscript>n</subscript><superscript>k</superscript>, whereas its fractional part
+            is the log10 of its mantissa. When the integer part of logb increases, the more digits
+            it takes among the 16 available ones in floating point encoding, the less remain
+            to encode the mantissa. Then,
+            knowing the mantissa of C<subscript>n</subscript><superscript>k</superscript> with
+            at least one significant digit requires
+            C<subscript>n</subscript><superscript>k</superscript>
+            &lt; 10<superscript>10<superscript>14</superscript></superscript>.
+       </para>
+       <para>
+            <emphasis role="bold">Relative accuracy on logb</emphasis>: As a rule of thumb,
+            except for special k cases, the relative uncertainty on logb is of the order of
+            %eps * n / k.
+       </para>
+    </refsect3>
+</refsection>
+
+<refsection>
+   <title>Examples</title>
+   <para>
+        Simple examples
+   </para>
+   <programlisting role="example"><![CDATA[
+c = nchoosek(4, 1)
+c = nchoosek(5, 0)
+c = nchoosek([4 5], [1 0])
+   ]]></programlisting>
+   <screen><![CDATA[
+--> c = nchoosek(4, 1)
+ c  =
+   4.
+
+--> c = nchoosek(5, 0)
+ c  =
+   1.
+
+--> c = nchoosek([4 5], [1 0])
+ c  =
+   4.   1.
+]]></screen>
+   <programlisting role="example"><![CDATA[
+nchoosek(10, 0:10)
+nchoosek(4:12, 4)
+   ]]></programlisting>
+   <screen><![CDATA[
+--> nchoosek(10, 0:10)
+ ans  =
+   1.  10.  45.  120.  210.  252.  210.  120.  45.  10.  1.
+
+--> nchoosek(4:12, 4)
+ ans  =
+   1.  5.  15.  35.  70.  126.  210.  330.  495.
+]]></screen>
+    <para>
+        For big values:
+    </para>
+   <programlisting role="example"><![CDATA[
+exact = 2.050083865033972676e307;
+c = nchoosek(10000, 134)
+relerror = abs(c-exact)/exact
+   ]]></programlisting>
+   <screen><![CDATA[
+--> c = nchoosek(10000, 134)
+ c  =
+   2.05D+307
+
+--> relerror = abs(c-exact)/exact
+ relerror  =
+   7.959D-14
+]]></screen>
+    <para>
+        The exact result for
+        <literal>C<subscript>n</subscript><superscript>k</superscript>(n, 2)</literal> is
+        <literal>n.(n-1)/2</literal>. Now, for values <literal>n > 1/%eps = 4e15</literal>,
+        <literal>(n-1)</literal> is numerically identical to <literal>n</literal>. In no way
+        we can expect an exact result below, but rather <literal>n^2/2</literal> :
+    </para>
+   <programlisting role="example"><![CDATA[
+n = 1e20;
+c = nchoosek(n, 2)
+c / (n*n/2) - 1    // Relative error wrt the expected numerical value
+   ]]></programlisting>
+   <screen><![CDATA[
+--> c = nchoosek(n, 2)
+ c  =
+   5.000D+39
+
+--> c /(n*n/2) - 1
+ ans  =
+  -6.661D-15
+]]></screen>
+    <para>
+    In logarithmic formats:
+    </para>
+   <programlisting role="example"><![CDATA[
+[logb, b] = nchoosek(10000, 1234);  [b, logb]
+logb = nchoosek(10000, 1234, "10.mant")
+
+[logb, b] = nchoosek(1000, 500);
+logb2 = nchoosek(1000, 500, "10.mant");
+[b, logb, logb2]
+logb = nchoosek(1000, 500, "log10")
+   ]]></programlisting>
+   <screen><![CDATA[
+--> [logb, b] = nchoosek(10000, 1234);  [b, logb]
+ ans  =
+   Inf   1620.803261
+
+--> logb = nchoosek(10000, 1234, "10.mant")
+ logb  =
+   1620.635713      // 6.35713D+1620
+
+--> [logb, b] = nchoosek(1000, 500);
+--> logb2 = nchoosek(1000, 500, "10.mant");
+--> [b, logb, logb2]
+ ans  =
+   2.7029D+299   299.4318272   299.2702882
+
+--> logb = nchoosek(1000, 500, "log10")
+ logb  =
+   299.4318272
+]]></screen>
+   <para>
+        Drawing <literal>nchoosek(n,k)</literal> on the main floating point domain, up to
+        10<superscript>300</superscript> :
+   </para>
+   <programlisting role="example"><![CDATA[
+function ax= drawCnk(n)
+    // Preparing data
+    [N,K] = meshgrid(n);
+    noOp = K>N;
+    K(noOp) = 0;
+    [logC, C] = nchoosek(N, K);
+    C(noOp) = %nan;
+    if max(n)<2000, logC = log10(C), end
+
+    // Surface of Cnk data
+    surf(N, K, logC);
+    gce().color_mode = -2; // hidding the mesh
+    plot2d([0 n],[0 n/2]);
+    gce().children.data(:,3) = max(logC(logC<>%inf))
+
+    // Axes settings
+    xgrid(25,1,7)
+    ax = gca();
+    set(ax, "view","2d", "tight_limits",["on" "on" "off"], "grid_position","foreground");
+    xtitle("","","","");
+    ax.margins(2) = 0.05;
+
+    // Color bar
+    colorbar();
+    cb = gce();
+    cb.y_ticks.labels = msprintf("$\\Large 10^{%s}$\n", cb.y_ticks.labels);
+    title(cb,"$C_n^k$", "fontsize",3)
+endfunction
+
+
+clf
+drawlater
+
+// Figure settings
+f = gcf();
+f.color_map = jetcolormap(100);
+//f.axes_size = [840 570];
+
+// Main plot
+ax = drawCnk(0:10:1500); sca(ax);
+xtitle("nchoosek(n, k)","n","k","$C_n^k$");
+set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
+xstring(1250, 450, "Overflow")
+gce().font_size = 4;
+
+// Inset
+xsetech([0.11 0.11 0.37 0.42]);
+ax = drawCnk(0:106);
+ax.sub_ticks = [3 3];
+gce().sub_ticks(2) = 4;
+
+drawnow
+   ]]></programlisting>
+   <scilab:image><![CDATA[
+function ax= drawCnk(n)
+    // Preparing data
+    [N,K] = meshgrid(n);
+    noOp = K>N;
+    K(noOp) = 0;
+    [logC, C] = nchoosek(N, K);
+    C(noOp) = %nan;
+    if max(n)<2000, logC = log10(C), end
+
+    // Surface of Cnk data
+    surf(N, K, logC);
+    gce().color_mode = -2; // hidding the mesh
+    plot2d([0 n],[0 n/2]);
+    gce().children.data(:,3) = max(logC(logC<>%inf))
+
+    // Axes settings
+    xgrid(25,1,7)
+    ax = gca();
+    set(ax, "view","2d", "tight_limits",["on" "on" "off"], "grid_position","foreground");
+    xtitle("","","","");
+    ax.margins(2) = 0.05;
+
+    // Color bar
+    colorbar();
+    cb = gce();
+    cb.y_ticks.labels = msprintf("$\\Large 10^{%s}$\n", cb.y_ticks.labels);
+    title(cb,"$C_n^k$", "fontsize",3)
+endfunction
+
+clf
+drawlater
+// Figure settings
+f = gcf();
+f.color_map = jetcolormap(100);
+f.axes_size = [840 570];
+
+// Main plot
+ax = drawCnk(0:10:1500); sca(ax);
+xtitle("nchoosek(n, k)","n","k","$C_n^k$");
+set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
+xstring(1250, 450, "Overflow")
+gce().font_size = 4;
+
+// Inset
+xsetech([0.11 0.11 0.37 0.42]);
+ax = drawCnk(0:106);
+ax.sub_ticks = [3 3];
+gce().sub_ticks(2) = 4;
+drawnow
+   ]]></scilab:image>
+   <para>
+        Going beyond, in logarithmic mode :
+   </para>
+   <programlisting role="example"><![CDATA[
+// /!\ : The drawCnk() function used here is defined in the previous example.
+
+clf
+drawlater
+
+// Figure settings
+f = gcf();
+f.color_map = jetcolormap(100);
+//f.axes_size = [840 570];
+
+// Main plot
+ax = drawCnk(0:10000:1e6); sca(ax);
+xtitle("nchoosek(n, k)","n","k","$C_n^k$");
+set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
+
+// Inset
+xsetech([0.12 0.11 0.37 0.42]);
+ax = drawCnk(0:1000:1e5);
+ax.sub_ticks = [3 3];
+gce().parent.sub_ticks(2) = 4;
+
+drawnow
+   ]]></programlisting>
+   <scilab:image><![CDATA[
+function ax= drawCnk(n)
+    // Preparing data
+    [N,K] = meshgrid(n);
+    noOp = K>N;
+    K(noOp) = 0;
+    [logC, C] = nchoosek(N, K);
+    C(noOp) = %nan;
+    if max(n)<2000, logC = log10(C), end
+
+    // Surface of Cnk data
+    surf(N, K, logC);
+    gce().color_mode = -2; // hidding the mesh
+    plot2d([0 n],[0 n/2]);
+    gce().children.data(:,3) = max(logC(logC<>%inf))
+
+    // Axes settings
+    xgrid(25,1,7)
+    ax = gca();
+    set(ax, "view","2d", "tight_limits",["on" "on" "off"], "grid_position","foreground");
+    xtitle("","","","");
+    ax.margins(2) = 0.05;
+
+    // Color bar
+    colorbar();
+    cb = gce().parent;
+    cb.y_ticks.labels = msprintf("$\\Large 10^{%s}$\n", cb.y_ticks.labels);
+    title(cb,"$C_n^k$", "fontsize",3)
+endfunction
+
+clf
+drawlater
+
+// Figure settings
+f = gcf();
+f.color_map = jetcolormap(100);
+f.axes_size = [840 570];
+
+// Main plot
+ax = drawCnk(0:10000:1e6+1e4); sca(ax);
+xtitle("nchoosek(n, k)","n","k","$C_n^k$");
+set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
+
+// Inset
+xsetech([0.12 0.11 0.37 0.42]);
+ax = drawCnk(0:1000:1e5);
+ax.sub_ticks = [3 3];
+gce().parent.sub_ticks(2) = 4;
+
+drawnow
+   ]]></scilab:image>
+   <para>
+        Top line <literal>C<subscript>n</subscript><superscript>n/2</superscript></literal> :
+        <literal>nchoosek(n, n/2)</literal> and its known close approximation
+        <latex alt="2^n /sqrt(pi.n/2)">2^n / \sqrt{\pi n/2}</latex>
+   </para>
+   <programlisting role="example"><![CDATA[
+n = round(10^(1:0.1:8))*2;
+logb = nchoosek(n, n/2, "log10");
+clf
+plot2d("ll", n, logb)
+ax = gca();
+xgrid(color("grey70"), 1, 7);
+//ax.sub_ticks = [ 8 0];
+ax.tight_limits = "on";
+ax.y_ticks.labels = msprintf("$\\Large 10^{%d}$\n", ax.y_ticks.locations);
+xlabel("n", "fontsize",4);
+xstring(60, 1000, "nchoosek(n, n/2)", 270)
+set(gce(), "clip_state", "off", "font_size",3);
+
+// Relative difference with the 2^n / sqrt(pi.n/2) approximation:
+e = abs(10.^(n*log10(2) - (log10(%pi)+log10(n/2))/2 - logb) - 1);
+
+ax = newaxes();
+ax.filled = "off";
+ax.y_location = "right";
+ax.tight_limits = ["on" "off"];
+ax.font_color = color("magenta");
+plot2d("ll", n, e, style=color("magenta"))
+ax.axes_visible(1) = "off";
+
+legend("$\large \left|{\frac{2^n}{\sqrt{\pi n/2}} / nchoosek(n, n/2)}-1\right|$", ..
+        "in_upper_left", %f);
+   ]]></programlisting>
+   <scilab:image><![CDATA[
+        n = round(10^(1:0.1:8))*2;
+        logb = nchoosek(n, n/2, "log10");
+        clf
+        plot2d("ll", n, logb)
+        ax = gca();
+        xgrid(color("grey70"), 1, 7);
+        //ax.sub_ticks = [ 8 0];
+        ax.tight_limits = "on";
+        ax.y_ticks.labels = msprintf("$\\Large 10^{%d}$\n", ax.y_ticks.locations);
+        xlabel("n", "fontsize",4);
+        xstring(60, 1000, "nchoosek(n, n/2)", 270)
+        set(gce(), "clip_state", "off", "font_size",3);
+
+        // Relative difference with the 2^n / sqrt(pi.n/2) approximation:
+        e = abs(10.^(n*log10(2) - (log10(%pi)+log10(n/2))/2 - logb)-1);
+
+        ax = newaxes();
+        ax.filled = "off";
+        ax.y_location = "right";
+        ax.tight_limits = ["on" "off"];
+        ax.font_color = color("magenta");
+        plot2d("ll", n, e, style=color("magenta"))
+        ax.axes_visible(1) = "off";
+
+        legend("$\large \left|{\frac{2^n}{\sqrt{\pi n/2}} / nchoosek(n, n/2)}-1\right|$", ..
+                "in_upper_left", %f);
+   ]]></scilab:image>
+</refsection>
+    <refsection>
+        <title>Bibliography</title>
+            <itemizedlist>
+                <listitem>
+                    <ulink url="http://en.wikipedia.org/wiki/Binomial_coefficients">Binomial coefficients (Wikipedia)</ulink>
+                </listitem>
+                <listitem>
+                    Boost C++ librairies, Binomial Coefficients, 2006 , 2007, 2008, 2009, 2010 :
+                    John Maddock, Paul A. Bristow, Hubert Holin, Xiaogang Zhang, Bruno Lalande,
+                    Johan Råde, Gautam Sewani and Thijs van den Berg.
+                </listitem>
+                <listitem>
+                    "Introduction to discrete probabilities with Scilab", Michael Baudin, 2010
+                </listitem>
+            </itemizedlist>
+    </refsection>
+    <refsection role="see also">
+        <title>See also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="factorial">factorial</link>
+            </member>
+            <member>
+                <link linkend="gamma">gamma</link>
+            </member>
+            <member>
+                <link linkend="perms">perms</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1.0</revnumber>
+                <revdescription>
+                    <literal>nchoosek()</literal> introduced.
+                </revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
diff --git a/scilab/modules/elementary_functions/macros/nchoosek.sci b/scilab/modules/elementary_functions/macros/nchoosek.sci
new file mode 100644 (file)
index 0000000..7cc6db0
--- /dev/null
@@ -0,0 +1,121 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+//
+// Copyright (C) 2009 - John Burkardt
+// Copyright (C) 2009 - 2010 - DIGITEO - Michael Baudin
+// Copyright (C) 2012 - Michael Baudin
+// Copyright (C) 2019 - Samuel GOUGEON - Le Mans Université
+//
+// This file is hereby licensed under the terms of the GNU GPL v2.0,
+// pursuant to article 5.3.4 of the CeCILL v.2.1.
+// This file was originally licensed under the terms of the CeCILL v2.1,
+// and continues to be available under such terms.
+// For more information, see the COPYING file which you should have received
+// along with this program.
+
+function varargout = nchoosek(n, k, logFormat)
+
+    //   Returns the binomial number (n,k) and/or its log10.
+    //
+    // Parameters
+    //   n : a matrix of floating point integers, must be positive
+    //   k : a matrix of floating point integers, must be positive
+    //
+    // Syntax
+    //   b = nchoosek(n, k)
+    //   [logb, b] = nchoosek(n, k)
+    //   [logb, b] = nchoosek(n, k, logFormat)
+    //   logb = nchoosek(n, k, logFormat)
+
+    fname = "nchoosek"
+
+    // CHECKING ARGUMENTS
+    // ==================
+    // Number of input argument
+    [lhs, rhs] = argn()
+    if rhs < 2
+        msg = _("%s: Wrong number of input arguments: %d or %d expected.\n")
+        error(msprintf(msg, fname, 2, 3))
+    end
+    // Type and value of n
+    msg1 = _("%s: Argument #%d: Decimal integer expected.\n")
+    msg2 = _("%s: Argument #%d: Non-negative integers expected.\n")
+    if type(n) <> 1 | or(round(n)<>n) then
+        error(msprintf(msg1, fname, 1))
+    end
+    if or(n<0)
+        error(msprintf(msg2, fname, 1))
+    end
+    // Type and value of k
+    if type(k) <> 1 | or(round(k)<>k) then
+        error(msprintf(msg1, fname, 2))
+    end
+    if or(k<0)
+        error(msprintf(msg2, fname, 2))
+    end
+    if length(n)>1 & length(k)>1 & or(size(n)<>size(k)) then
+        msg = _("%s: Arguments #%d and #%d: Incompatible sizes.\n")
+        error(msprintf(msg, fname, 1, 2))
+    end
+
+    i = find(n<k, 1)
+    if i <> [] then
+        msg = _("%s: n(%d) < k(%d) is forbidden.\n")
+        error(msprintf(msg, fname, i, i))
+    end
+
+    if length(k)==1 then
+        k = ones(n) * k
+    elseif length(n)==1
+        n = ones(k) * n
+    end
+    // logFormat
+    if isdef("logFormat","l") then
+        if type(logFormat) <> 10
+            msg = _("%s: Argument #%d: Text expected.\n")
+            error(msprintf(msg, fname, 3))
+        end
+        logFormat = convstr(logFormat(1))
+        if and(logFormat<>["log10" "10.mant"])
+            msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
+            error(msprintf(msg, fname, 3, "''log10'',''10.mant''"))
+        end
+    end
+
+    // PROCESSING
+    // ==========
+    b = ones(n)
+    i = find(n>k)
+    if (i<>[]) then
+        b(i) = 1 ./ ( (n(i)-k(i)) .* beta(k(i)+1, n(i)-k(i)) )
+    end
+    b(k==1) = n(k==1)  // C(n,1) = n
+    b(n==k) = 1        // C(n,n) = 1
+    b = round ( b )
+
+    if argn(1)>1 | isdef("logFormat","l")then
+        p = isinf(b)
+        logb = b;
+        logb(~p) = log10(b(~p));
+        if or(p)
+            logb(p) = (gammaln(n(p)+1) - gammaln(k(p)+1) - gammaln(n(p)-k(p))  - log(n(p)-k(p)))/log(10)
+        end
+        if isdef("logFormat","l") & logFormat=="10.mant" then
+            logb = int(logb) + 10^(logb-int(logb)-1)
+        end
+        logb(logb<0) = %nan // Too big n. Example: [c, l] = nchoosek(1e110,2)
+    end
+
+    // OUTPUT
+    // ======
+    //   b = nchoosek(n, k)
+    //   [logb, b] = nchoosek(n, k)
+    //   [logb, b] = nchoosek(n, k, logFormat)
+    //   logb = nchoosek(n, k, logFormat)
+    varargout = list()
+    if lhs==2 | isdef("logFormat","l") then
+        varargout(1) = logb
+    end
+    if lhs==1 & ~isdef("logFormat","l") | lhs==2 then
+        varargout($+1) = b
+    end
+endfunction
diff --git a/scilab/modules/elementary_functions/tests/unit_tests/nchoosek.tst b/scilab/modules/elementary_functions/tests/unit_tests/nchoosek.tst
new file mode 100644 (file)
index 0000000..39210f8
--- /dev/null
@@ -0,0 +1,113 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+//
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+// Copyright (C) 2019 - Samuel GOUGEON - Le Mans Université
+//
+// This file is hereby licensed under the terms of the GNU GPL v2.0,
+// pursuant to article 5.3.4 of the CeCILL v.2.1.
+// This file was originally licensed under the terms of the CeCILL v2.1,
+// and continues to be available under such terms.
+// For more information, see the COPYING file which you should have received
+// along with this program.
+// ===========================================================================
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+// <-- ENGLISH IMPOSED -->
+//
+// <-- Description -->
+// -------------------------
+// Unit tests for nchoosek()
+// -------------------------
+
+assert_checkequal( nchoosek( 4 , 1 ) , 4 );
+assert_checkequal( nchoosek( 5 , 0 ) , 1 );
+assert_checkequal( nchoosek( 5 , 1 ) , 5 );
+assert_checkequal( nchoosek( 5 , 2 ) , 10 );
+assert_checkequal( nchoosek( 5 , 3 ) , 10 );
+assert_checkequal( nchoosek( 5 , 4 ) , 5 );
+assert_checkequal( nchoosek( 5 , 5 ) , 1 );
+assert_checkalmostequal ( nchoosek( 10000 , 134 ) , 2.050083865033972676e307 , 1.e-13 );
+assert_checkequal( nchoosek(10,0:10) , [1,10,45,120,210,252,210,120,45,10,1] );
+assert_checkequal( nchoosek(1:6,0:5) , 1:6 );
+
+// Special cases
+// -------------
+// C(n,1) = n
+assert_checkalmostequal( nchoosek(1e20, 1), 1e20, 1e-15 );
+assert_checkalmostequal( nchoosek(1e200, 1), 1e200, 1e-13 );
+// C(n,n) = 1
+assert_checkalmostequal( nchoosek(1e20, 1e20), 1, %eps);
+assert_checkalmostequal( nchoosek(1e200, 1e200), 1, %eps);
+// C(n,n-1) = n : limit case
+assert_checkalmostequal( nchoosek(1/%eps, 1/%eps-1), 1/%eps, 20*%eps);
+// For any k < n*%eps, numerically C(n,k) = n^k/k!
+assert_checkalmostequal( nchoosek(1e20, 2), 1e40/2, 1e-13 );
+assert_checkalmostequal( nchoosek(1e100, 2),1e200/2, 1e-13 );
+assert_checkalmostequal( nchoosek(1e20, 3), 1e60/6, 1e-13 );
+assert_checkalmostequal( nchoosek(1e50, 3), 1e150/6, 1e-13 );
+
+// Log argout / Huge values
+// ------------------------
+// Refs from  W o l f r a m  alpha
+[logb, b] = nchoosek(1029, 514);
+ref = [1.42982068649890408D+308, 308.1552815761035049036];
+assert_checkalmostequal(b, ref(1), 1e-13);
+assert_checkalmostequal(logb, ref(2), 0);
+logb = nchoosek(1029, 514, "10.mant");
+assert_checkalmostequal(logb, 308.1429820686498904, 1e-15);
+[logb, b] = nchoosek(1030, 515);
+assert_checkalmostequal(b, %inf);
+assert_checkalmostequal(logb, 308.45631157176751324655, 1e-15);
+logb = nchoosek(1030, 515, "10.mant");
+assert_checkalmostequal(logb, 308.28596413729978081637, 1e-15);
+
+// 1e14, 1e4 : 33.51338111119899206051 × 10^104340
+logb = nchoosek(1e14, 1e4, "10.mant");
+assert_checkalmostequal(logb, 104340.35133826920627, 1e-6);
+// 1e13, 100 : 1.071510287595069331e1142
+logb = nchoosek(1e13, 100, "log10");
+assert_checkalmostequal(logb, 1142.0299963450692358, 1e-5);
+logb = nchoosek(1e12, 100, "log10");
+assert_checkalmostequal(logb, 1042.0299963431344539, 1e-6);
+// 1e11, 100  : 1.071510235085708956115 × 10^942
+logb = nchoosek(1e11, 100, "10.mant");
+assert_checkalmostequal(logb, 942.10715102350857089, 1e-8);
+// 1e11, 1000 : 2.485155729882870340071 × 10^8432
+logb = nchoosek(1e11, 1000, "10.mant");
+assert_checkalmostequal(logb, 8432.24851557298828703, 1e-8);
+// 1e11, 10000 : 3.51162679090294954661 × 10^74340
+logb = nchoosek(1e11, 10000, "10.mant");
+assert_checkalmostequal(logb, 74340.3511626790902950, 1e-8);
+// 1e11, 1e5 : 3.3681041688073051733304 × 10^643426
+logb = nchoosek(1e11, 1e5, "10.mant");
+assert_checkalmostequal(logb, 643426.336810416880730, 1e-9);
+// 1e11, 1e6 : 8.1533447321246152778592 × 10^5434288
+logb = nchoosek(1e11, 1e6, "10.mant");
+assert_checkalmostequal(logb, 5434288.8153344732124, 1e-10);
+
+// Check error cases
+// -----------------
+msg = "nchoosek: Wrong number of input arguments: 2 or 3 expected."
+assert_checkerror("nchoosek()", msg);
+assert_checkerror("nchoosek([])", msg );
+assert_checkerror("nchoosek(4)", msg );
+assert_checkerror("[logb, b]=nchoosek()", msg);
+assert_checkerror("[logb, b]=nchoosek(4)", msg );
+
+msg = "nchoosek: Arguments #1 and #2: Incompatible sizes.";
+assert_checkerror("nchoosek(10:12, 2:6)", msg );
+
+msg = "nchoosek: Argument #1: Decimal integer expected.";
+assert_checkerror("nchoosek( [4.5 1.5], [2 1])", msg);
+assert_checkerror("nchoosek( [4.5 1.5], [2 -1])", msg);
+msg = "nchoosek: Argument #2: Decimal integer expected.";
+assert_checkerror("nchoosek( [4  3], [2.5 1])", msg );
+assert_checkerror("nchoosek( [4  3], ""log10"")", msg );
+
+msg = "nchoosek: Argument #1: Non-negative integers expected.";
+assert_checkerror("nchoosek( [-4 3], [2  1])", msg);
+assert_checkerror("nchoosek( [-4 3], [2 -1])", msg);
+msg = "nchoosek: Argument #2: Non-negative integers expected.";
+assert_checkerror("nchoosek( [4 3], [2 -1])", msg);
+
+assert_checkerror("nchoosek( 17 , 18 )", "nchoosek: n(1) < k(1) is forbidden." );
index 271a26a..4ae4a58 100644 (file)
                 </ul>
             </li>
             <li><code>datafit</code> is now able to fit weighted data. It now supports any gap function vectorized for Data points, and so is much faster. It now accepts any matrix of parameters (not necessarily a colum). It now returns the average Mode-to-data distance, and the termination status for the quasi-Newton algo.</li>
+            <li><code>nchoosek</code> is introduced, to compute the binomial coefficients.</li>
         </ul>
         <h2 class="title">Help pages</h2>
         <ul>
-            <li>overhauled/rewritten: <code>bitget, edit, factorial</code></li>
+            <li>overhauled/rewritten: <code>bitget, edit, factorial, vectorfind, datafit</code></li>
             <li>fixed / improved:  <code>bench_run, M_SWITCH, comet, comet3d</code></li>
             <li>Rewritten: <code>weekday</code></li>
             <li>Translations added:
index 9bc17db..f0c2a89 100644 (file)
@@ -383,6 +383,8 @@ _LaTeX_lqr.xml_6.png=4f6dd247e1f511ea66298c4c2a76bb94
 _LaTeX_lqr.xml_7.png=c5cd7e3f5ea453e798ef0d5e0cbe734c
 _LaTeX_lqr.xml_8.png=1645328711c5c5406c0bff794c3afd40
 _LaTeX_lqr.xml_9.png=c76c57d5a1f43cdcb78069c14fac2e48
+_LaTeX_nchoosek.xml_1.png=135b3c20a094e91dbe86ffaf4db0ffcc
+_LaTeX_nchoosek.xml_2.png=5070bc471485a78578069bb4405225b4
 _LaTeX_number_properties.xml_1.png=e44429416209c43ee6736fc6e4fb475e
 _LaTeX_number_properties.xml_2.png=672a7aeac9254219b9609330a12e55e5
 _LaTeX_number_properties.xml_3.png=e0fec874e861bbf6bd95ae5bdd8df9e3
@@ -853,6 +855,9 @@ move_1.png=c4599a8693c9c099ac9c7f1a5dee3e91
 mrfit_1.png=d8d3736f01fc8fc707c53508cddd575f
 nanreglin_1.png=a2189dd0a4084c54b9c1e05b45db555b
 nanreglin_2.png=5a3c75b736030778ba4a1717051810cb
+nchoosek_1.png=f8a2e760e24aee2eef300957c9b5a0ab
+nchoosek_2.png=d59cb2852e6e414d9d14967de0ac3f2b
+nchoosek_3.png=f05fd5b52c25f76208d91dc4ca10af2f
 ndgrid_1.png=881b8c6e4e1ccf21fce8004766f2df8d
 ndgrid_2.png=02c9c5b309ed4b25ee57c8b35e6a1fe1
 ndgrid_ru_RU_2.png=4ea2b3c6a3035633dbbdf9ffc2c7a676
diff --git a/scilab/modules/helptools/images/_LaTeX_nchoosek.xml_1.png b/scilab/modules/helptools/images/_LaTeX_nchoosek.xml_1.png
new file mode 100644 (file)
index 0000000..de2d00f
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_nchoosek.xml_1.png differ
diff --git a/scilab/modules/helptools/images/_LaTeX_nchoosek.xml_2.png b/scilab/modules/helptools/images/_LaTeX_nchoosek.xml_2.png
new file mode 100644 (file)
index 0000000..c2dab2e
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_nchoosek.xml_2.png differ
diff --git a/scilab/modules/helptools/images/nchoosek_1.png b/scilab/modules/helptools/images/nchoosek_1.png
new file mode 100644 (file)
index 0000000..eea7ddf
Binary files /dev/null and b/scilab/modules/helptools/images/nchoosek_1.png differ
diff --git a/scilab/modules/helptools/images/nchoosek_2.png b/scilab/modules/helptools/images/nchoosek_2.png
new file mode 100644 (file)
index 0000000..14a2794
Binary files /dev/null and b/scilab/modules/helptools/images/nchoosek_2.png differ
diff --git a/scilab/modules/helptools/images/nchoosek_3.png b/scilab/modules/helptools/images/nchoosek_3.png
new file mode 100644 (file)
index 0000000..cf4af1d
Binary files /dev/null and b/scilab/modules/helptools/images/nchoosek_3.png differ