* 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.
* [#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.
<!--
* 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>
--- /dev/null
+<?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. & Wonacott, R.J.; Introductory Statistics, fifth edition, J.Wiley & Sons, 1990.
+ </para>
+ </refsection>
+</refentry>
+++ /dev/null
-<?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. & Wonacott, R.J.; Introductory Statistics, fifth edition, J.Wiley & Sons, 1990.
- </para>
- </refsection>
-</refentry>
-
// 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
--- /dev/null
+// =======================================================================
+// 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);