* Bug #14470 fixed - geomean often overflowed & did not check argins 68/18768/8
Samuel GOUGEON [Wed, 7 Dec 2016 20:41:34 +0000 (21:41 +0100)]
  http://bugzilla.scilab.org/14470

Change-Id: Ied3e2da4fa4e540ca886260b8834cb0ba5d59f60

scilab/CHANGES.md
scilab/modules/statistics/help/en_US/central_tendency/geomean.xml
scilab/modules/statistics/help/fr_FR/central_tendency/geomean.xml [new file with mode: 0644]
scilab/modules/statistics/help/ja_JP/central_tendency/geomean.xml [deleted file]
scilab/modules/statistics/macros/geomean.sci
scilab/modules/statistics/tests/unit_tests/geomean.tst [new file with mode: 0644]

index 1ca0cd8..361d3eb 100644 (file)
@@ -209,7 +209,7 @@ Help pages:
 * fixed / improved:  `members`, `part`, `ode`, `ode_optional_output`, `ode_roots`, `roots`,
   `printf`, `sprintf`, `iconvert`, `stdev`, `xlabel`, `and_op`, `or_op`, `%helps`
 * rewritten: `consolebox`, `double`, `isoview`, `householder`, `or`, `and`, `format`, `typeof`,
-`brackets`, `setlanguage`, `sleep`, `isinf`, `bitor`, `bitxor`, `bitand`
+`brackets`, `setlanguage`, `sleep`, `isinf`, `bitor`, `bitxor`, `bitand`, `geomean`
 * reorganized:
   - `else`, `elseif`, `end`, `try`, `sciargs`, `global`, `halt`, `empty`, `power`
   - CACSD and Signal Processing help pages have been sorted up.
@@ -365,6 +365,7 @@ Bug Fixes
 * [#14437](http://bugzilla.scilab.org/show_bug.cgi?id=14437): Problem with the affectation cmde "=" applied to a "list of struct"
 * [#14448](http://bugzilla.scilab.org/show_bug.cgi?id=14448): removed havewindow() was still documented
 * [#14461](http://bugzilla.scilab.org/show_bug.cgi?id=14461): Calling `grand(n, "markov", P, x0)` did not return all outputs.
+* [#14470](http://bugzilla.scilab.org/show_bug.cgi?id=14470): `geomean` often overflowed for easily computable entries, and did not check input arguments.
 * [#14513](http://bugzilla.scilab.org/show_bug.cgi?id=14513): `isqual` comparing two built-in functions yielded an error.
 * [#14527](http://bugzilla.scilab.org/show_bug.cgi?id=14527): Calling pathconvert function without parameters crashed Scilab.
 * [#14553](http://bugzilla.scilab.org/show_bug.cgi?id=14553): find(a=b) crashed Scilab.
index 35b19dd..9e07bf2 100644 (file)
@@ -2,6 +2,7 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2000 - INRIA - Carlos Klimann
+ * Copyright (C) 2016 - Samuel GOUGEON
  *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
  * 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:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="geomean">
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink"
+    xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML"
+    xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"
+    xml:lang="en" xml:id="geomean">
     <refnamediv>
         <refname>geomean</refname>
         <refpurpose>geometric mean</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
-        <synopsis>gm=geomean(x)
-            gm=geomean(x,'r')(or, equivalently, gm=geomean(x,1))
-            gm=geomean(x,'c')(or, equivalently, gm=geomean(x,2))
+        <synopsis>
+            gm = geomean(X)
+            GM = geomean(X, orien)   // orien: 'r'|1|'c'|2..ndims(X)
         </synopsis>
     </refsynopsisdiv>
     <refsection>
         <title>Arguments</title>
         <variablelist>
             <varlistentry>
-                <term>x</term>
+                <term>X</term>
                 <listitem>
-                    <para>real or complex vector or matrix</para>
+                    <para>
+                        Vector, matrix or hypermatrix of real or complex numbers.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>orien</term>
+                <listitem>
+                    <para>
+                        Dimension accross which the geometric average is computed. The value must be among <literal>'r', 1, 'c', 2, .. ndims(X)</literal>. Values <literal>'r'</literal> (rows) and <literal>1</literal> are equivalent, as <literal>'c'</literal> (columns) and <literal>2</literal> are.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>gm</term>
+                <listitem>
+                    <para>
+                        Scalar number: the geometric mean <literal>gm = prod(X)^(1/N)</literal>, where <literal>N = length(X)</literal> is the number of components in <literal>X</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>GM</term>
+                <listitem>
+                    <para>
+                        Vector, matrix or hypermatrix of numbers. <literal>s = size(GM)</literal> is equal to <literal>size(X)</literal>, except that <literal>s(orien)</literal> is set to 1 (due to the projected application of geomean() over components along the orien dimension).
+                    </para>
+                    <para>
+                        If <varname>X</varname> is a matrix, we have:
+                        <itemizedlist>
+                            <listitem>
+                                <literal>GM = geomean(X,1)  =>  GM(1,j) = geomean(X(:,j))</literal>
+                            </listitem>
+                            <listitem>
+                                <literal>GM = geomean(X,2)  =>  GM(i,1) = geomean(X(i,:))</literal>
+                            </listitem>
+                        </itemizedlist>
+                    </para>
                 </listitem>
             </varlistentry>
         </variablelist>
     <refsection>
         <title>Description</title>
         <para>
-            This  function computes the  geometric mean of a vector or
-            matrix <literal> x</literal>.  For a  vector or matrix  <literal> x</literal>,
-            <literal>gm=geomean(x) </literal>   returns  in  scalar  <literal>   gm</literal> the
-            geometric mean of all the entries of <literal> x</literal>.
+            <function>geomean(X,..)</function> computes the geometric mean of values stored in <literal>X</literal>.
         </para>
         <para>
-            <literal> gm=geomean(x,'r')    </literal>   (or,   equivalently,
-            <literal>gm=gmean(x,1)  </literal>  )  returns in  each   entry of the row
-            vector <literal> gm</literal> the geometric mean of each column of <literal> x</literal>.
-        </para>
-        <para>
-            <literal> gm=geomean(x,'c')    </literal>  (or,    equivalently,
-            <literal>gm=gmean(x,2) </literal>  ) returns in each  entry  of the column
-            vector <literal> gm</literal> the geometric mean of  each row of <literal> x </literal>.
+            If <varname>X</varname> stores only positive or null values, <varname>gm</varname> or <varname>GM</varname> are real. Otherwise they are most often complex.
         </para>
+        <note>
+            If <varname>X</varname> is sparse-encoded, then
+            <itemizedlist>
+                <listitem>it is reencoded in full format before being processed.</listitem>
+                <listitem>
+                    <varname>gm</varname> is always full-encoded.
+                </listitem>
+                <listitem>
+                    <varname>GM</varname> is sparse-encoded as well.
+                </listitem>
+            </itemizedlist>
+        </note>
     </refsection>
     <refsection>
         <title>Examples</title>
         <programlisting role="example"><![CDATA[
 geomean(1:10) // Returns factorial(10)^(1/10) = 4.5287286881167648
+
+// Projected geomean:
+// -----------------
+m = grand(4,5, "uin", 1, 100);
+m(3,2) = 0; m(2,4) = %inf; m(4,5) =  %nan
+geomean(m, "r")
+geomean(m, 2)
+h = grand(3,5,2, "uin",1,100)
+geomean(h,3)
+    ]]></programlisting>
+        <screen><![CDATA[
+--> m = grand(4,5, "uin", 1, 100);
+--> m(3,2) = 0; m(2,4) = %inf; m(4,5) =  %nan
+ m  =
+   13.   5.    99.   41.   20.
+   3.    92.   4.    Inf   5.
+   35.   0.    36.   40.   98.
+   86.   86.   66.   21.   Nan
+
+--> geomean(m, "r")
+ ans  =
+   18.510058   0.   31.14479   Inf   Nan
+
+--> geomean(m, 2)
+ ans  =
+   22.104082
+   Inf
+   0.
+   Nan
+
+--> h = grand(3,5,2, "uin",1,100)
+ h  =
+(:,:,1)
+   10.   40.   37.   72.   30.
+   10.   47.   54.   13.   19.
+   44.   27.   61.   10.   27.
+(:,:,2)
+   96.   88.   7.    98.   35.
+   54.   29.   96.   77.   8.
+   94.   45.   21.   46.   3.
+
+--> geomean(h,3)
+ ans  =
+   16.522712   43.150898   23.2379     36.91883    72.
+   14.142136   13.747727   64.311741   34.85685    35.79106
+   12.247449   30.983867   59.329588   16.093477   84.
+]]></screen>
+        <para></para>
+        <programlisting role="example"><![CDATA[
+// APPLICATION: Average growing rate
+// ---------------------------------
+// During 8 years, we measure the diameter D(i=1:8) of the trunc of a tree.
+D = [10 14 18 26 33 42 51 70];          // in mm
+
+// The growing rate gr(i) for year #i+1 wrt year #i is, in %:
+gr = (D(2:$)./D(1:$-1) - 1)*100
+
+// The average yearly growing rate is then, in %:
+mgr = (geomean(1+gr/100)-1)*100
+
+// If this tree had a constant growing rate, its diameter would have been:
+D(1)*(1+mgr/100)^(0:7)
  ]]></programlisting>
+        <screen><![CDATA[
+--> gr = (D(2:$)./D(1:$-1) - 1)*100
+ gr  =
+   40.   28.57   44.44   26.92   27.27   21.43   37.25
+
+--> mgr = (geomean(1+gr/100)-1)*100
+ mgr  =
+   32.05
+
+--> D(1)*(1+mgr/100)^(0:7)
+ ans  =
+   10.   13.2   17.44   23.02   30.4   40.15   53.01   70.
+]]></screen>
+    </refsection>
+    <refsection role="see also">
+        <title>See also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="prod">prod</link>
+            </member>
+            <member>
+                <link linkend="harmean">harmean</link>
+            </member>
+        </simplelist>
     </refsection>
     <refsection>
         <title>Bibliography</title>
diff --git a/scilab/modules/statistics/help/fr_FR/central_tendency/geomean.xml b/scilab/modules/statistics/help/fr_FR/central_tendency/geomean.xml
new file mode 100644 (file)
index 0000000..20259b5
--- /dev/null
@@ -0,0 +1,199 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2000 - INRIA - Carlos Klimann
+ * Copyright (C) 2016 - Samuel GOUGEON
+ *
+ * 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:mml="http://www.w3.org/1998/Math/MathML"
+    xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"
+    xml:lang="fr" xml:id="geomean">
+    <refnamediv>
+        <refname>geomean</refname>
+        <refpurpose>moyenne géométrique</refpurpose>
+    </refnamediv>
+    <refsynopsisdiv>
+        <title>Syntaxe</title>
+        <synopsis>
+            gm = geomean(X)
+            GM = geomean(X, orien)   // orien : 'r'|1|'c'|2..ndims(X)
+        </synopsis>
+    </refsynopsisdiv>
+    <refsection>
+        <title>Paramètres</title>
+        <variablelist>
+            <varlistentry>
+                <term>X</term>
+                <listitem>
+                    <para>
+                        Vecteur, matrice ou hypermatrice de réels ou de complexes.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>orien</term>
+                <listitem>
+                    <para>
+                        Dimension selon laquelle la moyenne géométrique est calculée. Sa valeur doit être parmi <literal>'r', 1, 'c', 2, .. ndims(X)</literal>. Les valeurs <literal>'r'</literal> (lignes) et <literal>1</literal> sont équivalentes, comme <literal>'c'</literal> (colonnes) et <literal>2</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>gm</term>
+                <listitem>
+                    <para>
+                        Réel scalaire : la moyenne géométrique <literal>gm = prod(X)^(1/N)</literal>, où <literal>N = length(X)</literal> est le nombre d'éléments de <literal>X</literal>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>GM</term>
+                <listitem>
+                    <para>
+                        Vecteur, matrice ou hypermatrice de nombres. <literal>s = size(GM)</literal> est égal à <literal>size(X)</literal>, sauf que <literal>s(orien)</literal> vaut 1 (dû à l'application projetée de geomean() selon la dimension orien).
+                    </para>
+                    <para>
+                        Si <varname>X</varname> est une matrice, on a :
+                        <itemizedlist>
+                            <listitem>
+                                <literal>GM = geomean(X,1)  =>  GM(1,j) = geomean(X(:,j))</literal>
+                            </listitem>
+                            <listitem>
+                                <literal>GM = geomean(X,2)  =>  GM(i,1) = geomean(X(i,:))</literal>
+                            </listitem>
+                        </itemizedlist>
+                    </para>
+                </listitem>
+            </varlistentry>
+        </variablelist>
+    </refsection>
+    <refsection>
+        <title>Description</title>
+        <para>
+            <function>geomean(X,..)</function> calcule la moyenne géométrique des valeurs de <literal>X</literal>.
+        </para>
+        <para>
+            Si <varname>X</varname> ne contient que des valeurs positives ou nulles, <varname>gm</varname> ou <varname>GM</varname> sont réels. Sinon ils sont souvent complexes.
+        </para>
+        <note>
+            Si <varname>X</varname> est une matrice creuse (sparse), alors
+            <itemizedlist>
+                <listitem>elle est convertie en matrice pleine avant d'être traitée.</listitem>
+                <listitem>
+                    <varname>gm</varname> est toujours pleine.
+                </listitem>
+                <listitem>
+                    <varname>GM</varname> est creuse.
+                </listitem>
+            </itemizedlist>
+        </note>
+    </refsection>
+    <refsection>
+        <title>Exemples</title>
+        <programlisting role="example"><![CDATA[
+geomean(1:10) // Retourne factorial(10)^(1/10) = 4.5287286881167648
+
+// Moyenne géométrique selon une direction choisie :
+// -----------------------------------------------
+m = grand(4,5, "uin", 1, 100);
+m(3,2) = 0; m(2,4) = %inf; m(4,5) =  %nan
+geomean(m, "r")
+geomean(m, 2)
+h = grand(3,5,2, "uin",1,100)
+geomean(h,3)
+    ]]></programlisting>
+        <screen><![CDATA[
+--> m = grand(4,5, "uin", 1, 100);
+--> m(3,2) = 0; m(2,4) = %inf; m(4,5) =  %nan
+ m  =
+   13.   5.    99.   41.   20.
+   3.    92.   4.    Inf   5.
+   35.   0.    36.   40.   98.
+   86.   86.   66.   21.   Nan
+
+--> geomean(m, "r")
+ ans  =
+   18.510058   0.   31.14479   Inf   Nan
+
+--> geomean(m, 2)
+ ans  =
+   22.104082
+   Inf
+   0.
+   Nan
+
+--> h = grand(3,5,2, "uin",1,100)
+ h  =
+(:,:,1)
+   10.   40.   37.   72.   30.
+   10.   47.   54.   13.   19.
+   44.   27.   61.   10.   27.
+(:,:,2)
+   96.   88.   7.    98.   35.
+   54.   29.   96.   77.   8.
+   94.   45.   21.   46.   3.
+
+--> geomean(h,3)
+ ans  =
+   30.983867   59.329588   16.093477   84.         32.403703
+   23.2379     36.91883    72.         31.638584   12.328828
+   64.311741   34.85685    35.79106    21.447611   9.
+]]></screen>
+        <para></para>
+        <programlisting role="example"><![CDATA[
+// APPLICATION : Taux de croissance moyen
+// --------------------------------------
+// Pendant 8 ans, on mesure le diamètre D(i=1:8) du tronc d'un arbre.
+D = [10 14 18 26 33 42 51 70];          // en mm
+
+// Le taux de croissance gr(i) pour l'année #i+1 par rapport à l'année #i est, en % :
+gr = (D(2:$)./D(1:$-1) - 1)*100
+
+// Le taux de croissance moyen est donc, en % :
+mgr = (geomean(1+gr/100)-1)*100
+
+// Si ce tronc avait un taux de croissance constant, son diamètre aurait été :
+D(1)*(1+mgr/100)^(0:7)
+ ]]></programlisting>
+        <screen><![CDATA[
+--> gr = (D(2:$)./D(1:$-1) - 1)*100
+ gr  =
+   40.   28.57   44.44   26.92   27.27   21.43   37.25
+
+--> mgr = (geomean(1+gr/100)-1)*100
+ mgr  =
+   32.05
+
+--> D(1)*(1+mgr/100)^(0:7)
+ ans  =
+   10.   13.2   17.44   23.02   30.4   40.15   53.01   70.
+]]></screen>
+    </refsection>
+    <refsection role="see also">
+        <title>Voir aussi</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="prod">prod</link>
+            </member>
+            <member>
+                <link linkend="harmean">harmean</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection>
+        <title>Bibliographie</title>
+        <para>
+            Wonacott, T.H. &amp; Wonacott, R.J.; Introductory Statistics, fifth edition, J.Wiley &amp; Sons, 1990.
+        </para>
+    </refsection>
+</refentry>
diff --git a/scilab/modules/statistics/help/ja_JP/central_tendency/geomean.xml b/scilab/modules/statistics/help/ja_JP/central_tendency/geomean.xml
deleted file mode 100644 (file)
index 684ffe0..0000000
+++ /dev/null
@@ -1,68 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2000 - INRIA - Carlos Klimann
- *
- * 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:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="ja" xml:id="geomean">
-    <refnamediv>
-        <refname>geomean</refname>
-        <refpurpose>幾何学的平均</refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>呼び出し手順</title>
-        <synopsis>gm=geomean(x)
-            gm=geomean(x,'r')(or, equivalently, gm=geomean(x,1))
-            gm=geomean(x,'c')(or, equivalently, gm=geomean(x,2))
-        </synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>パラメータ</title>
-        <variablelist>
-            <varlistentry>
-                <term>x</term>
-                <listitem>
-                    <para>実数または複素数のベクトルまたは行列</para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>説明</title>
-        <para>
-            この関数は,ベクトルまたは行列<literal>x</literal>の幾何学的平均
-            を計算します.
-            ベクトルまたは行列<literal>x</literal>に関して,
-            <literal>gm=geomean(x)</literal>は,
-            <literal>x</literal>の全てのエントリの幾何学的平均を
-            スカラー<literal>gm</literal>に返します.
-        </para>
-        <para>
-            <literal> gm=geomean(x,'r')    </literal>   (または等価な
-            <literal>gm=gmean(x,1)  </literal>  )  は,
-            行ベクトル<literal>gm</literal>の各エントリに,<literal>x</literal>の
-            各列の幾何学的平均を返します.
-        </para>
-        <para>
-            <literal> gm=geomean(x,'c')    </literal>  (または等価な
-            <literal>gm=gmean(x,2) </literal>  ) は,
-            列ベクトル<literal>gm</literal>の各エントリに,<literal>x</literal>の
-            各行の幾何学的平均を返します.
-        </para>
-    </refsection>
-    <refsection>
-        <title>参考文献</title>
-        <para>
-            Wonacott, T.H. &amp; Wonacott, R.J.; Introductory Statistics, fifth edition, J.Wiley &amp; Sons, 1990.
-        </para>
-    </refsection>
-</refentry>
index 8e2c34b..dc85216 100644 (file)
@@ -1,7 +1,7 @@
-
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 1999 - INRIA - Carlos Klimann
-// Copyright (C) 2010 - Samuel GOUGEON : bug 7204, code style
+// Copyright (C) 2010, 2016 - Samuel GOUGEON
+// Copyright (C) 2016 - Michael Baudin
 //
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
 //
 // along with this program.
 //
 
-function gm = geomean(x,orien)
-    //
-    //This function computes the geometric mean of a vector or matrix x.
-    //
-    //For a vector  or  matrix  x,  gm=geomean(x) returns in  scalar  gm the
-    //geometric mean of all the entries of x.
-    //
-    //gm=geomean(x,'r')(or,  equivalently,   gm=gmean(x,1)) returns  in each
-    //entry of the row vector gm the geometric mean of each column of x.
-    //
-    //gm=geomean(x,'c')(or,   equivalently,  gm=gmean(x,2)) returns in  each
-    //entry of the column vector gm the geometric mean of each row of x.
-    //
-    //References:  Wonacott, T.H. & Wonacott, R.J.; Introductory
-    //Statistics, J.Wiley & Sons, 1990.
-    //
-    //
-    [lhs, rhs] = argn(0)
-    if rhs == 0 then
-        msg = gettext("%s: Wrong number of input arguments: %d to %d expected.\n")
-        error(msprintf(msg, "geomean",1,2))
+function gm = geomean(x, orien)
+
+    rhs = argn(2)
+    if rhs==0 | rhs>2 then
+        msg = _("%s: Wrong number of input arguments: %d to %d expected.\n")
+        error(msprintf(msg, "geomean", 1, 2))
+    end
+    sp = typeof(x)=="sparse"
+    if sp then
+        x = full(x)
+    end
+    if typeof(x)~="constant" then
+        msg = _("%s: Argument #%d: Decimal or complex number(s) expected.\n")
+        error(msprintf(msg, "geomean", 1))
     end
-    if x == [] then
+    if isdef("orien", "l") then
+        if type(orien)==1
+            orien = int(orien)
+        end
+        if ~((type(orien)==10 && size(orien,"*")==1 && or(convstr(orien)==["r" "c"])) | ..
+            (type(orien)==1  && size(orien,"*")==1 && orien>0 && orien<=ndims(x)) )
+            msg = _("%s: Argument #%d: Must be in the set {%s}.\n")
+            tmp = msprintf("""r"",""c"",1..%d", ndims(x))
+            error(msprintf(msg, "geomean", 2, tmp))
+        end
+    end
+    if (x == []) then
         gm = %nan
         return
     end
+
+    mod = ieee()
+    ieee(2)     // To avoid procedural warning or error on log(0)
     if rhs == 1 then
-        gm = prod(x)^(1/size(x,"*"))
+        gm = exp(mean(log(x)))
     elseif rhs==2
-        gm = prod(x,orien).^(1/size(x,orien))
-    else
-        msg = gettext("%s: Wrong number of input arguments: %d to %d expected.\n")
-        error(msprintf(msg, "geomean",1,2))
+        gm = exp(mean(log(x),orien))
+    end
+    ieee(mod)   // Restoring initial mode
+
+    if sp & rhs==2 then
+        gm = sparse(gm)
     end
 endfunction
diff --git a/scilab/modules/statistics/tests/unit_tests/geomean.tst b/scilab/modules/statistics/tests/unit_tests/geomean.tst
new file mode 100644 (file)
index 0000000..b475cdd
--- /dev/null
@@ -0,0 +1,82 @@
+// =======================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2016 - Samuel GOUGEON
+//
+//  This file is distributed under the same license as the Scilab package.
+// =======================================================================
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+//
+// No specific expected console output to be checked against:
+// <-- NO CHECK REF -->
+// <-- CLI SHELL MODE -->
+//
+
+x = grand(50,40,"uin", 1, 1000);
+y = geomean(x);
+yi = geomean(1 ./ x);
+assert_checkequal(y*yi, 1);
+yc = geomean([x, 1../x], "c");
+assert_checkalmostequal(yc, ones(x(:,1)), 10*%eps);
+yr = geomean([x ; 1../x], "r");
+assert_checkalmostequal(yr, ones(x(1,:)), 10*%eps);
+
+// With a special spreading value somewhere
+sv = [0 %inf -%inf %nan];
+for v = sv
+    xs = x;
+    xs(20,30) = v;
+    if v==-%inf
+        r = %inf*(1+%i);
+    else
+        r = v;
+    end
+    assert_checkequal(geomean(xs), r);
+    ysc = geomean(xs, "c");
+    ysc(20) = r;
+    assert_checkalmostequal(geomean(xs, "c"), ysc, 10*%eps);
+    ysr = geomean(xs, "r");
+    ysr(30) = r;
+    assert_checkalmostequal(geomean(xs, "r"), ysr, 10*%eps);
+end
+
+// With an homogeneous matrix, possibly negative or complex or huge value
+sv = [1 2 -2 %i -%i 1.e234];
+for v = sv
+    xs = ones(50, 40)*v;
+    assert_checkalmostequal(clean(geomean(xs)), v, 1e-11);
+    assert_checkalmostequal(clean(geomean(xs, "c")), ones(x(:,1))*v, 1e-11);
+    assert_checkalmostequal(clean(geomean(xs, "r")), ones(x(1,:))*v, 1e-11);
+end
+
+// Hypermatrices
+x = grand(30,20,10, "uin", 1, 1000);
+y = geomean(x);
+yi = geomean(1 ./ x);
+assert_checkequal(y*yi, 1);
+yc = geomean([x, 1../x], "c");
+assert_checkalmostequal(yc, ones(x(:,1,:)), 10*%eps);
+yr = geomean([x ; 1../x], "r");
+assert_checkalmostequal(yr, ones(x(1,:,:)), 10*%eps);
+xe = x;
+xe(:,:,11:20) = 1../x;
+y3 = geomean(xe, 3);
+assert_checkalmostequal(y3, ones(x(:,:,1)), 10*%eps);
+
+// Forbidden data types
+assert_checkfalse(execstr("geomean(x>500)", "errcatch")==0);
+assert_checkfalse(execstr("geomean(int16(x))", "errcatch")==0);
+assert_checkfalse(execstr("geomean(uint16(x))", "errcatch")==0);
+assert_checkfalse(execstr("geomean(x+0*%z)", "errcatch")==0);
+assert_checkfalse(execstr("geomean(sparse(x))", "errcatch")==0);
+// Forbidden orientations
+assert_checkfalse(execstr("geomean(x,''g'')", "errcatch")==0);
+assert_checkfalse(execstr("geomean(x,[''c'' ''r''])", "errcatch")==0);
+assert_checkfalse(execstr("geomean(x,[1 2])", "errcatch")==0);
+assert_checkfalse(execstr("geomean(x,-2)", "errcatch")==0);
+assert_checkfalse(execstr("geomean(x,0)", "errcatch")==0);
+assert_checkfalse(execstr("geomean(x,4)", "errcatch")==0);