* Bugs 7562 + 15517 fixed: factorial() improved 40/19940/8
Samuel GOUGEON [Sat, 7 Apr 2018 19:12:04 +0000 (21:12 +0200)]
  http://bugzilla.scilab.org/7562
  http://bugzilla.scilab.org/15517
  help page (PDF): http://bugzilla.scilab.org/attachment.cgi?id=4702
  SEP & discussion: http://mailinglists.scilab.org/Scilab-users-SEP-factorial-extension-on-171-10-4-tp4037916.html

Change-Id: I55dbfa752f16d45705fb445092c966755e905d25

scilab/CHANGES.md
scilab/modules/elementary_functions/help/en_US/discrete/factorial.xml
scilab/modules/elementary_functions/help/ja_JP/discrete/factorial.xml [deleted file]
scilab/modules/elementary_functions/help/ru_RU/discrete/factorial.xml
scilab/modules/elementary_functions/macros/factorial.sci
scilab/modules/elementary_functions/tests/unit_tests/factorial.dia.ref [deleted file]
scilab/modules/elementary_functions/tests/unit_tests/factorial.tst
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/factorial_1.png [new file with mode: 0644]

index 362f585..25d3e8c 100644 (file)
@@ -141,13 +141,16 @@ Feature changes and additions
 * `sciargs()` returns a column instead of formerly a row.
 * Booleans and encoded integers can now be concatenated together, as in `[%f int8(-5)]`.
 * `gsort` can now perform multilevel sorting. This noticeably allows to sort completely complex numbers.
+* `factorial(n)` can be used now from n=171 up to n=10^14.
+
 
 Help pages:
 -----------
 
-* overhauled / rewritten: `bitget`, `edit`
+* overhauled / rewritten: `bitget`, `edit`, `factorial`
 * fixed / improved:  `bench_run` `M_SWITCH`, `comet`, `comet3d`
 
+
 User Interface improvements:
 ----------------------------
 
@@ -206,6 +209,7 @@ Bug Fixes
 ### Bugs fixed in 6.1.0:
 * [#2694](http://bugzilla.scilab.org/show_bug.cgi?id=2694): `bitget` did not accept positive integers of types int8, int16 or int32.
 * [#6070](http://bugzilla.scilab.org/show_bug.cgi?id=6070): Making multiscaled plots was not documented.
+* [#7562](http://bugzilla.scilab.org/show_bug.cgi?id=7562): `factorial` could use a huge memory amount even for a scalar argument.
 * [#7657](http://bugzilla.scilab.org/show_bug.cgi?id=7657): `lstsize` was a duplicate of `size` and should be removed.
 * [#7724](http://bugzilla.scilab.org/show_bug.cgi?id=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.
 * [#7765](http://bugzilla.scilab.org/show_bug.cgi?id=7765): `champ1()` is useless. `champ().colored` is available for a long time.
@@ -267,6 +271,7 @@ Bug Fixes
 * [#15431](http://bugzilla.scilab.org/show_bug.cgi?id=15431): The empty matrix `[]` and its non-convertibility were poorly documented.
 * [#15451](http://bugzilla.scilab.org/show_bug.cgi?id=15451): The code was not adapted to use `lucene 4.10` in Debian.
 * [#15514](http://bugzilla.scilab.org/show_bug.cgi?id=15514): The `set()` documentation page needed to be overhauled.
+* [#15517](http://bugzilla.scilab.org/show_bug.cgi?id=15517): `factorial(n)` could be actually used up to only n=170.
 * [#15522](http://bugzilla.scilab.org/show_bug.cgi?id=15522): `unique()` was not able to consider all Nan values as the same value. A `uniqueNan` option now allows it.
 * [#15523](http://bugzilla.scilab.org/show_bug.cgi?id=15523): `%ODEOPTIONS(1)=2` didn't work with solvers 'rk' and 'rkf'
 * [#15534](http://bugzilla.scilab.org/show_bug.cgi?id=15534): Booleans and encoded integers could not be concatenated together.
index 2309171..b7a3543 100644 (file)
@@ -1,10 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <!--
- *
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2011 - DIGITEO - Michael Baudin
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2018 - Samuel GOUGEON
  *
  * 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.
  * 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:ns3="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="factorial" xml:lang="en">
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+        xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns3="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="factorial" xml:lang="en">
     <refnamediv>
         <refname>factorial</refname>
-        <refpurpose>The factorial function</refpurpose>
+        <refpurpose>factorial function : product of the n first positive integers</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
         <synopsis>
-            f = factorial ( n )
-
+            f = factorial(n)
+            [f, p] = factorial(n)
+            [f, p, m] = factorial(n)
         </synopsis>
     </refsynopsisdiv>
-    <refsection>
-        <title>Parameters</title>
+    <refsection role="arguments">
+        <title>Arguments</title>
         <variablelist>
             <varlistentry>
-                <term>n :</term>
+                <term>n</term>
+                <listitem>
+                    <para>
+                        scalar, vector, matrix or hypermatrix of positive integers
+                        &lt;= 10<superscript>14</superscript>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>f</term>
                 <listitem>
-                    <para> a matrix of doubles. Must contain positive integers.</para>
+                    <para>
+                        array of doubles, of the sizes of n: <literal>f(i) = n(i)!</literal>.
+                        <varname>f</varname> values are exact up to <literal>n=21</literal> included.
+                    </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>f :</term>
+                <term>p</term>
                 <listitem>
-                    <para> a matrix of doubles.</para>
+                    <para>
+                        array of doubles, of the sizes of n: power of 10 of <varname>f</varname>:
+                        <literal>p(i) = int(log10(f(i)!))</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>m</term>
+                <listitem>
+                    <para>
+                        array of doubles in [1,10[, of the sizes of n: Mantissae of <varname>f</varname>,
+                        such that <literal>n(i)! = m(i) * 10^p(i)</literal>.
+                    </para>
                 </listitem>
             </varlistentry>
         </variablelist>
     </refsection>
-    <refsection>
+    <refsection role="description">
         <title>Description</title>
         <para>
-            Returns the factorial of n, that is, the product of all
-            integers 1 * 2 * ... * n.
-        </para>
-        <para>
-            This function overflows as soon as n&gt;170.
+            Returns the factorial of n, that is the product 1 * 2 * ... * n.
         </para>
+        <warning>
+            <varname>f</varname> overflows as soon as n&gt;170 and always returns %inf for any
+            bigger n.
+        </warning>
+        <note>
+            <itemizedlist>
+                <listitem>
+                    For <literal>n in [22, 170]</literal>, the relative accuracy of
+                    <varname>f</varname> is roughly <literal>%eps ~ 2e-16</literal>.
+                </listitem>
+                <listitem>
+                    For <literal>n in [171, 1.0x10<superscript>14</superscript>]</literal>, the
+                    power <varname>p</varname> value is exact, and the relative accuracy of the
+                    mantissa <varname>m</varname> goes roughly as
+                    <literal>n*%eps ~ n * 1e-16</literal> (see the last example).
+                </listitem>
+                <listitem>
+                    Beyond <literal>n > 10.0<superscript>14</superscript></literal>,
+                    <varname>p</varname> becomes > 1/%eps and gets truncated. It is then no longer
+                    possible to retrieve a reliable mantissa.
+                </listitem>
+            </itemizedlist>
+        </note>
+               <para/>
     </refsection>
     <refsection>
+        <title>Graph</title>
+        <scilab:image>
+            x = (10^(0:13)).*.(1:9); x(1)=[]; x($)=1e14;
+            [f, p, m] = factorial(x);
+            plot2d("ll", x, p+log10(m))
+            xlabel("n", "fontsize",3)
+            title("$\mathsf{log_{10}(n!)}$", "fontsize", 4)
+            xgrid(color("grey70"), 1, 7)
+            set(gca(), "sub_ticks",[8 2], "tight_limits","on");
+            gca().data_bounds([1 4]) = [1 2e15];
+            gcf().axes_size = [840 470];
+        </scilab:image>
+    </refsection>
+    <refsection role="examples">
         <title>Examples</title>
+        <para>Table of the first n! exact values :</para>
         <programlisting role="example"><![CDATA[
-// Make a table of factorial
-n = (0:30)';
+format(22);
+n = (0:21)';
 [n factorial(n)]
-
-// See the limits of factorial: f(171)=%inf
+format(10);
+   ]]></programlisting>
+   <screen><![CDATA[
+--> [n factorial(n)]
+ ans  =
+   0.    1.
+   1.    1.
+   2.    2.
+   3.    6.
+   4.    24.
+   5.    120.
+   6.    720.
+   7.    5040.
+   8.    40320.
+   9.    362880.
+   10.   3628800.
+   11.   39916800.
+   12.   479001600.
+   13.   6227020800.
+   14.   87178291200.
+   15.   1307674368000.
+   16.   20922789888000.
+   17.   355687428096000.
+   18.   6402373705728000.
+   19.   121645100408832000.
+   20.   2432902008176640000.
+   21.   51090942171709440000.
+]]></screen>
+        <para>Ceiling of factorial() in floating point representation:</para>
+      <programlisting role="example"><![CDATA[
 factorial(170) // 7.257415615307998967e306
 factorial(171) // %inf
-
-// Plot the function on all its range.
-scf();
-plot ( 1:170 , factorial , "b-o" )
-h = gcf();
-h.children.log_flags="nln";
-
    ]]></programlisting>
+        <para>Plot the function on its whole range:</para>
+      <programlisting role="example"><![CDATA[
+x = (10^(0:13)).*.(1:9); x(1)=[]; x($)=1e14;
+[f, p, m] = factorial(x);
+clf
+plot2d("ll", x, p+log10(m))
+xlabel("n", "fontsize",3)
+title("$\mathsf{log_{10}(n!)}$", "fontsize", 4)
+xgrid(color("grey70"), 1, 7)
+set(gca(), "sub_ticks",[8 2], "tight_limits","on");
+gca().data_bounds([1 4]) = [1 2e15];
+gcf().axes_size = [850 480];
+   ]]></programlisting>
+        <para>Relative factorial() errors:</para>
+      <programlisting role="example"><![CDATA[
+n = 10^(1:14)';
+[f, p, m] = factorial(n);
+// Exact (truncated) mantissae for n = 10^(1:14) :
+m0 = [
+   3.6288000000000000  9.3326215443944153  4.0238726007709377 ..  // n = 10     100   1000
+   2.8462596809170545  2.8242294079603479                     ..  // n = 10000  100000
+   8.2639316883312401  1.2024233400515904  1.6172037949214624 ..  // n = 1e6    1e7   1e8
+   9.9046265792229937  2.3257962056730834  3.7489285991050270 ..  // n = 1e9    1e10  1e11
+   1.4036611603737561  2.4033300843401153  1.6456020559872979     // n = 1e12   1e13  1e14
+   ]';
+r_err = m./m0 - 1;
+[n r_err]
+   ]]></programlisting>
+   <screen><![CDATA[
+--> [n r_err]
+ ans  =
+   10.          0.
+   100.        -5.551D-16
+   1000.        1.132D-13
+   10000.       1.918D-12
+   100000.      6.611D-12
+   1000000.     9.962D-11
+   10000000.    5.048D-08
+   100000000.   1.050D-08
+   1.000D+09   0.0000001
+   1.000D+10   0.0000019
+   1.000D+11   0.0000062
+   1.000D+12   0.0001327
+   1.000D+13   0.0004839
+   1.000D+14   0.0071116
+]]></screen>
     </refsection>
-    <refsection>
-        <title>Bibliography</title>
-        <para>
-            <ulink url="http://www.scilab.org/en/support/documentation/tutorials">Introduction to discrete probabilities</ulink>, Michael Baudin
-        </para>
+    <refsection role="see also">
+        <title>See Also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="cumprod">cumprod</link>
+            </member>
+            <member>
+                <link linkend="gamma">gamma</link>
+            </member>
+            <member>
+                <link linkend="gammaln">gammaln</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection role="history">
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1</revnumber>
+                <revdescription>
+                    Extension up to n = 10<superscript>14</superscript>. <varname>p</varname>
+                    10-power and <varname>m</varname> mantissa output added.
+                </revdescription>
+            </revision>
+        </revhistory>
     </refsection>
 </refentry>
diff --git a/scilab/modules/elementary_functions/help/ja_JP/discrete/factorial.xml b/scilab/modules/elementary_functions/help/ja_JP/discrete/factorial.xml
deleted file mode 100644 (file)
index ec834b1..0000000
+++ /dev/null
@@ -1,76 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- *
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2011 - DIGITEO - Michael Baudin
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * 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:ns3="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="factorial" xml:lang="ja">
-    <refnamediv>
-        <refname>factorial</refname>
-        <refpurpose>階乗関数</refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>呼び出し手順</title>
-        <synopsis>
-            f = factorial ( n )
-        </synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>パラメータ</title>
-        <variablelist>
-            <varlistentry>
-                <term>n :</term>
-                <listitem>
-                    <para>doubleの行列. 値は正の整数である必要があります.</para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>f :</term>
-                <listitem>
-                    <para> doubleの行列.</para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>説明</title>
-        <para>
-            nの階乗,つまり,整数 1 * 2 * ... * nの積を返します.
-        </para>
-        <para>
-            この関数は n &gt; 170 となるとオーバーフローします.
-        </para>
-    </refsection>
-    <refsection>
-        <title>例</title>
-        <programlisting role="example"><![CDATA[
-// Make a table of factorial
-n = (0:30)';
-[n factorial(n)]
-// See the limits of factorial: f(171)=%inf
-factorial(170) // 7.257415615307998967e306
-factorial(171) // %inf
-// Plot the function on all its range.
-scf();
-plot ( 1:170 , factorial , "b-o" )
-h = gcf();
-h.children.log_flags="nln";
-   ]]></programlisting>
-    </refsection>
-    <refsection>
-        <title>参考文献</title>
-        <para>
-            <ulink url="http://www.scilab.org/en/support/documentation/tutorials">Introduction to discrete probabilities</ulink>, Michael Baudin
-        </para>
-    </refsection>
-</refentry>
index b2d743d..f3edcd2 100644 (file)
@@ -1,10 +1,7 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <!--
- *
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2011 - DIGITEO - Michael Baudin
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2018 - Samuel GOUGEON
  *
  * 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.
  * 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:ns3="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="factorial" xml:lang="ru">
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+        xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns3="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="factorial" xml:lang="ru">
     <refnamediv>
         <refname>factorial</refname>
-        <refpurpose>Функция факториал</refpurpose>
+        <refpurpose>функция факториала : произведение первых n положительных целых чисел</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Синтаксис</title>
         <synopsis>
-            f = factorial ( n )
-
+            f = factorial(n)
+            [f, p] = factorial(n)
+            [f, p, m] = factorial(n)
         </synopsis>
     </refsynopsisdiv>
-    <refsection>
-        <title>Параметры</title>
+    <refsection role="arguments">
+        <title>Аргументы</title>
         <variablelist>
             <varlistentry>
-                <term>n :</term>
+                <term>n</term>
+                <listitem>
+                    <para>
+                        скаляр, вектор, матрица или гиперматрица положительных целых чисел
+                        &lt;= 10<superscript>14</superscript>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>f</term>
+                <listitem>
+                    <para>
+                        массив чисел типа double того же размера, что и <literal>n</literal>:
+                        <literal>f(i) = n(i)!</literal>.
+                        Значения <varname>f</varname> могут быть вплоть до <literal>n=21</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>p</term>
                 <listitem>
                     <para>
-                        Матрица чисел удвоенной точности (double). Числа должны быть положительными целыми.
+                        массив чисел типа double того же размера, что и <literal>n</literal>:
+                        10 в степени <varname>f</varname>: <literal>p(i) = int(log10(f(i)!))</literal>.
                     </para>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>f :</term>
+                <term>m</term>
                 <listitem>
-                    <para> матрица чисел удвоенной точности.</para>
+                    <para>
+                        массив чисел типа double в интервале <literal>[1,10)</literal> того же размера, что и <literal>n</literal>:
+                        мантисса у <varname>f</varname> такая, что <literal>n(i)! = m(i) * 10^p(i)</literal>.
+                    </para>
                 </listitem>
             </varlistentry>
         </variablelist>
     </refsection>
-    <refsection>
+    <refsection role="description">
         <title>Описание</title>
         <para>
-            Возвращает факториал <literal>n</literal>, то есть, произведение всех целых чисел
-            1 * 2 * ... * n.
-        </para>
-        <para>
-            Эта функция переполняет разрядную сетку при n&gt;170.
+            Возвращает факториал от <literal>n</literal>, равный произведению <literal>1 * 2 * ... * n</literal>.
         </para>
+        <warning>
+            <varname>f</varname> превышает разрядную сетку при <literal>n&gt;170</literal> и всегда возвращает <literal>%inf</literal> для любого большего <literal>n</literal>.
+        </warning>
+        <note>
+            <itemizedlist>
+                <listitem>
+                    Для <literal>n</literal> в <literal>[22, 170]</literal> относительная точность
+                    <varname>f</varname> примерно равна <literal>%eps ~ 2e-16</literal>.
+                </listitem>
+                <listitem>
+                    Для <literal>n</literal> в <literal>[171, 1.0x10<superscript>14</superscript>]</literal> значение
+                    степени <varname>p</varname> равна точно и относительная точность мантиссы <varname>m</varname> 
+                    примерно равно <literal>n*%eps ~ n * 1e-16</literal> (см. последний пример).
+                </listitem>
+                <listitem>
+                    После <literal>n > 10.0<superscript>14</superscript></literal>,
+                    <varname>p</varname> становится > 1/%eps и ограничивается. Тогда нельзя более
+                    получить надёжную мантиссу.
+                </listitem>
+            </itemizedlist>
+        </note>
+               <para/>
     </refsection>
     <refsection>
+        <title>Граф</title>
+        <scilab:image>
+            x = (10^(0:13)).*.(1:9); x(1)=[]; x($)=1e14;
+            [f, p, m] = factorial(x);
+            plot2d("ll", x, p+log10(m))
+            xlabel("n", "fontsize",3)
+            title("$\mathsf{log_{10}(n!)}$", "fontsize", 4)
+            xgrid(color("grey70"), 1, 7)
+            set(gca(), "sub_ticks",[8 2], "tight_limits","on");
+            gca().data_bounds([1 4]) = [1 2e15];
+            gcf().axes_size = [840 470];
+        </scilab:image>
+    </refsection>
+    <refsection role="examples">
         <title>Примеры</title>
+        <para>Таблица первых точных значений факториала n! :</para>
         <programlisting role="example"><![CDATA[
-// Делаем таблицу факториала
-n = (0:30)';
+format(22);
+n = (0:21)';
 [n factorial(n)]
-
-// Смотрим пределы факториала: f(171)=%inf
+format(10);
+   ]]></programlisting>
+   <screen><![CDATA[
+--> [n factorial(n)]
+ ans  =
+   0.    1.
+   1.    1.
+   2.    2.
+   3.    6.
+   4.    24.
+   5.    120.
+   6.    720.
+   7.    5040.
+   8.    40320.
+   9.    362880.
+   10.   3628800.
+   11.   39916800.
+   12.   479001600.
+   13.   6227020800.
+   14.   87178291200.
+   15.   1307674368000.
+   16.   20922789888000.
+   17.   355687428096000.
+   18.   6402373705728000.
+   19.   121645100408832000.
+   20.   2432902008176640000.
+   21.   51090942171709440000.
+]]></screen>
+        <para>Округление значения factorial() в большую сторону в представлении с плавающей запятой:</para>
+      <programlisting role="example"><![CDATA[
 factorial(170) // 7.257415615307998967e306
 factorial(171) // %inf
-
-// Построим на графике функцию во всём её диапазоне.
-scf();
-plot ( 1:170 , factorial , "b-o" );xgrid
-h = gcf();
-h.children.log_flags="nln";
-
    ]]></programlisting>
+        <para>Построение графика функции во всём её интервале:</para>
+      <programlisting role="example"><![CDATA[
+x = (10^(0:13)).*.(1:9); x(1)=[]; x($)=1e14;
+[f, p, m] = factorial(x);
+clf
+plot2d("ll", x, p+log10(m))
+xlabel("n", "fontsize",3)
+title("$\mathsf{log_{10}(n!)}$", "fontsize", 4)
+xgrid(color("grey70"), 1, 7)
+set(gca(), "sub_ticks",[8 2], "tight_limits","on");
+gca().data_bounds([1 4]) = [1 2e15];
+gcf().axes_size = [850 480];
+   ]]></programlisting>
+        <para>Относительные ошибки factorial():</para>
+      <programlisting role="example"><![CDATA[
+n = 10^(1:14)';
+[f, p, m] = factorial(n);
+// Точная (ограниченная) мантисса для n = 10^(1:14) :
+m0 = [
+   3.6288000000000000  9.3326215443944153  4.0238726007709377 ..  // n = 10     100   1000
+   2.8462596809170545  2.8242294079603479                     ..  // n = 10000  100000
+   8.2639316883312401  1.2024233400515904  1.6172037949214624 ..  // n = 1e6    1e7   1e8
+   9.9046265792229937  2.3257962056730834  3.7489285991050270 ..  // n = 1e9    1e10  1e11
+   1.4036611603737561  2.4033300843401153  1.6456020559872979     // n = 1e12   1e13  1e14
+   ]';
+r_err = m./m0 - 1;
+[n r_err]
+   ]]></programlisting>
+   <screen><![CDATA[
+--> [n r_err]
+ ans  =
+   10.          0.
+   100.        -5.551D-16
+   1000.        1.132D-13
+   10000.       1.918D-12
+   100000.      6.611D-12
+   1000000.     9.962D-11
+   10000000.    5.048D-08
+   100000000.   1.050D-08
+   1.000D+09   0.0000001
+   1.000D+10   0.0000019
+   1.000D+11   0.0000062
+   1.000D+12   0.0001327
+   1.000D+13   0.0004839
+   1.000D+14   0.0071116
+]]></screen>
     </refsection>
-    <refsection>
-        <title>Литература</title>
-        <para>"Introduction to discrete probabilities", Michael Baudin, 2011</para>
+    <refsection role="see also">
+        <title>Смотрите также</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="cumprod">cumprod</link>
+            </member>
+            <member>
+                <link linkend="gamma">gamma</link>
+            </member>
+            <member>
+                <link linkend="gammaln">gammaln</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection role="history">
+        <title>История</title>
+        <revhistory>
+            <revision>
+                <revnumber>6.1</revnumber>
+                <revdescription>
+                    Расширение вплоть до n = 10<superscript>14</superscript>.
+                    Добавлено <varname>p</varname> в 10-й степени и вывод мантиссы <varname>m</varname>.
+                </revdescription>
+            </revision>
+        </revhistory>
     </refsection>
-</refentry>
+</refentry>
\ No newline at end of file
index 59b5517..bd36f09 100644 (file)
@@ -1,5 +1,6 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) INRIA - Farid BELAHCENE
+// Copyright (C) 2018 - Samuel GOUGEON
 //
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
 //
 // For more information, see the COPYING file which you should have received
 // along with this program.
 
-function y=factorial(n)
+function [f, p, m] = factorial(n)
 
-    // This function returns the factorial n. If n is a vector, matrix or hypermatrix it returns the element wise factorial
+    // This function returns the factorial n. If n is a vector, matrix or
+    //  hypermatrix it returns the element wise factorial
     // Input : n, a scalar/vector/matrix/hypermat of positives integers.
-    // Output : y, a scalar/vector/matrix/hypermat
+    // Output : 
+    //  - f: the factorials
+    //  - p: the powers of 10
+    //  - m: the mantissae, in [1, 10[
 
-    rhs = argn(2);
+    [lhs, rhs] = argn();
 
+    // CHECKING INPUT ARGUMENTS
+    // ------------------------
     if rhs <> 1 then
-        error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"factorial",1));
+        msg = gettext("%s: Wrong number of input argument(s): %d expected.\n")
+        error(msprintf(msg, "factorial", 1));
     end
 
-    if (type(n) <> 1) | (~isempty(n) & (or((n-floor(n)<>0)) | or(n<0))) then
-        error(msprintf(gettext("%s: Wrong value for input argument #%d: Scalar/vector/matrix/hypermatrix of positive integers expected.\n"),"factorial",1));
-    elseif isempty(n)
-        y=n
+    if (type(n) <> 1) | (n~=[] & (or((n-floor(n)<>0)) | or(n<0))) then
+        msg = gettext("%s: Argument #%d: Non-negative integers expected.\n")
+        error(msprintf(msg, "factorial", 1));
+    end
+
+    // TRIVIAL CASE
+    if n==[]
+        f = []
         return
-    else
-        n(n==0)=1
-        ntemp=cumprod(1:max(n))
-        y=matrix(ntemp(n(:)),size(n))
     end
 
+    // PROCESSING
+    // ----------
+    f = n
+    f(n==0) = 1
+    f(n>170) = %inf;
+    k = find(n>0 & n<171);
+    if k~=[] then
+        ntemp = cumprod(1:max(n(k)))
+        f(k) = ntemp(n(k))
+        //f(k) = gamma(n(k)+1)  // slower
+    end
+
+    if lhs>1 then
+        p = n
+        m = n
+        p(n==0) = 0
+        m(n==0) = 1
+        if k~=[]
+            p(k) = int(log10(f(k)));
+            m(k) = f(k) ./ (10 .^(p(k)))
+        end
+
+        // General Stirling formula:
+        // n! ~ sqrt(2*%pi*n)*(n/%e)^n * A
+        // with A = sum(C ./ (n.^0:5))
+        // and  C = [1, 1/12, 1/288, -139/51840, -571/2488320, 163879/209018880]
+        // log(n!) = (log(2*%pi)+log(n))/2 + n*[log(n)-1] + log(A)
+        //         =  log(2*%pi)/2 + log(n)*(n+1/2) - n + log(A)
+        k = find(n>170)
+        if k ~= []
+            // (a straightforward processing with gammaln(n+1) is less accurate
+            // than the following, by a factor of 1 to 100 (see unit tests))
+            num = [1,  1,   1,  -139,    -571,    163879,     5246819,   -534703531]
+            den = [1, 12, 288, 51840, 2488320, 209018880, 75246796800, 902961561600]
+            // Next terms are listed herebelow. They don't increase the accuracy within %eps
+            //num = [num,    -4483131259,    432261921612371,     6232523202521089]
+            //den = [den, 86684309913600, 514904800886784000, 86504006548979712000]
+            c = num ./ den
+            q = length(c);
+            A = sum((ones(length(k),1)*c) ./ (n(k)(:)*ones(1,q)).^(ones(length(k),1)*(0:q-1)),"c");
+            // Careful processing of the  fractional part: we try to always
+            // keep decimals with an integer part as small as possible (and so
+            // the biggest number of fractional digits as possible):
+            tmp = log10(n(k)(:))
+            in = int(tmp).*n(k)(:)
+            tmp2 = (tmp - int(tmp)).*n(k)(:);
+            in = in + int(tmp2);
+            fr = tmp2 - int(tmp2);
+            
+            tmp2 = tmp/2;
+            in = in + int(tmp2);
+            fr = fr + (tmp2 - int(tmp2)); 
+            
+            tmp = n(k)(:)/log(10);
+            in = in - int(tmp);
+            fr = fr - (tmp-int(tmp)) + log10(2*%pi)/2 + log10(A);
+            in = in + floor(fr);
+            fr = fr - floor(fr);
+            p(k) = in;
+            m(k) = 10 .^ fr
+            m = matrix(m, size(n));
+        end
+    end
 endfunction
diff --git a/scilab/modules/elementary_functions/tests/unit_tests/factorial.dia.ref b/scilab/modules/elementary_functions/tests/unit_tests/factorial.dia.ref
deleted file mode 100644 (file)
index 3508cd2..0000000
+++ /dev/null
@@ -1,42 +0,0 @@
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
-//
-// Copyright (C) 2012 - 2016 - Scilab Enterprises
-//
-// 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 -->
-x = factorial ( 1 );
-assert_checkequal ( x , 1 );
-x = factorial ( 2 );
-assert_checkequal ( x , 2 );
-x = factorial ( [1 2 3 4] );
-assert_checkequal ( x , [1 2 6 24] );
-x = factorial ( 170 );
-assert_checkalmostequal ( x , 7.25741561530799896739e306 , 10 * %eps );
-// Test with a matrix
-n = [
-1 2 3
-4 5 6
-7 8 9
-];
-x = factorial ( n );
-e = [
-    1.       2.        6.       
-    24.      120.      720.     
-    5040.    40320.    362880.  
-];
-assert_checkequal ( x , e );
-// Test with an hypermatrix
-clear n;
-clear e;
-n(1,1,1,1:2)=[1 2];
-x = factorial ( n );
-e(1,1,1,1:2)=[1 2];
-assert_checkequal ( x , e );
-clear n;
-clear e;
index 215936f..4d7ab3c 100644 (file)
@@ -1,5 +1,6 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+// Copyright (C) 2018 - Samuel GOUGEON
 //
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
 //
 // 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 -->
 
 x = factorial ( 1 );
 assert_checkequal ( x , 1 );
-x = factorial ( 2 );
-assert_checkequal ( x , 2 );
 x = factorial ( [1 2 3 4] );
 assert_checkequal ( x , [1 2 6 24] );
+x = factorial ( [1 2 3 4]' );
+assert_checkequal ( x , [1 2 6 24]' );
 x = factorial ( 170 );
 assert_checkalmostequal ( x , 7.25741561530799896739e306 , 10 * %eps );
 // Test with a matrix
@@ -29,18 +30,109 @@ n = [
 ];
 x = factorial ( n );
 e = [
-    1.       2.        6.       
-    24.      120.      720.     
-    5040.    40320.    362880.  
+    1.       2.        6.
+    24.      120.      720.
+    5040.    40320.    362880.
 ];
 assert_checkequal ( x , e );
+
 // Test with an hypermatrix
-clear n;
-clear e;
-n(1,1,1,1:2)=[1 2];
-x = factorial ( n );
-e(1,1,1,1:2)=[1 2];
-assert_checkequal ( x , e );
-clear n;
-clear e;
+clear n expected
+n(2,2,2) = 1;
+n(1:8) = 1:8;
+expected(2,2,2) = 1;
+expected(1:8) = e'(1:8)
+assert_checkequal(factorial(n), expected);
+
+// -----------
+// Accurate values from Wolfram Alpha
+expected = [
+   1.
+   1.
+   2.
+   6.
+   24.
+   120.
+   720.
+   5040.
+   40320.
+   362880.
+   3628800.
+   39916800.
+   479001600.
+   6227020800.
+   87178291200.
+   1307674368000.
+   20922789888000.
+   355687428096000.
+   6402373705728000.
+   121645100408832000.
+   2432902008176640000.
+   51090942171709440000.
+   1124000727777607700000.
+   25852016738884978000000.
+ ];
+assert_checkequal(factorial(0:23)', expected);
+
+expected = [
+   2.652528598121910D+32     8.1591528324789774D+47   3.04140932017133780D+64 ..
+   1.19785716699698918D+100  9.3326215443944153D+157  5.7133839564458546D+262
+   ];
+assert_checkalmostequal(factorial([30 40 50 70 100 150]), expected, 5*%eps);
+
+[f, p, m] = factorial(171);
+assert_checkequal(f, %inf);
+assert_checkequal(p, 309);
+assert_checkalmostequal(m, 1.2410180702176678, 1.5e-13); // 2e-14
+
+[f, p, m] = factorial(300);
+assert_checkequal(f, %inf);
+assert_checkequal(p, 614);
+assert_checkalmostequal(m, 3.060575122164406360, 2e-13); // 1e-13
+
+[f, p, m] = factorial(200:100:900);
+m_exp = [7.886578673647905036  3.06057512216440636  6.40345228466238953 ..
+         1.220136825991110068  1.26557231622543074  2.42204012475027218 ..
+         7.710530113353860041  6.75268022096458416
+        ];
+p_exp = [ 374  614  868  1134  1408  1689  1976  2269  ];
+assert_checkequal(p, p_exp);
+assert_checkalmostequal(m, m_exp, 5e-13);
+
+[f, p, m] = factorial(1000:1000:1e4);
+m_exp = [
+   4.02387260077093774  ..
+   3.31627509245063324   4.149359603437854086   1.828801951514065013 ..
+   4.22857792660554352   2.683999765726739596   8.842007956963112247 ..
+   5.18418106048087694   8.099589986687190858   2.846259680917054519
+    ];
+p_exp = [ 2567  5735  9130  12673  16325  20065  23877  27752  31681  35659 ];
+assert_checkequal(p, p_exp);
+assert_checkalmostequal(m, m_exp, 6e-12);  // 1.2e-11
+
+[f, p, m] = factorial(2e4:1e4:1e5);
+m_exp = [
+   1.81920632023034513   2.75953724621938460   2.09169242212132363  ..
+   3.34732050959714483   1.56413770802606692   1.17681241537969008 ..
+   3.09772225166224929   1.58011915409779537   2.82422940796034787
+    ];
+p_exp = [ 77337  121287  166713  213236  260634  308759  357506  406798  456573 ];
+assert_checkequal(p, p_exp);
+assert_checkalmostequal(m, m_exp, 1.5e-10);  //3e-10
 
+[f, p, m] = factorial(round(10^(6:14)));
+m_exp = [
+   8.26393168833124006  1.202423340051590346  1.61720379492146239  ..
+   9.90462657922299373  2.325796205673083365  3.74892859910502696  ..
+   1.40366116037375609  2.403330084340115344  1.64560205598729788
+   ];
+p_exp = [5565708  65657059  756570556  8565705522  95657055186 ..
+         1056570551815  11565705518103  125657055180974  1356570551809682
+        ];
+//r_err = [3e-9 3e-8 3e-7 2e-6 5e-6 4e-4 5e-3 7e-2 2e-1]; // with 1st implementation
+//r_err = [3e-9 4e-8 2e-7 2e-6 3e-5 2e-4 2e-3 5e-2 8e-2]; // with gammaln(n+1)
+r_err  = [1e-10 6e-8 2e-8 2e-7 2e-6 7e-6 2e-4 5e-4 8e-3];
+for i = 1:length(m_exp)
+    assert_checkequal(p(i), p_exp(i));
+    assert_checkalmostequal(m(i), m_exp(i), r_err(i));
+end
index dc18f73..1fbd61c 100644 (file)
@@ -302,7 +302,7 @@ _LaTeX_floor.xml_ru_RU_1.png=1d5ba78bbbafd3226f371146bc348363
 _LaTeX_grand.xml_1.png=dd59088e24bed7a6af5a6ccd16e58616
 _LaTeX_grand.xml_2.png=4065036eed5d60beaa7f246c013cbff0
 _LaTeX_hank.xml_1.png=fc6c604bc8c86af20a8f0673047332db
-_LaTeX_histc.xml_1.png=1db429279b06e231c9795e1e83601db1
+_LaTeX_histc.xml_1.png=f1c5acc5939d55326dfab4af50e13f97
 _LaTeX_histplot.xml_1.png=f1c5acc5939d55326dfab4af50e13f97
 _LaTeX_implicitplot.xml_1.png=43ca5ad9e1f094a31392f860ef481e5c
 _LaTeX_interp.xml_1.png=b99c07a3557a83033fdeedd84352b082
@@ -661,6 +661,7 @@ evans_fr_FR_1.png=91ce84e7714915348624b03889256895
 evans_ja_JP_1.png=91ce84e7714915348624b03889256895
 evans_pt_BR_1.png=91ce84e7714915348624b03889256895
 evans_ru_RU_1.png=91ce84e7714915348624b03889256895
+factorial_1.png=f1b25752b2eafe563d1ae461774334c4
 fchamp_1.png=ee67b6aca421a82b83da821e4d7121de
 fchamp_2.png=2a078aea05fbcd35df81170239376b2a
 fec_1.png=e0a7716ead7fd2e608cc0945d0c4331c
diff --git a/scilab/modules/helptools/images/factorial_1.png b/scilab/modules/helptools/images/factorial_1.png
new file mode 100644 (file)
index 0000000..082c2e1
Binary files /dev/null and b/scilab/modules/helptools/images/factorial_1.png differ