* Bug 16617 fixed: gamma() extended to incomplete cases 52/21652/4
Samuel GOUGEON [Mon, 21 Dec 2020 23:09:00 +0000 (00:09 +0100)]
  http://bugzilla.scilab.org/16617

  gamma page (en_US PDF): http://bugzilla.scilab.org/attachment.cgi?id=5200

Change-Id: Ied34fae9a75f0bdc5cea60906bdba586feb7365f

21 files changed:
scilab/CHANGES.md
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/_LaTeX_gamma.xml_1.png [new file with mode: 0644]
scilab/modules/helptools/images/_LaTeX_gamma.xml_2.png [new file with mode: 0644]
scilab/modules/helptools/images/_LaTeX_gamma.xml_3.png [new file with mode: 0644]
scilab/modules/helptools/images/_LaTeX_gamma.xml_4.png [new file with mode: 0644]
scilab/modules/helptools/images/_LaTeX_gamma.xml_5.png [new file with mode: 0644]
scilab/modules/helptools/images/_LaTeX_gamma.xml_6.png [new file with mode: 0644]
scilab/modules/helptools/images/gamma_1.png
scilab/modules/helptools/images/gamma_2.png [new file with mode: 0644]
scilab/modules/helptools/images/gamma_3.png [new file with mode: 0644]
scilab/modules/helptools/images/gamma_equation1.mml.png [deleted file]
scilab/modules/special_functions/help/en_US/gamma.xml
scilab/modules/special_functions/help/ja_JP/gamma.xml [deleted file]
scilab/modules/special_functions/help/mml/gamma_equation1.mml [deleted file]
scilab/modules/special_functions/help/pt_BR/gamma.xml [deleted file]
scilab/modules/special_functions/macros/%s_gamma.sci
scilab/modules/special_functions/sci_gateway/cpp/sci_gamma.cpp
scilab/modules/special_functions/tests/unit_tests/gamma.tst [new file with mode: 0644]
scilab/modules/statistics/tests/unit_tests/cdfgam.dia.ref [deleted file]
scilab/modules/statistics/tests/unit_tests/cdfgam.tst

index 8d2a5cf..6c5623e 100644 (file)
@@ -201,7 +201,6 @@ Feature changes and additions on 6.1.1
   - It can now sort any sparse 2D matrix, in all `g, r, c, lr, lc` methods, including sparse booleans and in multi-level mode. It was formerly limited to sparse real or complex vectors and only to the `g` mode.
   - Any hypermatrix can be sorted along a dimension > 2.
 * `unique` is enabled for any 2D sparse arrays, in simple, 'c' and 'r' modes.
-<<<<<<< HEAD
 * `%chars` constant added, to easily access to some selected sets of unicode symbols.
 * Lists are displayed in a more compact and comprehensive way.
 * `interp1` is upgraded:
@@ -221,9 +220,9 @@ Feature changes and additions on 6.1.1
 * `setdiff` now supports input booleans and sparse matrices (boolean or numeric).
 * `intersect` is extended to any sparse boolean or numeric arrays, in all simple, 'c' or 'r' modes.
 * `union` now support boolean, sparse boolean, and sparse numerical matrices.
-=======
 * `kron` and `.*.` are extended to boolean arrays.
->>>>>>> 0e05946c3d7 (* Bug 16613 fixed: .*. and kron extended to booleans)
+* `gamma` is extended to incomplete gamma integrals.
+
 
 Help pages:
 -----------
@@ -413,13 +412,13 @@ Bug Fixes
 * [#16609](https://bugzilla.scilab.org/16609): `bitcmp` needed to be upgraded for Scilab 6.
 * [#16612](https://bugzilla.scilab.org/16612): Unlike the `.*.` operator, `kron()` was not defined for sparse numeric matrices.
 * [#16613](https://bugzilla.scilab.org/16613): Neither `.*.` nor `kron()` worked with boolean or sparse boolean matrices.
+* [#16617](https://bugzilla.scilab.org/16617): `gamma` did not support incomplete gamma integrals.
 * [#16622](https://bugzilla.scilab.org/16622): `inv` could no longer be overloaded for hypermatrices of decimal or complex numbers.
 * [#16623](https://bugzilla.scilab.org/16623): `rand(2,2,2)^2` yielded a wrong result instead of trying to call the `%s_p_s` overload for input hypermatrices.
 * [#16629](https://bugzilla.scilab.org/16629): `interp1`'s documentation did not tell the spline edges conditions ; extrapolation modes were poorly explained. ; the description of the result's size was completely wrong ; x as an option was not documented. A wrong extrapolation value could silently return a wrong result. There was some dead code like `if varargin(5)==%nan`. A bugged error message yielded its own error. When x is implicit, the argument index in error messages could be wrong. `periodic` and `edgevalue` extrapolation modes were not available. `linear` extrapolation was not available for splines. When `xp` is an hypermatrix with `size(xp,1)==1`, the size of the result was irregular/wrong.
 * [#16644](https://bugzilla.scilab.org/16644): `input("message:")` yielded a wrong error message about `mprintf` in case of non-interpretable input.
 * [#16654](https://bugzilla.scilab.org/16654): `interp` was leaking memory.
 
-
 ### Bugs fixed in 6.1.0:
 * [#2694](https://bugzilla.scilab.org/2694): `bitget` did not accept positive integers of types int8, int16 or int32.
 * [#5824](https://bugzilla.scilab.org/5824): The `datafit` algorithm was not documented.
index 63a1a8e..c26157f 100644 (file)
@@ -303,6 +303,12 @@ _LaTeX_fft.xml_1.png=b16d8d250dade64a45e1e557c4673442
 _LaTeX_fft.xml_2.png=fa418662faf728120edcfff2a838bc8e
 _LaTeX_filter.xml_1.png=f5eab130f2e7fd10b80a402952421a05
 _LaTeX_floor.xml_ru_RU_1.png=1d5ba78bbbafd3226f371146bc348363
+_LaTeX_gamma.xml_1.png=960350f434595f88477814ea577e1359
+_LaTeX_gamma.xml_2.png=3c98a6fd836c18394e05939979158ebf
+_LaTeX_gamma.xml_3.png=7aeff606c2524edc54dd085a53d673af
+_LaTeX_gamma.xml_4.png=92502c47ca464eaea97178b3a176595d
+_LaTeX_gamma.xml_5.png=073a5bfe2552bdd14720b09e0ba75ff3
+_LaTeX_gamma.xml_6.png=328f28d7f22c29a647b1b93c6843b443
 _LaTeX_grand.xml_1.png=dd59088e24bed7a6af5a6ccd16e58616
 _LaTeX_grand.xml_2.png=4065036eed5d60beaa7f246c013cbff0
 _LaTeX_hank.xml_1.png=fc6c604bc8c86af20a8f0673047332db
@@ -715,7 +721,9 @@ gainplot_fr_FR_1.png=d213684f7b1a3ed02b63da2b67d4ac33
 gainplot_ja_JP_1.png=d213684f7b1a3ed02b63da2b67d4ac33
 gainplot_pt_BR_1.png=d213684f7b1a3ed02b63da2b67d4ac33
 gainplot_ru_RU_1.png=d213684f7b1a3ed02b63da2b67d4ac33
-gamma_1.png=edb640abc2cbb2d381e6e59a3add4992
+gamma_1.png=f394755ab046e37f0c580b182a669f46
+gamma_2.png=1ada3c5257f5a417f2b68ae90a54bd30
+gamma_3.png=f7217206a73638c2fe8e5cfd96fb6cee
 genfac3d_1.png=82ced89cadf20315ddd7c7d026c81443
 geom3d_1.png=4bd6a5c64616f25637e575edfe070cd9
 geom3d_2.png=0896ca44de9f622ef3a366a18b4bea53
diff --git a/scilab/modules/helptools/images/_LaTeX_gamma.xml_1.png b/scilab/modules/helptools/images/_LaTeX_gamma.xml_1.png
new file mode 100644 (file)
index 0000000..cfd6439
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_gamma.xml_1.png differ
diff --git a/scilab/modules/helptools/images/_LaTeX_gamma.xml_2.png b/scilab/modules/helptools/images/_LaTeX_gamma.xml_2.png
new file mode 100644 (file)
index 0000000..f584ed6
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_gamma.xml_2.png differ
diff --git a/scilab/modules/helptools/images/_LaTeX_gamma.xml_3.png b/scilab/modules/helptools/images/_LaTeX_gamma.xml_3.png
new file mode 100644 (file)
index 0000000..bcf5ebe
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_gamma.xml_3.png differ
diff --git a/scilab/modules/helptools/images/_LaTeX_gamma.xml_4.png b/scilab/modules/helptools/images/_LaTeX_gamma.xml_4.png
new file mode 100644 (file)
index 0000000..4ce1b06
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_gamma.xml_4.png differ
diff --git a/scilab/modules/helptools/images/_LaTeX_gamma.xml_5.png b/scilab/modules/helptools/images/_LaTeX_gamma.xml_5.png
new file mode 100644 (file)
index 0000000..f437f7e
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_gamma.xml_5.png differ
diff --git a/scilab/modules/helptools/images/_LaTeX_gamma.xml_6.png b/scilab/modules/helptools/images/_LaTeX_gamma.xml_6.png
new file mode 100644 (file)
index 0000000..018d241
Binary files /dev/null and b/scilab/modules/helptools/images/_LaTeX_gamma.xml_6.png differ
index f8ecd29..7153a0b 100644 (file)
Binary files a/scilab/modules/helptools/images/gamma_1.png and b/scilab/modules/helptools/images/gamma_1.png differ
diff --git a/scilab/modules/helptools/images/gamma_2.png b/scilab/modules/helptools/images/gamma_2.png
new file mode 100644 (file)
index 0000000..297cc8d
Binary files /dev/null and b/scilab/modules/helptools/images/gamma_2.png differ
diff --git a/scilab/modules/helptools/images/gamma_3.png b/scilab/modules/helptools/images/gamma_3.png
new file mode 100644 (file)
index 0000000..43555cf
Binary files /dev/null and b/scilab/modules/helptools/images/gamma_3.png differ
diff --git a/scilab/modules/helptools/images/gamma_equation1.mml.png b/scilab/modules/helptools/images/gamma_equation1.mml.png
deleted file mode 100644 (file)
index 8f0b231..0000000
Binary files a/scilab/modules/helptools/images/gamma_equation1.mml.png and /dev/null differ
index 56f68f3..0865793 100644 (file)
@@ -2,8 +2,8 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
- *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ * Copyright (C) 2020 - 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: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="gamma" 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:ns4="http://www.w3.org/1999/xhtml"
+          xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"
+          xml:id="gamma" xml:lang="en">
     <refnamediv>
         <refname>gamma</refname>
-        <refpurpose>The gamma function.</refpurpose>
+        <refpurpose>gamma function, complete or incomplete normalized</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
-        <synopsis>y = gamma(x)</synopsis>
+        <synopsis>
+            y = gamma(u)
+            y = gamma(x, a)
+            y = gamma(x, a, b)
+            y = gamma(x, .., "upper")
+        </synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Arguments</title>
         <variablelist>
             <varlistentry>
-                <term>x</term>
+                <term>u</term>
                 <listitem>
-                    <para>
-                       scalar, vector, matrix, or hypermatrix of real numbers.
-                    </para>
-                    <note>
-                      <literal>gamma</literal> can be overloaded for complex numbers or
-                      of lists, tlists or mlists.
-                    </note>
+                    array of positive or negative real numbers
+                    <para/>
+                    <literal>gamma(u)</literal> and <literal>gamma(x,…)</literal> can be overloaded
+                    for complex numbers with <literal>%s_gamma_user()</literal>, and for other
+                    <varname>a</varname> types with the usual overload naming rule.
+                    <para/>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>x, a, b</term>
+                <listitem>
+                    arrays of positive real numbers.
+                    If at least one input is not scalar, scalar ones are expanded to its size.
+                    If several inputs are not scalar, they must have the same size.
+                    <para/>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>y</term>
                 <listitem>
-                    <para>real vector or matrix with same size than x.</para>
+                    array of real numbers, with the size of <varname>u</varname>
+                    or of (the non-scalar) <varname>x</varname>, <varname>a</varname>,
+                    or <varname>b</varname>.
+                    <para/>
                 </listitem>
             </varlistentry>
         </variablelist>
     <refsection>
         <title>Description</title>
         <para>
-            <literal>gamma(x)</literal> evaluates the gamma function at all the
-            elements of <literal>x</literal>. The gamma function is defined by
-            :
+            <literal>gamma(…)</literal> computes and yields the complete or incomplete gamma
+            function for each element of its input(s), in an element-wise way. The complete
+            gamma function extends the factorial one to non-integer real positive or negative
+            numbers, as <literal>gamma(u+1)=u*gamma(u)</literal>.
         </para>
-        <informalequation>
-            <mediaobject>
-                <imageobject>
-                    <imagedata align="center" fileref="../mml/gamma_equation1.mml"/>
-                </imageobject>
-            </mediaobject>
-        </informalequation>
-        <para>and generalizes the factorial function for real numbers
-            (<literal>gamma(u+1) = u*gamma(u)</literal>).
+        <para>
+            <emphasis role="bold">gamma(u)</emphasis> computes
+            <latex style="display" fontsize="18" alt="Γ(u)= ∫_0→∞ t^{u-1}.exp(-t).dt">
+                \Gamma(u)=\int_0^\infty\! t^{u-1}e^{-t}\,dt
+            </latex>
         </para>
-    </refsection>
+        <refsect3>
+            <title>Incomplete normalized integrals</title>
+            <para>
+                <emphasis role="bold">gamma(x, a)</emphasis> computes the integral
+                <latex style="display" fontsize="18" alt="P(x,a)= ∫_0→x t^{a-1}.exp(-t).dt / Γ(a)">
+                    P(x,a)=\frac{1}{\Gamma(a)}\int_0^x\! t^{a-1}e^{-t}\,dt
+                </latex>
+            </para>
+            <para>
+                <emphasis role="bold">gamma(x, a, b)</emphasis> computes the generalized integral
+                <latex style="display" fontsize="18" alt="P(x,a,b)= ∫_0→x t^{a-1}.exp(-bt).dt . b^a / Γ(a)">
+                    P(x,a,b)=\frac{b^a}{\Gamma(a)}\int_0^x\! t^{a-1}e^{-b\,t}\,dt
+                </latex>
+            </para>
+            <para>
+                <emphasis role="bold">gamma(x, a, "upper")</emphasis> computes accurately the
+                complementary integral
+                <latex style="display" fontsize="18" alt="Q(x,a)= ∫_x→∞ t^{a-1}.exp(-t).dt / Γ(a) = 1-P(x,a)">
+                    Q(x,a)=1-P(x,a)=\frac{1}{\Gamma(a)}\int_x^\infty\! t^{a-1}e^{-t}\,dt
+                </latex>
+                even for big x and P(x,a)→1. Finally,
+            </para>
+            <para>
+                <emphasis role="bold">gamma(x, a, b, "upper")</emphasis> computes the generalized
+                complementary integral
+                <latex style="display" fontsize="18" alt="Q(x,a,b)= ∫_x→∞ t^{a-1}.exp(-bt).dt . b^a / Γ(a)">
+                    Q(x,a,b)=\frac{b^a}{\Gamma(a)}\int_x^\infty\! t^{a-1}e^{-b\,t}\,dt
+                </latex>
+            </para>
+            <note>
+                <para>
+                    The inverse incomplete normalized gamma function can be computed with
+                    <emphasis role="bold">x = cdfgam("X", a, b, y, 1-y)</emphasis>,
+                    that is the <varname>x</varname> bound such that
+                    <latex alt="y=∫_0→x t^{a-1}.exp(-bt).dt . b^a / Γ(a)">
+                    y=\frac{b^a}{\Gamma(a)}\int_0^x\! t^{a-1}e^{-b\,t}\,dt
+                    </latex>
+                </para>
+                <para>
+                    Calling <emphasis role="bold">x = cdfgam("X", a, b, z-1, z)</emphasis> with
+                    <literal>z=1-y</literal> will be preferred when 0.5 &lt; y &lt; 1, to get a full
+                    accuracy on <varname>x</varname>.
+                </para>
+            </note>
+        </refsect3>
+   </refsection>
     <refsection>
         <title>Examples</title>
+        <para>
+            Gamma as the extension of the factorial function to non-integer numbers:
+        </para>
         <programlisting role="example"><![CDATA[
-// simple examples
-gamma(0.5)
-gamma(6)-prod(1:5)
- ]]></programlisting>
+[gamma(2:7) ; factorial(1:6)]
+gamma(1.5:7)
+gamma(1.5:7) ./ gamma(0.5:6)
+     ]]></programlisting>
+        <screen><![CDATA[
+--> [gamma(2:7) ; factorial(1:6)]
+ ans  =
+   1.   2.   6.   24.   120.   720.
+   1.   2.   6.   24.   120.   720.
+
+--> gamma(1.5:7)
+ ans  =
+   0.8862269   1.3293404   3.323351   11.631728   52.342778   287.88528
+
+--> gamma(1.5:7) ./ gamma(0.5:6)
+ ans  =
+   0.5   1.5   2.5   3.5   4.5   5.5
+]]></screen>
+        <para/>
+        <para>
+            Graph of the Gamma function around 0:
+        </para>
         <programlisting role="example"><![CDATA[
-// the graph of the Gamma function on [a,b]
-a = -3; b = 5;
+[a, b] = (-3, 5);
 x = linspace(a,b,40000);
 y = gamma(x);
-clf()
-plot2d(x, y, style=0, axesflag=5, rect=[a, -10, b, 10])
-xtitle("The gamma function on ["+string(a)+","+string(b)+"]")
-show_window() ]]></programlisting>
+clf
+plot2d(x, y, style=0, axesflag=5, rect=[a,-10,b,10])
+title("$\Gamma(u)$", "fontsize",3.5)
+xgrid(color("grey60"))
+     ]]></programlisting>
         <scilab:image>
-            a = -3; b = 5;
-            x = linspace(a,b,40000)';
+            [a, b] = (-3, 5);
+            x = linspace(a,b,40000);
             y = gamma(x);
-            plot2d(x, y, style=0, axesflag=5, rect=[a, -10, b, 10])
-            xtitle("The gamma function on ["+string(a)+","+string(b)+"]")
+            clf
+            plot2d(x, y, style=0, axesflag=5, rect=[a,-10,b,10])
+            title("$\Gamma(u)$", "fontsize",3.5)
+            xgrid(color("grey60"))
+            gcf().axes_size = [550,400];
+        </scilab:image>
+        <para/>
+        <para>
+            Incomplete normalized P(x,a) gamma function:
+        </para>
+        <programlisting role="example"><![CDATA[
+x = 0.1:0.2:8;
+a = 0.1:0.2:7;
+[X, A] = ndgrid(x, a);
+P = gamma(X,A);
+clf
+gcf().color_map = coolcolormap(100);
+surf(a,x,P)
+title("$P(x,a)=\frac{1}{\Gamma(a)}\int_0^x\! t^{a-1}e^{-t}\,dt$","fontsize",3.5)
+xlabel(["" "a"], "fontsize",2)
+ylabel("x", "fontsize",2)
+zlabel("P(x,a)", "fontsize",2)
+xgrid
+     ]]></programlisting>
+        <scilab:image>
+            x = 0.1:0.2:8;
+            a = 0.1:0.2:7;
+            [X, A] = ndgrid(x, a);
+            P = gamma(X,A);
+            clf
+            gcf().color_map = coolcolormap(100);
+            surf(a,x,P)
+            title("$P(x,a)=\frac{1}{\Gamma(a)}\int_0^x\! t^{a-1}e^{-t}\,dt$","fontsize",3.5)
+            xlabel(["" "a"], "fontsize",2)
+            ylabel("x", "fontsize",2)
+            zlabel("P(x,a)", "fontsize",2)
+            xgrid
+            gca().rotation_angles = [64 18];
+        </scilab:image>
+        <para/>
+        <para>
+            Incomplete generalized normalized P(x,a,b) function:
+        </para>
+        <programlisting role="example"><![CDATA[
+a = 0.1:0.2:8;
+b = 0.1:0.2:7;
+[A, B] = ndgrid(a, b);
+P = gamma(1,A,B);
+clf
+gcf().color_map = parulacolormap(100);
+surf(b,a,P)
+title("$P(x,a,b)=\frac{b^a}{\Gamma(a)}\int_0^x\! t^{a-1}e^{-b\,t}\,dt\quad for\quad x=1$","fontsize",3.7)
+xlabel("b", "fontsize",2)
+ylabel("a", "fontsize",2)
+zlabel("")
+gca().rotation_angles = [58 75];
+xgrid
+     ]]></programlisting>
+        <scilab:image>
+            a = 0.1:0.2:8;
+            b = 0.1:0.2:7;
+            [A, B] = ndgrid(a, b);
+            P = gamma(1,A,B);
+            clf
+            gcf().color_map = parulacolormap(100);
+            surf(b,a,P)
+            title("$P(x,a,b)=\frac{b^a}{\Gamma(a)}\int_0^x\! t^{a-1}e^{-b\,t}\,dt\quad for\quad x=1$","fontsize",3.7)
+            xlabel("b", "fontsize",2)
+            ylabel("a", "fontsize",2)
+            zlabel("")
+            gca().rotation_angles = [58 75];
+            xgrid
         </scilab:image>
     </refsection>
     <refsection role="see also">
@@ -97,6 +249,9 @@ show_window() ]]></programlisting>
                 <link linkend="dlgamma">dlgamma</link>
             </member>
             <member>
+                <link linkend="cdfgam">cdfgam</link>
+            </member>
+            <member>
                 <link linkend="factorial">factorial</link>
             </member>
         </simplelist>
@@ -121,6 +276,10 @@ show_window() ]]></programlisting>
                 </itemizedlist>
               </revremark>
             </revision>
+            <revision>
+                <revnumber>6.1.1</revnumber>
+                <revremark>gamma(x,..) incomplete versions added.</revremark>
+            </revision>
         </revhistory>
     </refsection>
 </refentry>
diff --git a/scilab/modules/special_functions/help/ja_JP/gamma.xml b/scilab/modules/special_functions/help/ja_JP/gamma.xml
deleted file mode 100644 (file)
index 44c7ce2..0000000
+++ /dev/null
@@ -1,130 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA
- *
- * 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: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="gamma" xml:lang="ja">
-    <refnamediv>
-        <refname>gamma</refname>
-        <refpurpose>ガンマ関数.</refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>呼び出し手順</title>
-        <synopsis>y = gamma(x)</synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>引数</title>
-        <variablelist>
-            <varlistentry>
-                <term>x</term>
-                <listitem>
-                    <para>
-                       scalar, vector, matrix, or hypermatrix of real numbers.
-                    </para>
-                    <note>
-                      <literal>gamma</literal> can be overloaded for complex numbers or
-                      of lists, tlists or mlists.
-                    </note>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>y</term>
-                <listitem>
-                    <para>(xと同じ大きさの)実数ベクトルまたは行列.</para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>説明</title>
-        <para>
-            <literal>gamma(x)</literal> は,
-            <literal>x</literal>の全要素についてガンマ関数を計算します.
-            ガンマ関数は以下のように定義されます:
-        </para>
-        <informalequation>
-            <mediaobject>
-                <imageobject>
-                    <imagedata align="center" fileref="../mml/gamma_equation1.mml"/>
-                </imageobject>
-            </mediaobject>
-        </informalequation>
-        <para>そして,階乗関数を実数に一般化します.
-            (<literal>gamma(u+1) = u*gamma(u)</literal>).
-        </para>
-    </refsection>
-    <refsection>
-        <title>例</title>
-        <programlisting role="example"><![CDATA[
-// 簡単な例
-gamma(0.5)
-gamma(6)-prod(1:5)
- ]]></programlisting>
-        <programlisting role="example"><![CDATA[
-// [a,b]のガンマ関数のグラフ
-a = -3; b = 5;
-x = linspace(a,b,40000);
-y = gamma(x);
-clf()
-plot2d(x, y, style=0, axesflag=5, rect=[a, -10, b, 10])
-xtitle("The gamma function on ["+string(a)+","+string(b)+"]")
-show_window()
- ]]></programlisting>
-        <scilab:image>
-            a = -3; b = 5;
-            x = linspace(a,b,40000)';
-            y = gamma(x);
-            plot2d(x, y, style=0, axesflag=5, rect=[a, -10, b, 10])
-            xtitle("The gamma function on ["+string(a)+","+string(b)+"]")
-        </scilab:image>
-    </refsection>
-    <refsection role="see also">
-        <title>参照</title>
-        <simplelist type="inline">
-            <member>
-                <link linkend="gammaln">gammaln</link>
-            </member>
-            <member>
-                <link linkend="dlgamma">dlgamma</link>
-            </member>
-            <member>
-                <link linkend="factorial">factorial</link>
-            </member>
-        </simplelist>
-    </refsection>
-    <refsection>
-        <title>履歴</title>
-        <revhistory>
-            <revision>
-                <revnumber>5.4.0</revnumber>
-                <revremark>
-                    list, mlist, tlistおよびハイパー行列型のオーバーロードが
-                    可能となりました.
-                </revremark>
-            </revision>
-            <revision>
-              <revnumber>6.0.2</revnumber>
-              <revremark>
-                <itemizedlist>
-                  <listitem>
-                    The input can now be an hypermatrix.
-                  </listitem>
-                  <listitem>
-                    <literal>gamma</literal> can now be overloaded for complex numbers.
-                  </listitem>
-                </itemizedlist>
-              </revremark>
-            </revision>
-        </revhistory>
-    </refsection>
-</refentry>
diff --git a/scilab/modules/special_functions/help/mml/gamma_equation1.mml b/scilab/modules/special_functions/help/mml/gamma_equation1.mml
deleted file mode 100644 (file)
index 85df255..0000000
+++ /dev/null
@@ -1,50 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!DOCTYPE math:math PUBLIC "-//OpenOffice.org//DTD Modified W3C MathML 1.01//EN" "math.dtd">
-<math:math xmlns:math="http://www.w3.org/1998/Math/MathML">
- <math:semantics>
-  <math:mrow>
-   <math:mo math:stretchy="false">Γ</math:mo>
-   <math:mrow>
-    <math:mrow>
-     <math:mo math:stretchy="false">(</math:mo>
-     <math:mi>x</math:mi>
-     <math:mo math:stretchy="false">)</math:mo>
-    </math:mrow>
-    <math:mo math:stretchy="false">=</math:mo>
-    <math:mrow>
-     <math:mrow>
-      <math:mrow>
-       <math:munderover>
-        <math:mo math:stretchy="false">∫</math:mo>
-        <math:mn>0</math:mn>
-        <math:mrow>
-         <math:mo math:stretchy="false">+</math:mo>
-         <math:mo math:stretchy="false">∞</math:mo>
-        </math:mrow>
-       </math:munderover>
-       <math:msup>
-        <math:mi>t</math:mi>
-        <math:mrow>
-         <math:mi>x</math:mi>
-         <math:mo math:stretchy="false">−</math:mo>
-         <math:mn>1</math:mn>
-        </math:mrow>
-       </math:msup>
-      </math:mrow>
-      <math:mo math:stretchy="false">⋅</math:mo>
-      <math:msup>
-       <math:mi>e</math:mi>
-       <math:mrow>
-        <math:mo math:stretchy="false">−</math:mo>
-        <math:mi>t</math:mi>
-       </math:mrow>
-      </math:msup>
-     </math:mrow>
-     <math:mo math:stretchy="false">⋅</math:mo>
-     <math:mi math:fontstyle="italic">dt</math:mi>
-    </math:mrow>
-   </math:mrow>
-  </math:mrow>
-  <math:annotation math:encoding="StarMath 5.0">%GAMMA(x) = int from 0 to {+infty} t^{x-1} cdot e^{-t} cdot dt</math:annotation>
- </math:semantics>
-</math:math>
\ No newline at end of file
diff --git a/scilab/modules/special_functions/help/pt_BR/gamma.xml b/scilab/modules/special_functions/help/pt_BR/gamma.xml
deleted file mode 100644 (file)
index 839a1c0..0000000
+++ /dev/null
@@ -1,118 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA
- *
- * 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:ns5="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="gamma" xml:lang="pt">
-    <refnamediv>
-        <refname>gamma</refname>
-        <refpurpose>função gama </refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>Seqüência de Chamamento</title>
-        <synopsis>y = gamma(x)</synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>Parâmetros</title>
-        <variablelist>
-            <varlistentry>
-                <term>x</term>
-                <listitem>
-                    <para>
-                       scalar, vector, matrix, or hypermatrix of real numbers.
-                    </para>
-                    <note>
-                      <literal>gamma</literal> can be overloaded for complex numbers or
-                      of lists, tlists or mlists.
-                    </note>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>y</term>
-                <listitem>
-                    <para>vetor ou matriz de reais ou complexos de mesmo tamanho que
-                        x
-                    </para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>Descrição</title>
-        <para>
-            <literal>gamma(x)</literal> avalia a função gama em todos os
-            elementos de <literal>x</literal>. A função gama é defininda por :
-        </para>
-        <programlisting role=""><![CDATA[
-            /+oo
-            |   x-1  -t
-gamma(x) =  |  t    e    dt
-            /0
- ]]></programlisting>
-        <para>e generaliza a função fatorial para os números reais
-            (<literal>gamma(u+1) = u*gamma(u)</literal>).
-        </para>
-    </refsection>
-    <refsection>
-        <title>Exemplos</title>
-        <programlisting role="example"><![CDATA[
-// exemplos simples
-gamma(0.5)
-gamma(6)-prod(1:5)
-
-// o gráfico da função gama em [a,b]
-a = -3; b = 5;
-x = linspace(a,b,40000);
-y = gamma(x);
-clf()
-plot2d(x, y, style=0, axesflag=5, rect=[a, -10, b, 10])
-xtitle("A função gama em  ["+string(a)+","+string(b)+"]")
-show_window() ]]></programlisting>
-    </refsection>
-    <refsection role="see also">
-        <title>Ver Também</title>
-        <simplelist type="inline">
-            <member>
-                <link linkend="gammaln">gammaln</link>
-            </member>
-            <member>
-                <link linkend="dlgamma">dlgamma</link>
-            </member>
-            <member>
-                <link linkend="factorial">factorial</link>
-            </member>
-        </simplelist>
-    </refsection>
-    <refsection>
-        <title>Histórico</title>
-        <revhistory>
-            <revision>
-                <revnumber>5.4.0</revnumber>
-                <revremark>Overloading allowed for list, mlist, tlist and hypermatrix types.</revremark>
-            </revision>
-            <revision>
-              <revnumber>6.0.2</revnumber>
-              <revremark>
-                <itemizedlist>
-                  <listitem>
-                    The input can now be an hypermatrix.
-                  </listitem>
-                  <listitem>
-                    <literal>gamma</literal> can now be overloaded for complex numbers.
-                  </listitem>
-                </itemizedlist>
-              </revremark>
-            </revision>
-        </revhistory>
-    </refsection>
-</refentry>
index e3b6a23..3891801 100644 (file)
@@ -1,5 +1,5 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2018 - Samuel GOUGEON
+// Copyright (C) 2018 - 2020 - 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.
 // For more information, see the COPYING file which you should have received
 // along with this program.
 
-function y = %s_gamma(x)
+function y = %s_gamma(varargin)
     // call for complex numbers or/and for hypermatrix of reals
-    s = size(x);
-    x = x(:);
-    if isreal(x, 0)
-        y = gamma(real(x));
-    else
+    //      or/and for incomplete integral
+    // gamma(a)          // with a complex or real hypermat
+    // gamma(x, a)       // with real x >= 0
+    // gamma(x, a, b)
+    // gamma(x,.., option)  // "upper": complementary integral
+
+    msg = _("%s: Function not defined for the given argument type.\n  Check arguments or define function %s for overloading.\n")
+    rhs = argn(2)
+    a = varargin(min(2,rhs))
+    sa = size(a)
+    y = []
+
+    // gamma(..)  with complex numbers (all cases: complete, uncomplete, hypermat)
+    // -------------------------------
+    if or(type(a)==[1 5]) & ~isreal(a,0)
         if isdef("%s_gamma_user") & type(%s_gamma_user)==13
-            y = %s_gamma_user(x);
+            y = %s_gamma_user(varargin(:));
+            // Note: the overload must be able to process hypermatrices
+            return
         else
-            msg = _("%s: Function not defined for the given argument type.\n  Check arguments or define function %s for overloading.\n")
             error(msprintf(msg, "%s_gamma", "%s_gamma_user()"))
         end
     end
-    y = matrix(y, s);
+    if a==[] then
+        return
+    end
+    a = real(a)
+
+    // gamma(a) with hypermatrix of real numbers
+    // -----------------------------------------
+    if rhs==1 then
+        y = matrix(gamma(a(:)), sa)
+        return
+    end
+
+    // gamma(x,a,..) : uncomplete gamma (all cases with real numbers)
+    // -------------
+    x = varargin(1)
+    if type(x)<>1 | ~isreal(x,0) then
+        msg = _("%s: Argument #%d: Decimal numbers expected.\n")
+        error(msprintf(msg, "gamma", 1))
+    end
+    if x==[] then
+        return
+    end
+    x = real(x)
+    if min(x) <= 0 then
+        msg = _("%s: Argument #%d: Must be > %d.\n")
+        error(msprintf(msg, "gamma", 1, 0))
+    end
+    if ~isscalar(x) & ~isscalar(a) & or(size(x)<>size(a)) then
+        msg = _("%s: Arguments #%d and #%d: Same sizes expected.\n")
+        error(msprintf(msg, "gamma", 1, 2))
+    end
+    if isscalar(a) then
+        a = a * ones(x)
+    elseif isscalar(x)
+        x = x * ones(a)
+    end
+    // "upper"
+    upper = %f
+    u = varargin($)
+    if type(u)==10
+        if convstr(u(1))<>"upper" then
+            msg = _("%s: Argument #%d: ''%s'' expected .\n")
+            error(msprintf(msg, "gamma", rhs, "upper"))
+        end
+        upper = %t
+        varargin($) = null()
+    end
+    // b
+    if length(varargin) > 2 then
+        b = varargin(3)
+        if type(b)<>1 | ~isreal(b,0) then
+            msg = _("%s: Argument #%d: Decimal numbers expected.\n")
+            error(msprintf(msg, "gamma", 3))
+        end
+        b = real(b)
+        if min(b) <= 0 then
+            msg = _("%s: Argument #%d: Must be > %d.\n")
+            error(msprintf(msg, "gamma", 3, 0))
+        end
+        if ~isscalar(a) & ~isscalar(b) & or(size(a)<>size(b)) then
+            msg = _("%s: Arguments #%d and #%d: Same sizes expected.\n")
+            error(msprintf(msg, "gamma", 2, 3))
+        end
+    else
+        b = 1
+    end
+    if b==[] then
+        return
+    end
+    if isscalar(b) then
+        b = b * ones(a)
+    end
+
+    // PROCESSING
+    // ==========
+    n = ndims(x)
+    sa = size(a)
+    if n > 2 then
+        [x, a, b] = (x(:), a(:), b(:))
+    end
+    z = find(a==0)
+    a(z) = 1
+
+    if upper then
+        [?,y] = cdfgam("PQ", x, a, b)
+    else
+        y = cdfgam("PQ", x, a, b)
+    end
+
+    if z<>[] then
+        y(z) = 1
+    end
+    if n > 2 then
+        y = matrix(y, sa)
+    end
 endfunction
index c2144dc..8690a3d 100644 (file)
@@ -35,24 +35,12 @@ types::Function::ReturnValue sci_gamma(types::typed_list &in, int _iRetCount, ty
         return types::Function::Error;
     }
 
-    if (in[0]->isList() || in[0]->isTList() || in[0]->isMList())
+    if (in.size() > 1 || in[0]->isDouble() == false)
     {
         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_gamma";
         return Overload::call(wstFuncName, in, _iRetCount, out);
     }
 
-    if (in.size() > 1)
-    {
-        Scierror(77, _("%s: Wrong number of input argument: %d expected.\n"), "gamma", 1);
-        return types::Function::Error;
-    }
-
-    if (in[0]->isDouble() == false)
-    {
-        Scierror(999, _("%s: Argument #%d: Decimal number(s) expected.\n"), "gamma", 1);
-        return types::Function::Error;
-    }
-
     /***** get data *****/
     types::Double* pDblIn = in[0]->getAs<types::Double>();
 
diff --git a/scilab/modules/special_functions/tests/unit_tests/gamma.tst b/scilab/modules/special_functions/tests/unit_tests/gamma.tst
new file mode 100644 (file)
index 0000000..c554c0e
--- /dev/null
@@ -0,0 +1,119 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2020 - Samuel GOUGEON
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+// <-- Short Description -->
+// Unitary tests of the gamma() function
+//
+
+// Gamma complete
+// --------------
+assert_checkequal(gamma([]), []);
+
+M = [0.5 1.5 ; 2.5 3.5];
+H = matrix(M,2,1,2);        // 3D hypermatrix
+Ref = [1.77245385090551610 0.88622692545275805      // R 4.0.3
+       1.32934038817913702 3.32335097044784256 ];
+assert_checkequal(gamma(M), Ref);
+assert_checkequal(gamma(H), matrix(Ref,2,1,2));
+
+
+// Gamma incomplete
+// ----------------
+assert_checkequal(%s_gamma(0.5,[]), []);
+assert_checkequal(%s_gamma([],0.5), []);
+// x scalar, a array
+Ref = [0.68268949213708585 0.19874804309879918      // R 4.0.3
+       0.03743422675270363 0.0051714634834845166 ];
+assert_checkalmostequal(%s_gamma(0.5, M), Ref, 20*%eps);
+assert_checkalmostequal(%s_gamma(0.5, H), matrix(Ref,2,1,2), 20*%eps);
+assert_checkequal(%s_gamma(%inf, M), ones(2,2));
+assert_checkequal(%s_gamma(%inf, H), ones(2,1,2));
+// x array, a scalar
+Ref = [0.68268949213708585 0.91673548333644994      // R 4.0.3
+       0.97465268132253169 0.99184902840649725 ];
+assert_checkalmostequal(%s_gamma(M, 0.5), Ref, 2*%eps);
+assert_checkalmostequal(%s_gamma(H, 0.5), matrix(Ref,2,1,2), 2*%eps);
+// x and a arrays
+Ref = [0.68268949213708585 0.60837482372891105      // R 4.0.3
+       0.58411981300449223 0.57112014244694542 ];
+assert_checkalmostequal(%s_gamma(M, M), Ref, 10*%eps);
+assert_checkalmostequal(%s_gamma(H, H), matrix(Ref,2,1,2), 10*%eps);
+
+
+// Gamma incomplete COMPLEMENTARY
+// ------------------------------
+assert_checkequal(%s_gamma(0.5,[],"upper"), []);
+assert_checkequal(%s_gamma([],0.5,"upper"), []);
+// x scalar, a array
+Ref = 1- [0.68268949213708585 0.19874804309879918      // R 4.0.3
+          0.03743422675270363 0.0051714634834845166 ];
+assert_checkalmostequal(%s_gamma(0.5, M, "upper"), Ref, 5*%eps);
+assert_checkalmostequal(%s_gamma(0.5, H, "upper"), matrix(Ref,2,1,2), 5*%eps);
+assert_checkequal(%s_gamma(%inf, M, "upper"), zeros(2,2));
+assert_checkequal(%s_gamma(%inf, H, "upper"), zeros(2,1,2));
+// x array, a scalar
+Ref = 1 - [0.68268949213708585 0.91673548333644994      // R 4.0.3
+           0.97465268132253169 0.99184902840649725 ];
+assert_checkalmostequal(%s_gamma(M, 0.5, "upper"), Ref, 50*%eps);
+assert_checkalmostequal(%s_gamma(H, 0.5, "upper"), matrix(Ref,2,1,2), 50*%eps);
+// x and a arrays
+Ref = 1 - [0.68268949213708585 0.60837482372891105      // R 4.0.3
+           0.58411981300449223 0.57112014244694542 ];
+assert_checkalmostequal(%s_gamma(M, M, "upper"), Ref, 10*%eps);
+assert_checkalmostequal(%s_gamma(H, H, "upper"), matrix(Ref,2,1,2), 10*%eps);
+
+
+// Gamma incomplete generalized
+// ----------------------------
+assert_checkequal(%s_gamma([],0.5,2), []);
+assert_checkequal(%s_gamma(0.5,[],2), []);
+assert_checkequal(%s_gamma(0.5,0.5,[]), []);
+// x scalar, a and b arrays
+Ref = [0.52049987781304685   0.31772966966378619
+       0.22350492887667642   0.16477451738965726 ];
+assert_checkalmostequal(%s_gamma(0.5, M, M), Ref, %eps);
+assert_checkalmostequal(%s_gamma(0.5, H, H), matrix(Ref,2,1,2), %eps);
+// x and a arrays, b scalar
+assert_checkalmostequal(%s_gamma(M, M, 0.5), Ref, %eps);
+assert_checkalmostequal(%s_gamma(H, H, 0.5), matrix(Ref,2,1,2), %eps);
+// x array, a scalar, b array
+Ref = [0.52049987781304685   0.96610514647531076
+       0.99959304798255499   0.99999925690162761 ];
+assert_checkalmostequal(%s_gamma(M, 0.5, M), Ref, %eps);
+assert_checkalmostequal(%s_gamma(H, 0.5, H), matrix(Ref,2,1,2), %eps);
+// x, a, and b arrays
+Ref = [0.52049987781304685   0.78770971263986711
+       0.9714568766738326    0.99906979057712397 ];
+assert_checkalmostequal(%s_gamma(M, M, M), Ref, %eps);
+assert_checkalmostequal(%s_gamma(H, H, H), matrix(Ref,2,1,2), %eps);
+
+
+// Gamma incomplete generalized COMPLEMENTARY
+// ------------------------------------------
+assert_checkequal(%s_gamma([],0.5,2, "upper"), []);
+assert_checkequal(%s_gamma(0.5,[],2, "upper"), []);
+assert_checkequal(%s_gamma(0.5,0.5,[], "upper"), []);
+// x scalar, a and b arrays
+Ref = 1 - [0.52049987781304685   0.31772966966378619
+           0.22350492887667642   0.16477451738965726 ];
+assert_checkalmostequal(%s_gamma(0.5, M, M, "upper"), Ref, %eps);
+assert_checkalmostequal(%s_gamma(0.5, H, H, "upper"), matrix(Ref,2,1,2), %eps);
+// x and a arrays, b scalar
+assert_checkalmostequal(%s_gamma(M, M, 0.5, "upper"), Ref, %eps);
+assert_checkalmostequal(%s_gamma(H, H, 0.5, "upper"), matrix(Ref,2,1,2), %eps);
+// x array, a scalar, b array
+Ref = [0.47950012218695315   3.389485352468927D-02
+       4.069520174449587D-04 7.430983723414120D-07 ];
+assert_checkalmostequal(%s_gamma(M, 0.5, M, "upper"), Ref, %eps);
+assert_checkalmostequal(%s_gamma(H, 0.5, H, "upper"), matrix(Ref,2,1,2), %eps);
+// x, a, and b arrays
+Ref = [0.47950012218695315    0.21229028736013292
+       2.854312332616740D-02  9.302094228759846D-04 ];
+assert_checkalmostequal(%s_gamma(M, M, M, "upper"), Ref, %eps);
+assert_checkalmostequal(%s_gamma(H, H, H, "upper"), matrix(Ref,2,1,2), %eps);
diff --git a/scilab/modules/statistics/tests/unit_tests/cdfgam.dia.ref b/scilab/modules/statistics/tests/unit_tests/cdfgam.dia.ref
deleted file mode 100644 (file)
index e66bd0e..0000000
+++ /dev/null
@@ -1,127 +0,0 @@
-// =============================================================================
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2010 - DIGITEO - Michael Baudin
-// Copyright (C) 2011 - DIGITEO - Michael Baudin
-//
-//  This file is distributed under the same license as the Scilab package.
-// =============================================================================
-// <-- CLI SHELL MODE -->
-// <-- ENGLISH IMPOSED -->
-// Run with test_run('statistics','cdfgam',['no_check_error_output']);
-//
-// Assessing the quality of the Normal distribution function
-// References
-//   Yalta, A. T. 2008. The accuracy of statistical distributions in Microsoft®Excel 2007. Comput. Stat. Data Anal. 52, 10 (Jun. 2008), 4579-4586. DOI= http://dx.doi.org/10.1016/j.csda.2008.03.005
-//   Computation of Statistical Distributions (ELV), Leo Knüsel
-// Table 5
-// Check Gamma distribution with parameters (x, alpha, beta = 1, Sigma = 1)
-//
-// Table of inputs from Yalta, 2008
-// [x shape scale P ]
-table = [
-0.1 , 0.1 , 1 , 0.827552
-0.2 , 0.1 , 1 , 0.879420
-0.2 , 0.2 , 1 , 0.764435
-0.3 , 0.2 , 1 , 0.816527
-0.3 , 0.3 , 1 , 0.726957
-0.4 , 0.3 , 1 , 0.776381
-0.4 , 0.4 , 1 , 0.701441
-0.5 , 0.4 , 1 , 0.748019
-0.5 , 0.5 , 1 , 0.682689
-0.6 , 0.5 , 1 , 0.726678
-];
-precision = 1.e-5;
-ntests = size(table,"r");
-for i = 1 : ntests
-    x = table(i,1);
-    shape = table(i,2);
-    scale = table(i,3);
-    expected = table(i,4);
-    // Caution: this is the rate !
-    rate = 1/scale;
-    [computed,Q]=cdfgam("PQ",x,shape,rate);
-    assert_checkalmostequal ( computed , expected , precision );
-    assert_checkalmostequal ( Q , 1 - expected , precision );
-end
-// Table of inputs computed from R-2.8.1
-// [x shape scale PDF-P CDF-P CDF-Q]
-table = [
-1.000000000000000056D-01 1.000000000000000056D-01 1.000000000000000000D+00 7.554920138253073958D-01 8.275517595858505882D-01 1.724482404141494951D-01
-2.000000000000000111D-01 1.000000000000000056D-01 1.000000000000000000D+00 3.663307993056703071D-01 8.794196267900568076D-01 1.205803732099432063D-01
-2.000000000000000111D-01 2.000000000000000111D-01 1.000000000000000000D+00 6.462857778271943188D-01 7.644345975029189777D-01 2.355654024970809945D-01
-2.999999999999999889D-01 2.000000000000000111D-01 1.000000000000000000D+00 4.227875047602157044D-01 8.165267943336527168D-01 1.834732056663473110D-01
-2.999999999999999889D-01 2.999999999999999889D-01 1.000000000000000000D+00 5.752117576599179438D-01 7.269573437103662439D-01 2.730426562896338116D-01
-4.000000000000000222D-01 2.999999999999999889D-01 1.000000000000000000D+00 4.255407854753925911D-01 7.763805810166358734D-01 2.236194189833642099D-01
-4.000000000000000222D-01 4.000000000000000222D-01 1.000000000000000000D+00 5.236648604477927016D-01 7.014412706419403953D-01 2.985587293580597157D-01
-5.000000000000000000D-01 4.000000000000000222D-01 1.000000000000000000D+00 4.144555659263016167D-01 7.480185547260104206D-01 2.519814452739895239D-01
-5.000000000000000000D-01 5.000000000000000000D-01 1.000000000000000000D+00 4.839414490382866751D-01 6.826894921370858516D-01 3.173105078629140929D-01
-5.999999999999999778D-01 5.000000000000000000D-01 1.000000000000000000D+00 3.997355278034666060D-01 7.266783217077018575D-01 2.733216782922981980D-01
-5.000000000000000000D-01 5.000000000000000000D-01 2.000000000000000000D+00 4.393912894677223790D-01 5.204998778130465187D-01 4.795001221869534258D-01
-5.000000000000000000D-01 5.000000000000000000D-01 3.000000000000000000D+00 3.899393114454822729D-01 4.362971383492270094D-01 5.637028616507729906D-01
-5.000000000000000000D-01 5.000000000000000000D-01 4.000000000000000000D+00 3.520653267642995243D-01 3.829249225480261809D-01 6.170750774519737636D-01
-1.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 2.075537487102973866D-01 8.427007929497148941D-01 1.572992070502851891D-01
-2.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.399096651318804896D-02 9.544997361036415828D-01 4.550026389635838248D-02
-4.000000000000000000D+00 5.000000000000000000D-01 1.000000000000000000D+00 5.166746338523012412D-03 9.953222650189527121D-01 4.677734981047261889D-03
-1.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 8.099910956089122777D-06 9.999922557835689840D-01 7.744216431044085842D-06
-2.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 2.600281868827196957D-10 9.999999997460371493D-01 2.539628589470869077D-10
-4.000000000000000000D+01 5.000000000000000000D-01 1.000000000000000000D+00 3.789795640412981196D-19 1.000000000000000000D+00 3.744097384202895045D-19
-3.000000000000000000D+02 5.000000000000000000D-01 1.000000000000000000D+00 1.67694904029982009D-132 1.000000000000000000D+00 1.67416798469182012D-132
-1.000000000000000021D-02 5.000000000000000000D-01 1.000000000000000000D+00 5.585758033944684620D+00 1.124629160182848975D-01 8.875370839817151580D-01
-1.000000000000000048D-04 5.000000000000000000D-01 1.000000000000000000D+00 5.641331674102550409D+01 1.128341555584961957D-02 9.887165844441503371D-01
-1.000000000000000021D-08 5.000000000000000000D-01 1.000000000000000000D+00 5.641895779058606422D+03 1.128379163334249004D-04 9.998871620836665697D-01
-9.999999999999999452D-21 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477570534D+09 1.128379167095512970D-10 9.999999998871620388D-01
-9.999999999999999293D-41 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477568717D+19 1.128379167095512972D-20 1.000000000000000000D+00
-1.00000000000000002D-100 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477541988D+49 1.128379167095513082D-50 1.000000000000000000D+00
-9.99999999999999982D-201 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477511468D+99 1.12837916709551300D-100 1.000000000000000000D+00
-];
-// For the inversion of Shape, require only 8 digits, as
-// a consequence of bug #7569: http://bugzilla.scilab.org/show_bug.cgi?id=7569
-//
-// Some tests do not pass:
-// http://bugzilla.scilab.org/show_bug.cgi?id=8031
-// http://bugzilla.scilab.org/show_bug.cgi?id=8030
-//
-// Prints the number of accurate digits.
-precision = 1.e-12;
-precinverse = 1.e-8;
-ntests = size(table,"r");
-for i = 1 : ntests
-    x = table(i,1);
-    shape = table(i,2);
-    scale = table(i,3);
-    p = table(i,5);
-    q = table(i,6);
-    // Caution: this is the rate !
-    rate = 1/scale;
-    [p1,q1] = cdfgam("PQ",x,shape,rate);
-    x1      = cdfgam("X",shape,rate,p,q);
-    shape1  = cdfgam("Shape",rate,p,q,x);
-    rate1   = cdfgam("Rate",p,q,x,shape);
-    if ( %t ) then
-        assert_checkalmostequal ( p1 , p , precision );
-        assert_checkalmostequal ( q1 , q , precision );
-        assert_checkalmostequal ( x1 , x , precision );
-        assert_checkalmostequal ( shape1 , shape , precinverse );
-        assert_checkalmostequal ( rate1 , rate , precinverse );
-    end
-    if ( %f ) then
-        dp = assert_computedigits ( p1 , p );
-        dq = assert_computedigits ( q1 , q );
-        dx = assert_computedigits ( x1 , x );
-        ds = assert_computedigits ( shape1 , shape );
-        dr = assert_computedigits ( rate1 , rate );
-        mprintf("Test #%3d/%3d: Digits p1= %.1f, q1=%.1f, X= %.1f, S= %.1f, R= %.1f\n",i,ntests,dp,dq,dx,ds,dr);
-    end
-end
-// IEEE support
-// See http://bugzilla.scilab.org/show_bug.cgi?id=7296
-Shape = 0;
-Rate = 1;
-X = %inf; // Inf
-[P, Q] = cdfgam("PQ", X, Shape, Rate);
-assert_checkequal(P, 1);
-assert_checkequal(Q, 0);
-X = %nan; // NaN
-[P, Q] = cdfgam("PQ", X, Shape, Rate);
-assert_checkequal(P, %nan);
-assert_checkequal(Q, %nan);
index 4d3896d..0fa239b 100644 (file)
@@ -7,6 +7,7 @@
 // =============================================================================
 
 // <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
 // <-- ENGLISH IMPOSED -->
 
 // Run with test_run('statistics','cdfgam',['no_check_error_output']);
@@ -82,17 +83,13 @@ table = [
 9.99999999999999982D-201 5.000000000000000000D-01 1.000000000000000000D+00 5.641895835477511468D+99 1.12837916709551300D-100 1.000000000000000000D+00
 ];
 
-// For the inversion of Shape, require only 8 digits, as
-// a consequence of bug #7569: http://bugzilla.scilab.org/show_bug.cgi?id=7569
-//
 // Some tests do not pass:
-// http://bugzilla.scilab.org/show_bug.cgi?id=8031
-// http://bugzilla.scilab.org/show_bug.cgi?id=8030
+// http://bugzilla.scilab.org/8030
 //
 // Prints the number of accurate digits.
 
 precision = 1.e-12;
-precinverse = 1.e-8;
+precinverse = 1.e-14;
 
 ntests = size(table,"r");
 for i = 1 : ntests
@@ -125,7 +122,7 @@ for i = 1 : ntests
 end
 
 // IEEE support
-// See http://bugzilla.scilab.org/show_bug.cgi?id=7296
+// See http://bugzilla.scilab.org/7296
 Shape = 0;
 Rate = 1;