fix nchoosek example
[scilab.git] / scilab / modules / elementary_functions / help / en_US / discrete / nchoosek.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
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  -->
16  <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
17           xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml"
18           xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook"
19           xmlns:scilab="http://www.scilab.org" xml:id="nchoosek" xml:lang="en">
20   <refnamediv>
21     <refname>nchoosek</refname>
22     <refpurpose>
23         Computes binomial numbers (n,k) = numbers of combinations
24     </refpurpose>
25   </refnamediv>
26
27 <refsynopsisdiv>
28    <title>Syntax</title>
29    <synopsis>
30     b = nchoosek(n, k)
31     logb = nchoosek(n, k, logFormat)
32     [logb, b] = nchoosek(n, k, logFormat)
33     [logb, b] = nchoosek(n, k)
34    </synopsis>
35 </refsynopsisdiv>
36
37 <refsection>
38     <title>Arguments</title>
39     <variablelist>
40         <varlistentry>
41             <term>n, k</term>
42             <listitem>
43                 arrays of positive decimal integers. If both <varname>n</varname> and
44                 <varname>k</varname> are not scalars, they must have the same size.
45             </listitem>
46         </varlistentry>
47         <varlistentry>
48             <term>b</term>
49             <listitem>
50                 array of positive decimal integers, with the size of the biggest array
51                 <varname>n</varname> or <varname>k</varname> :
52                 C<subscript>n</subscript><superscript>k</superscript>
53             </listitem>
54         </varlistentry>
55         <varlistentry>
56             <term>logb</term>
57             <listitem>
58                 array of positive decimal numbers, with the size of <varname>b</varname>:
59                 <literal>log10(C<subscript>n</subscript><superscript>k</superscript>)</literal>
60             </listitem>
61         </varlistentry>
62         <varlistentry>
63             <term>logFormat</term>
64             <listitem>
65                 single word "log10" | "10.mant". "log10" by default.
66                 <itemizedlist>
67                     <listitem>
68                         If "log10" is used, <varname>logb</varname> returns <literal>log10(b)</literal>.
69                         <para/>
70                     </listitem>
71                     <listitem>
72                         If "10.mant" is used, then <varname>logb</varname> returns
73                         <literal>int(log10(b)) + 10^(log10(b)-int(log10(b)))/10</literal>:
74                         The 10-exponent of b is the logb integer part, while its fractional part
75                         directly shows the mantissa/10 of b, in [1.0, 10)/10.
76                     </listitem>
77                 </itemizedlist>
78             </listitem>
79         </varlistentry>
80     </variablelist>
81 </refsection>
82
83 <refsection>
84     <title>Description</title>
85     <para>
86         For every (n(i),k(i)) element-wise pair, nchoosek() computes the number of
87         k<subscript>i</subscript>-element subsets of an n<subscript>i</subscript>-element set.
88         It is mathematically defined by
89         <latex alt="C(n,k) = n! / [ k! (n-k)! ]" fontsize="18">C_n^k = {{n!} \over {k! \,(n-k)!}}
90         = {{(n-k+1)\cdot ..(n-1)\cdot n} \over {1\cdot 2\cdot ...(k-1)\cdot k}}</latex>.
91         The implementation, though, uses the βeta function.
92     </para>
93     <para>
94         If n &lt; 0, or k &lt; 0, or k &gt; n, then an error is generated.
95     </para>
96     <refsect3>
97         <title>Note about floating point accuracy</title>
98         <para>
99             Despite n! overflows for n>170, nchoosek(n,k) is computed within almost the %eps relative
100             accuracy for n much beyond 170.
101         </para>
102         <para>
103             For any given n value, we know that C<subscript>n</subscript><superscript>k</superscript>
104             is maximal for k=n/2. Scilab computes b = nchoosek(n,n/2) close to the %eps relative
105             accuracy for n up to 1029. Beyond this value, b=%inf is returned.
106         </para>
107         <para></para>
108         <para>
109             This narrow n=1029 limit can be overcome by computing
110             log10(C<subscript>n</subscript><superscript>k</superscript>) and returning
111             it through the second output argument logb. This allows to use still bigger n values,
112             and to get some valid information about the exponent and the mantissa of the true result.
113         </para>
114         <para>
115             However, we must keep in mind that beyond n = 1/%eps ~ 4.10<superscript>15</superscript>,
116             the numerical step between consecutive integers m and m+1 stored as floating point
117             numbers become > 1. Then <literal>n*(n-1)</literal> is numerically equal to
118             <literal>n*n</literal>, and getting reliable results becomes harder.
119        </para>
120         <para>
121             The integer part of logb represents the exponent (of 10) of
122             C<subscript>n</subscript><superscript>k</superscript>, whereas its fractional part
123             is the log10 of its mantissa. When the integer part of logb increases, the more digits
124             it takes among the 16 available ones in floating point encoding, the less remain
125             to encode the mantissa. Then,
126             knowing the mantissa of C<subscript>n</subscript><superscript>k</superscript> with
127             at least one significant digit requires
128             C<subscript>n</subscript><superscript>k</superscript>
129             &lt; 10<superscript>10<superscript>14</superscript></superscript>.
130        </para>
131        <para>
132             <emphasis role="bold">Relative accuracy on logb</emphasis>: As a rule of thumb,
133             except for special k cases, the relative uncertainty on logb is of the order of
134             %eps * n / k.
135        </para>
136     </refsect3>
137 </refsection>
138
139 <refsection>
140    <title>Examples</title>
141    <para>
142         Simple examples
143    </para>
144    <programlisting role="example"><![CDATA[
145 c = nchoosek(4, 1)
146 c = nchoosek(5, 0)
147 c = nchoosek([4 5], [1 0])
148    ]]></programlisting>
149    <screen><![CDATA[
150 --> c = nchoosek(4, 1)
151  c  =
152    4.
153
154 --> c = nchoosek(5, 0)
155  c  =
156    1.
157
158 --> c = nchoosek([4 5], [1 0])
159  c  =
160    4.   1.
161 ]]></screen>
162    <programlisting role="example"><![CDATA[
163 nchoosek(10, 0:10)
164 nchoosek(4:12, 4)
165    ]]></programlisting>
166    <screen><![CDATA[
167 --> nchoosek(10, 0:10)
168  ans  =
169    1.  10.  45.  120.  210.  252.  210.  120.  45.  10.  1.
170
171 --> nchoosek(4:12, 4)
172  ans  =
173    1.  5.  15.  35.  70.  126.  210.  330.  495.
174 ]]></screen>
175     <para>
176         For big values:
177     </para>
178    <programlisting role="example"><![CDATA[
179 exact = 2.050083865033972676e307;
180 c = nchoosek(10000, 134)
181 relerror = abs(c-exact)/exact
182    ]]></programlisting>
183    <screen><![CDATA[
184 --> c = nchoosek(10000, 134)
185  c  =
186    2.05D+307
187
188 --> relerror = abs(c-exact)/exact
189  relerror  =
190    7.959D-14
191 ]]></screen>
192     <para>
193         The exact result for
194         <literal>C<subscript>n</subscript><superscript>k</superscript>(n, 2)</literal> is
195         <literal>n.(n-1)/2</literal>. Now, for values <literal>n > 1/%eps = 4e15</literal>,
196         <literal>(n-1)</literal> is numerically identical to <literal>n</literal>. In no way
197         we can expect an exact result below, but rather <literal>n^2/2</literal> :
198     </para>
199    <programlisting role="example"><![CDATA[
200 n = 1e20;
201 c = nchoosek(n, 2)
202 c / (n*n/2) - 1    // Relative error wrt the expected numerical value
203    ]]></programlisting>
204    <screen><![CDATA[
205 --> c = nchoosek(n, 2)
206  c  =
207    5.000D+39
208
209 --> c /(n*n/2) - 1
210  ans  =
211   -6.661D-15
212 ]]></screen>
213     <para>
214     In logarithmic formats:
215     </para>
216    <programlisting role="example"><![CDATA[
217 [logb, b] = nchoosek(10000, 1234);  [b, logb]
218 logb = nchoosek(10000, 1234, "10.mant")
219
220 [logb, b] = nchoosek(1000, 500);
221 logb2 = nchoosek(1000, 500, "10.mant");
222 [b, logb, logb2]
223 logb = nchoosek(1000, 500, "log10")
224    ]]></programlisting>
225    <screen><![CDATA[
226 --> [logb, b] = nchoosek(10000, 1234);  [b, logb]
227  ans  =
228    Inf   1620.803261
229
230 --> logb = nchoosek(10000, 1234, "10.mant")
231  logb  =
232    1620.635713      // 6.35713D+1620
233
234 --> [logb, b] = nchoosek(1000, 500);
235 --> logb2 = nchoosek(1000, 500, "10.mant");
236 --> [b, logb, logb2]
237  ans  =
238    2.7029D+299   299.4318272   299.2702882
239
240 --> logb = nchoosek(1000, 500, "log10")
241  logb  =
242    299.4318272
243 ]]></screen>
244    <para>
245         Drawing <literal>nchoosek(n,k)</literal> on the main floating point domain, up to
246         10<superscript>300</superscript> :
247    </para>
248    <programlisting role="example"><![CDATA[
249 function ax= drawCnk(n)
250     // Preparing data
251     [N,K] = meshgrid(n);
252     noOp = K>N;
253     K(noOp) = 0;
254     [logC, C] = nchoosek(N, K);
255     C(noOp) = %nan;
256     if max(n)<2000, logC = log10(C), end
257
258     // Surface of Cnk data
259     surf(N, K, logC);
260     gce().color_mode = -2; // hidding the mesh
261     plot2d([0 n],[0 n/2]);
262     gce().children.data(:,3) = max(logC(logC<>%inf))
263
264     // Axes settings
265     xgrid(25,1,7)
266     ax = gca();
267     set(ax, "view","2d", "tight_limits",["on" "on" "off"], "grid_position","foreground");
268     xtitle("","","","");
269     ax.margins(2) = 0.05;
270
271     // Color bar
272     colorbar();
273     cb = gce();
274     cb.y_ticks.labels = msprintf("$\\Large 10^{%s}$\n", cb.y_ticks.labels);
275     title(cb,"$C_n^k$", "fontsize",3)
276 endfunction
277
278
279 clf
280 drawlater
281
282 // Figure settings
283 f = gcf();
284 f.color_map = jetcolormap(100);
285 //f.axes_size = [840 570];
286
287 // Main plot
288 ax = drawCnk(0:10:1500); sca(ax);
289 xtitle("nchoosek(n, k)","n","k","$C_n^k$");
290 set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
291 xstring(1250, 450, "Overflow")
292 gce().font_size = 4;
293
294 // Inset
295 xsetech([0.11 0.11 0.37 0.42]);
296 ax = drawCnk(0:106);
297 ax.sub_ticks = [3 3];
298 gce().sub_ticks(2) = 4;
299
300 drawnow
301    ]]></programlisting>
302    <scilab:image><![CDATA[
303 function ax= drawCnk(n)
304     // Preparing data
305     [N,K] = meshgrid(n);
306     noOp = K>N;
307     K(noOp) = 0;
308     [logC, C] = nchoosek(N, K);
309     C(noOp) = %nan;
310     if max(n)<2000, logC = log10(C), end
311
312     // Surface of Cnk data
313     surf(N, K, logC);
314     gce().color_mode = -2; // hidding the mesh
315     plot2d([0 n],[0 n/2]);
316     gce().children.data(:,3) = max(logC(logC<>%inf))
317
318     // Axes settings
319     xgrid(25,1,7)
320     ax = gca();
321     set(ax, "view","2d", "tight_limits",["on" "on" "off"], "grid_position","foreground");
322     xtitle("","","","");
323     ax.margins(2) = 0.05;
324
325     // Color bar
326     colorbar();
327     cb = gce();
328     cb.y_ticks.labels = msprintf("$\\Large 10^{%s}$\n", cb.y_ticks.labels);
329     title(cb,"$C_n^k$", "fontsize",3)
330 endfunction
331
332 clf
333 drawlater
334 // Figure settings
335 f = gcf();
336 f.color_map = jetcolormap(100);
337 f.axes_size = [840 570];
338
339 // Main plot
340 ax = drawCnk(0:10:1500); sca(ax);
341 xtitle("nchoosek(n, k)","n","k","$C_n^k$");
342 set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
343 xstring(1250, 450, "Overflow")
344 gce().font_size = 4;
345
346 // Inset
347 xsetech([0.11 0.11 0.37 0.42]);
348 ax = drawCnk(0:106);
349 ax.sub_ticks = [3 3];
350 gce().sub_ticks(2) = 4;
351 drawnow
352    ]]></scilab:image>
353    <para>
354         Going beyond, in logarithmic mode :
355    </para>
356    <programlisting role="example"><![CDATA[
357 // /!\ : The drawCnk() function used here is defined in the previous example.
358
359 clf
360 drawlater
361
362 // Figure settings
363 f = gcf();
364 f.color_map = jetcolormap(100);
365 //f.axes_size = [840 570];
366
367 // Main plot
368 ax = drawCnk(0:10000:1e6); sca(ax);
369 xtitle("nchoosek(n, k)","n","k","$C_n^k$");
370 set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
371
372 // Inset
373 xsetech([0.12 0.11 0.37 0.42]);
374 ax = drawCnk(0:1000:1e5);
375 ax.sub_ticks = [3 3];
376 gce().sub_ticks(2) = 4;
377
378 drawnow
379    ]]></programlisting>
380    <scilab:image><![CDATA[
381 function ax= drawCnk(n)
382     // Preparing data
383     [N,K] = meshgrid(n);
384     noOp = K>N;
385     K(noOp) = 0;
386     [logC, C] = nchoosek(N, K);
387     C(noOp) = %nan;
388     if max(n)<2000, logC = log10(C), end
389
390     // Surface of Cnk data
391     surf(N, K, logC);
392     gce().color_mode = -2; // hidding the mesh
393     plot2d([0 n],[0 n/2]);
394     gce().children.data(:,3) = max(logC(logC<>%inf))
395
396     // Axes settings
397     xgrid(25,1,7)
398     ax = gca();
399     set(ax, "view","2d", "tight_limits",["on" "on" "off"], "grid_position","foreground");
400     xtitle("","","","");
401     ax.margins(2) = 0.05;
402
403     // Color bar
404     colorbar();
405     cb = gce();
406     cb.y_ticks.labels = msprintf("$\\Large 10^{%s}$\n", cb.y_ticks.labels);
407     title(cb,"$C_n^k$", "fontsize",3)
408 endfunction
409
410 clf
411 drawlater
412
413 // Figure settings
414 f = gcf();
415 f.color_map = jetcolormap(100);
416 f.axes_size = [840 570];
417
418 // Main plot
419 ax = drawCnk(0:10000:1e6+1e4); sca(ax);
420 xtitle("nchoosek(n, k)","n","k","$C_n^k$");
421 set([ax.title ax.x_label ax.y_label ax.z_label], "font_size",4);
422
423 // Inset
424 xsetech([0.12 0.11 0.37 0.42]);
425 ax = drawCnk(0:1000:1e5);
426 ax.sub_ticks = [3 3];
427 gce().sub_ticks(2) = 4;
428
429 drawnow
430    ]]></scilab:image>
431    <para>
432         Top line <literal>C<subscript>n</subscript><superscript>n/2</superscript></literal> :
433         <literal>nchoosek(n, n/2)</literal> and its known close approximation
434         <latex alt="2^n /sqrt(pi.n/2)">2^n / \sqrt{\pi n/2}</latex>
435    </para>
436    <programlisting role="example"><![CDATA[
437 n = round(10^(1:0.1:8))*2;
438 logb = nchoosek(n, n/2, "log10");
439 clf
440 plot2d("ll", n, logb)
441 ax = gca();
442 xgrid(color("grey70"), 1, 7);
443 //ax.sub_ticks = [ 8 0];
444 ax.tight_limits = "on";
445 ax.y_ticks.labels = msprintf("$\\Large 10^{%d}$\n", ax.y_ticks.locations);
446 xlabel("n", "fontsize",4);
447 xstring(60, 1000, "nchoosek(n, n/2)", 270)
448 set(gce(), "clip_state", "off", "font_size",3);
449
450 // Relative difference with the 2^n / sqrt(pi.n/2) approximation:
451 e = abs(10.^(n*log10(2) - (log10(%pi)+log10(n/2))/2 - logb) - 1);
452
453 ax = newaxes();
454 ax.filled = "off";
455 ax.y_location = "right";
456 ax.tight_limits = ["on" "off"];
457 ax.font_color = color("magenta");
458 plot2d("ll", n, e, style=color("magenta"))
459 ax.axes_visible(1) = "off";
460
461 legend("$\large \left|{\frac{2^n}{\sqrt{\pi n/2}} / nchoosek(n, n/2)}-1\right|$", ..
462         "in_upper_left", %f);
463    ]]></programlisting>
464    <scilab:image><![CDATA[
465         n = round(10^(1:0.1:8))*2;
466         logb = nchoosek(n, n/2, "log10");
467         clf
468         plot2d("ll", n, logb)
469         ax = gca();
470         xgrid(color("grey70"), 1, 7);
471         //ax.sub_ticks = [ 8 0];
472         ax.tight_limits = "on";
473         ax.y_ticks.labels = msprintf("$\\Large 10^{%d}$\n", ax.y_ticks.locations);
474         xlabel("n", "fontsize",4);
475         xstring(60, 1000, "nchoosek(n, n/2)", 270)
476         set(gce(), "clip_state", "off", "font_size",3);
477
478         // Relative difference with the 2^n / sqrt(pi.n/2) approximation:
479         e = abs(10.^(n*log10(2) - (log10(%pi)+log10(n/2))/2 - logb)-1);
480
481         ax = newaxes();
482         ax.filled = "off";
483         ax.y_location = "right";
484         ax.tight_limits = ["on" "off"];
485         ax.font_color = color("magenta");
486         plot2d("ll", n, e, style=color("magenta"))
487         ax.axes_visible(1) = "off";
488
489         legend("$\large \left|{\frac{2^n}{\sqrt{\pi n/2}} / nchoosek(n, n/2)}-1\right|$", ..
490                 "in_upper_left", %f);
491    ]]></scilab:image>
492 </refsection>
493     <refsection>
494         <title>Bibliography</title>
495             <itemizedlist>
496                 <listitem>
497                     <ulink url="http://en.wikipedia.org/wiki/Binomial_coefficients">Binomial coefficients (Wikipedia)</ulink>
498                 </listitem>
499                 <listitem>
500                     Boost C++ librairies, Binomial Coefficients, 2006 , 2007, 2008, 2009, 2010 :
501                     John Maddock, Paul A. Bristow, Hubert Holin, Xiaogang Zhang, Bruno Lalande,
502                     Johan Råde, Gautam Sewani and Thijs van den Berg.
503                 </listitem>
504                 <listitem>
505                     "Introduction to discrete probabilities with Scilab", Michael Baudin, 2010
506                 </listitem>
507             </itemizedlist>
508     </refsection>
509     <refsection role="see also">
510         <title>See also</title>
511         <simplelist type="inline">
512             <member>
513                 <link linkend="factorial">factorial</link>
514             </member>
515             <member>
516                 <link linkend="gamma">gamma</link>
517             </member>
518             <member>
519                 <link linkend="perms">perms</link>
520             </member>
521         </simplelist>
522     </refsection>
523     <refsection>
524         <title>History</title>
525         <revhistory>
526             <revision>
527                 <revnumber>6.1.0</revnumber>
528                 <revdescription>
529                     <literal>nchoosek()</literal> introduced.
530                 </revdescription>
531             </revision>
532         </revhistory>
533     </refsection>
534 </refentry>