CMSCOPE changed

32 files changed, 4233 insertions, 0 deletions

diff --git a/modules/randlib/tests/nonreg_tests/bug_12606.dia.ref b/modules/randlib/tests/nonreg_tests/bug_12606.dia.ref
new file mode 100755
index 000000000..d02d4404f
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_12606.dia.ref
@@ -0,0 +1,25 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 12606 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12606
+//
+// <-- Short Description -->
+// overloads names for grand() are not standard: %#_Rand() instead of %#_grand()
+// Save warning mode and turn off warnings
+status = warning("query");
+warning("off");
+function r=%b_grand(arg)
+    r="boolean overload for grand is %b_grand() as expected";
+endfunction
+assert_checkequal(grand(%t),"boolean overload for grand is %b_grand() as expected");
+// Revert to the previous warning mode
+warning(status);
diff --git a/modules/randlib/tests/nonreg_tests/bug_12606.tst b/modules/randlib/tests/nonreg_tests/bug_12606.tst
new file mode 100755
index 000000000..849c8d293
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_12606.tst
@@ -0,0 +1,28 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 12606 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12606
+//
+// <-- Short Description -->
+// overloads names for grand() are not standard: %#_Rand() instead of %#_grand()
+
+// Save warning mode and turn off warnings
+status = warning("query");
+warning("off");
+
+function r=%b_grand(arg)
+    r="boolean overload for grand is %b_grand() as expected";
+endfunction
+assert_checkequal(grand(%t),"boolean overload for grand is %b_grand() as expected");
+
+// Revert to the previous warning mode
+warning(status);
diff --git a/modules/randlib/tests/nonreg_tests/bug_1568.dia.ref b/modules/randlib/tests/nonreg_tests/bug_1568.dia.ref
new file mode 100755
index 000000000..6fda586f5
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_1568.dia.ref
@@ -0,0 +1,36 @@
+// <-- Non-regression test for bug 1568 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=1568
+//
+// <-- Short Description -->
+//    The random number generator grand() does not always produce the same
+//    result starting from the same seed, as shown by the following
+//    transcript of a Scilab session:
+//
+//    $ scilab -nw
+//            -------------------------------------------
+//                             scilab-3.1.1
+//
+//                      Copyright (c) 1989-2005
+//                  Consortium Scilab (INRIA, ENPC)
+//            -------------------------------------------
+//
+//
+//    Startup execution:
+//      loading initial environment
+//
+//    -->grand('setsd',12)
+//
+//    -->grand(1,'prm',[1:5]')
+//     ans  =
+// ...
+// Copyright INRIA
+// Scilab Project - Pierre MARECHAL
+// Copyright INRIA 2005
+// Date : 21 octobre 2005
+grand('setsd',12);
+A = grand(1,'prm',[1:5]');
+grand('setsd',12);
+B = grand(1,'prm',[1:5]');
+if or(A<>B) then bugmes();quit;end
diff --git a/modules/randlib/tests/nonreg_tests/bug_1568.tst b/modules/randlib/tests/nonreg_tests/bug_1568.tst
new file mode 100755
index 000000000..2cdbe0351
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_1568.tst
@@ -0,0 +1,41 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2005-2008 - INRIA - Pierre MARECHAL <pierre.marechal@inria.fr>
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- Non-regression test for bug 1568 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=1568
+//
+// <-- Short Description -->
+//    The random number generator grand() does not always produce the same 
+//    result starting from the same seed, as shown by the following 
+//    transcript of a Scilab session:
+//
+//    $ scilab -nw
+//            -------------------------------------------
+//                             scilab-3.1.1
+//
+//                      Copyright (c) 1989-2005
+//                  Consortium Scilab (INRIA, ENPC)
+//            -------------------------------------------
+//
+//
+//    Startup execution:
+//      loading initial environment
+//
+//    -->grand('setsd',12)
+//
+//    -->grand(1,'prm',[1:5]')
+//     ans  =
+// ...
+
+grand('setsd',12);
+A = grand(1,'prm',[1:5]');
+grand('setsd',12);
+B = grand(1,'prm',[1:5]');
+
+if or(A<>B) then pause,end
diff --git a/modules/randlib/tests/nonreg_tests/bug_2964.dia.ref b/modules/randlib/tests/nonreg_tests/bug_2964.dia.ref
new file mode 100755
index 000000000..5ad73091a
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_2964.dia.ref
@@ -0,0 +1,18 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2005-2008 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- Non-regression test for bug 2964 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=2964
+//
+// <-- Short Description -->
+//    the function grand with option 'mul' does not accept probability one as an argument.
+y=grand(1,'mul',1,1);
+if or(y<>[1;0]) then bugmes();quit;end
+y=grand(1,'mul',1,ones(10,1)/10);
+if y($)<>0 then bugmes();quit;end
+if sum(y)<>1 then bugmes();quit;end
diff --git a/modules/randlib/tests/nonreg_tests/bug_2964.tst b/modules/randlib/tests/nonreg_tests/bug_2964.tst
new file mode 100755
index 000000000..1c8a0327b
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_2964.tst
@@ -0,0 +1,20 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2005-2008 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- Non-regression test for bug 2964 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=2964
+//
+// <-- Short Description -->
+//    the function grand with option 'mul' does not accept probability one as an argument. 
+y=grand(1,'mul',1,1);
+if or(y<>[1;0]) then pause,end
+
+y=grand(1,'mul',1,ones(10,1)/10);
+if y($)<>0 then pause,end
+if sum(y)<>1 then pause,end
diff --git a/modules/randlib/tests/nonreg_tests/bug_4839.dia.ref b/modules/randlib/tests/nonreg_tests/bug_4839.dia.ref
new file mode 100755
index 000000000..9d59ba4b2
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_4839.dia.ref
@@ -0,0 +1,16 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- Non-regression test for bug 4839 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=4839
+//
+// <-- Short Description -->
+//    grand does not perform well with multivariate gaussian
+A = [5 3 1;3 5 1; 1 1 5];
+computed = grand(10,'mn',[0  0 0]',A);
+assert_checkalmostequal ( size(computed) , [3 10] );
diff --git a/modules/randlib/tests/nonreg_tests/bug_4839.tst b/modules/randlib/tests/nonreg_tests/bug_4839.tst
new file mode 100755
index 000000000..fd4722948
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_4839.tst
@@ -0,0 +1,20 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- Non-regression test for bug 4839 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=4839
+//
+// <-- Short Description -->
+//    grand does not perform well with multivariate gaussian
+
+
+
+A = [5 3 1;3 5 1; 1 1 5];
+computed = grand(10,'mn',[0  0 0]',A);
+assert_checkalmostequal ( size(computed) , [3 10] );
+
diff --git a/modules/randlib/tests/nonreg_tests/bug_5207.dia.ref b/modules/randlib/tests/nonreg_tests/bug_5207.dia.ref
new file mode 100755
index 000000000..fa6f1423a
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_5207.dia.ref
@@ -0,0 +1,45 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Antoine ELIAS
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 5207 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=5207
+//
+// <-- Short Description -->
+// grand did not managed 3-D inputs/outputs
+//dimension as input parameters
+//init mt generator with seed = 1337
+grand("setsd", 1337);
+a = grand(3, 4, 5, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5, "def");
+ref = matrix(ref, [3, 4, 5]);
+assert_checkequal(a, ref);
+grand("setsd", 1337);
+a = grand(3, 4, 5, 6, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5 * 6, "def");
+ref = matrix(ref, [3, 4, 5, 6]);
+assert_checkequal(a, ref);
+//dimension from size of input parameter
+grand("setsd", 1337);
+Z = zeros(3,4,5);
+a = grand(Z, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5, "def");
+ref = matrix(ref, [3, 4, 5]);
+assert_checkequal(a, ref);
+grand("setsd", 1337);
+Z = zeros(3,4,5,6);
+a = grand(Z, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5 * 6, "def");
+ref = matrix(ref, [3, 4, 5, 6]);
+assert_checkequal(a, ref);
diff --git a/modules/randlib/tests/nonreg_tests/bug_5207.tst b/modules/randlib/tests/nonreg_tests/bug_5207.tst
new file mode 100755
index 000000000..ea54c2fef
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_5207.tst
@@ -0,0 +1,51 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Antoine ELIAS
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 5207 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=5207
+//
+// <-- Short Description -->
+// grand did not managed 3-D inputs/outputs
+
+//dimension as input parameters
+
+//init mt generator with seed = 1337
+grand("setsd", 1337);
+a = grand(3, 4, 5, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5, "def");
+ref = matrix(ref, [3, 4, 5]);
+assert_checkequal(a, ref);
+
+grand("setsd", 1337);
+a = grand(3, 4, 5, 6, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5 * 6, "def");
+ref = matrix(ref, [3, 4, 5, 6]);
+assert_checkequal(a, ref);
+
+//dimension from size of input parameter
+grand("setsd", 1337);
+Z = zeros(3,4,5);
+a = grand(Z, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5, "def");
+ref = matrix(ref, [3, 4, 5]);
+assert_checkequal(a, ref);
+
+grand("setsd", 1337);
+Z = zeros(3,4,5,6);
+a = grand(Z, "def");
+grand("setsd", 1337);
+ref = grand(1, 3 * 4 * 5 * 6, "def");
+ref = matrix(ref, [3, 4, 5, 6]);
+assert_checkequal(a, ref);
+
diff --git a/modules/randlib/tests/nonreg_tests/bug_6690.dia.ref b/modules/randlib/tests/nonreg_tests/bug_6690.dia.ref
new file mode 100755
index 000000000..b5afacd00
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_6690.dia.ref
@@ -0,0 +1,178 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// Set the seed to always get the same random numbers
+rand("seed", 0);
+grand("setsd", 0);
+//
+// 1. Permute some column vectors, and check that the output is basically correct.
+//
+// With Complexes.
+//
+X = (2:10)';
+P = grand(5,"prm",X*%i);
+assert_checkequal ( isreal(P) , %f );
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , X*%i*ones(1,5) );
+//
+// With Integers.
+//
+P = grand(5,"prm",int8(X));
+assert_checkequal ( typeof(P) , "int8" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , int8(X)*ones(1,5) );
+//
+P = grand(5,"prm",int16(X));
+assert_checkequal ( typeof(P) , "int16" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , int16(X)*ones(1,5) );
+//
+P = grand(5,"prm",int32(X));
+assert_checkequal ( typeof(P) , "int32" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , int32(X)*ones(1,5) );
+//
+// With Strings.
+//
+X = string(X);
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "string" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , ["3" "9" "5" "7" "10" "6" "8" "2" "4"]');
+//
+// With Booleans.
+//
+X = [%f %t %t %t %t %f]';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "boolean" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , [%t %f %t %t %f %t]' );
+//
+// With Polynomials.
+//
+s = poly(0,"s");
+X = [0 s 1+s 1+s^2 s^4]';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "polynomial" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , [1+s^2 1+s s s^4 0]' );
+//
+// With Sparses.
+//
+X = sparse([1 2 3 0 0 0 0 5 6])';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "sparse" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , sparse([0 0 0 3 5 6 2 1 0])' );
+//
+// 2. With row vectors / matrices.
+//
+X = 2:10;
+P = grand(1,"prm",X);
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+X(2, :) = 12:20;
+P = grand(2,"prm",X*%i);
+refP = [
+13  9   20  3   10  2   4  5  15
+19  14  6   18  12  16  7  8  17 ];
+refP(:, :, 2) = [
+14  7  8  9   13  10  15  19  4
+12  5  2  17  6   16  18  3   20 ];
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( P, refP*%i );
+//
+// With Integers.
+//
+X_int = int32(X);
+P = grand(2,"prm",X_int);
+refP = int32([
+5  17  20  7   15  10  2   13  4
+3  16  8   12  19  9   18  14  6 ]);
+refP(:, :, 2) = int32([
+16  13  6  2  20  19  18  14  7
+9   12  3  8  15  17  10  4   5 ]);
+assert_checkequal ( size(P) , [size(X_int) 2] );
+assert_checkequal ( P, refP );
+//
+// With Strings.
+//
+X_str = string(X);
+P = grand(2,"prm",X_str);
+refP = string([
+19  7   14  10  18  17  5  3  4
+9   15  2   20  13  12  6  8  16 ]);
+refP(:, :, 2) = string([
+13  20  9   6  8   14  3  4   16
+15  5   19  2  18  12  7  10  17 ]);
+assert_checkequal ( size(P) , [size(X_str) 2] );
+assert_checkequal ( P, refP );
+//
+// With Hypermatrices.
+//
+X(:, :, 2) = X+20;
+P = grand(2,"prm",X);
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( members(P, 2:40), ones(P) );
+//
+// Of Complexes.
+//
+P = grand(2,"prm",X*%i);
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( members(P, (2:40)*%i), ones(P) );
+//
+// Of Integers.
+//
+P = grand(2,"prm",int16(X));
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( members(P, int16(2:40)), ones(P) );
+//
+// Of Strings.
+//
+Xs = string(X);
+Ps = grand(2,"prm",Xs);
+assert_checkequal ( size(Ps) , [size(Xs) 2] );
+assert_checkequal ( members(Ps, string(2:40)), ones(Ps) );
+//
+// Of Booleans.
+//
+Xb = floor(10*rand(2, 9, 2));
+Xb = Xb<5;
+Pb = grand(2,"prm",Xb);
+refPb = gsort(Pb+0,"g","i");
+assert_checkequal ( size(Pb) , [size(Xb) 2] );
+assert_checkequal ( refPb(1:36), zeros(36, 1) );
+assert_checkequal ( refPb(37:72), ones(36, 1) );
+//
+// Of Ploynomials.
+//
+Xp = matrix(X(:)+s, 2, 9, 2);
+Pp = grand(2,"prm",Xp);
+refPp = matrix(Pp(:)-s, 2, 9, 2, 2);
+refPp = gsort(coeff(refPp(:)),"g","i");
+assert_checkequal ( size(Pp) , [size(Xp) 2] );
+assert_checkequal ( members(2:40, refPp), 2*([ones(1,9) 0 ones(1,9) 0 ones(1,9) 0 ones(1,9)]) );
+//
+// Sparse hypermatrices do not exist yet.
+//
+//
+// Bug #6689
+//
+X = [%i 1-%i 2+3*%i].';
+P = grand(4, "prm", X);
+refP = [
+2+3*%i  %i      2+3*%i  1-%i
+%i      2+3*%i  %i      %i
+1-%i    1-%i    1-%i    2+3*%i ];
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , [3 4] );
+assert_checkequal ( P , refP );
+refMsg = msprintf(_("%s: Wrong type for input argument: Matrix (full or sparse) or Hypermatrix of Reals, Complexes, Integers, Booleans, Strings or Polynomials expected.\n"), "grand");
+assert_checkerror("grand(2, ""prm"", list())", refMsg);
diff --git a/modules/randlib/tests/nonreg_tests/bug_6690.tst b/modules/randlib/tests/nonreg_tests/bug_6690.tst
new file mode 100755
index 000000000..331fa85b3
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_6690.tst
@@ -0,0 +1,182 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// Set the seed to always get the same random numbers
+rand("seed", 0);
+grand("setsd", 0);
+
+//
+// 1. Permute some column vectors, and check that the output is basically correct.
+//
+// With Complexes.
+//
+X = (2:10)';
+P = grand(5,"prm",X*%i);
+assert_checkequal ( isreal(P) , %f );
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , X*%i*ones(1,5) );
+//
+// With Integers.
+//
+P = grand(5,"prm",int8(X));
+assert_checkequal ( typeof(P) , "int8" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , int8(X)*ones(1,5) );
+//
+P = grand(5,"prm",int16(X));
+assert_checkequal ( typeof(P) , "int16" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , int16(X)*ones(1,5) );
+//
+P = grand(5,"prm",int32(X));
+assert_checkequal ( typeof(P) , "int32" );
+assert_checkequal ( size(P) , [9 5] );
+assert_checkequal ( gsort(P,"r","i") , int32(X)*ones(1,5) );
+//
+// With Strings.
+//
+X = string(X);
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "string" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , ["3" "9" "5" "7" "10" "6" "8" "2" "4"]');
+//
+// With Booleans.
+//
+X = [%f %t %t %t %t %f]';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "boolean" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , [%t %f %t %t %f %t]' );
+//
+// With Polynomials.
+//
+s = poly(0,"s");
+X = [0 s 1+s 1+s^2 s^4]';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "polynomial" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , [1+s^2 1+s s s^4 0]' );
+//
+// With Sparses.
+//
+X = sparse([1 2 3 0 0 0 0 5 6])';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "sparse" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( P , sparse([0 0 0 3 5 6 2 1 0])' );
+
+//
+// 2. With row vectors / matrices.
+//
+X = 2:10;
+P = grand(1,"prm",X);
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+X(2, :) = 12:20;
+P = grand(2,"prm",X*%i);
+refP = [
+13  9   20  3   10  2   4  5  15
+19  14  6   18  12  16  7  8  17 ];
+refP(:, :, 2) = [
+14  7  8  9   13  10  15  19  4
+12  5  2  17  6   16  18  3   20 ];
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( P, refP*%i );
+//
+// With Integers.
+//
+X_int = int32(X);
+P = grand(2,"prm",X_int);
+refP = int32([
+5  17  20  7   15  10  2   13  4
+3  16  8   12  19  9   18  14  6 ]);
+refP(:, :, 2) = int32([
+16  13  6  2  20  19  18  14  7
+9   12  3  8  15  17  10  4   5 ]);
+assert_checkequal ( size(P) , [size(X_int) 2] );
+assert_checkequal ( P, refP );
+//
+// With Strings.
+//
+X_str = string(X);
+P = grand(2,"prm",X_str);
+refP = string([
+19  7   14  10  18  17  5  3  4
+9   15  2   20  13  12  6  8  16 ]);
+refP(:, :, 2) = string([
+13  20  9   6  8   14  3  4   16
+15  5   19  2  18  12  7  10  17 ]);
+assert_checkequal ( size(P) , [size(X_str) 2] );
+assert_checkequal ( P, refP );
+//
+// With Hypermatrices.
+//
+X(:, :, 2) = X+20;
+P = grand(2,"prm",X);
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( members(P, 2:40), ones(P) );
+//
+// Of Complexes.
+//
+P = grand(2,"prm",X*%i);
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( members(P, (2:40)*%i), ones(P) );
+//
+// Of Integers.
+//
+P = grand(2,"prm",int16(X));
+assert_checkequal ( size(P) , [size(X) 2] );
+assert_checkequal ( members(P, int16(2:40)), ones(P) );
+//
+// Of Strings.
+//
+Xs = string(X);
+Ps = grand(2,"prm",Xs);
+assert_checkequal ( size(Ps) , [size(Xs) 2] );
+assert_checkequal ( members(Ps, string(2:40)), ones(Ps) );
+//
+// Of Booleans.
+//
+Xb = floor(10*rand(2, 9, 2));
+Xb = Xb<5;
+Pb = grand(2,"prm",Xb);
+refPb = gsort(Pb+0,"g","i");
+assert_checkequal ( size(Pb) , [size(Xb) 2] );
+assert_checkequal ( refPb(1:36), zeros(36, 1) );
+assert_checkequal ( refPb(37:72), ones(36, 1) );
+//
+// Of Ploynomials.
+//
+Xp = matrix(X(:)+s, 2, 9, 2);
+Pp = grand(2,"prm",Xp);
+refPp = matrix(Pp(:)-s, 2, 9, 2, 2);
+refPp = gsort(coeff(refPp(:)),"g","i");
+assert_checkequal ( size(Pp) , [size(Xp) 2] );
+assert_checkequal ( members(2:40, refPp), 2*([ones(1,9) 0 ones(1,9) 0 ones(1,9) 0 ones(1,9)]) );
+//
+// Sparse hypermatrices do not exist yet.
+//
+
+//
+// Bug #6689
+//
+X = [%i 1-%i 2+3*%i].';
+P = grand(4, "prm", X);
+refP = [
+2+3*%i  %i      2+3*%i  1-%i
+%i      2+3*%i  %i      %i
+1-%i    1-%i    1-%i    2+3*%i ];
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , [3 4] );
+assert_checkequal ( P , refP );
+
+refMsg = msprintf(_("%s: Wrong type for input argument: Matrix (full or sparse) or Hypermatrix of Reals, Complexes, Integers, Booleans, Strings or Polynomials expected.\n"), "grand");
+assert_checkerror("grand(2, ""prm"", list())", refMsg);
diff --git a/modules/randlib/tests/nonreg_tests/bug_8597.dia.ref b/modules/randlib/tests/nonreg_tests/bug_8597.dia.ref
new file mode 100755
index 000000000..2e5ffa3a2
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_8597.dia.ref
@@ -0,0 +1,27 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- ENGLISH IMPOSED -->
+//
+// <-- Non-regression test for bug 8597 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8597
+//
+// <-- Short Description -->
+// Uncontrolled message of grand/clcg4 should be displayed as warning
+warning("off");
+grand('setgen',"clcg4");
+grand('setsd',123456,123456,123456,123456);
+warning("on");
+grand('setgen',"clcg4");
+grand('setsd',123456,123456,123456,123456);
+WARNING: be aware that you may have lost synchronization
+    between the virtual generator %d and the others.
+    use grand("setall", s1, s2, s3, s4) if you want to recover it.
diff --git a/modules/randlib/tests/nonreg_tests/bug_8597.tst b/modules/randlib/tests/nonreg_tests/bug_8597.tst
new file mode 100755
index 000000000..e861fc04e
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_8597.tst
@@ -0,0 +1,26 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- ENGLISH IMPOSED -->
+//
+// <-- Non-regression test for bug 8597 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8597
+//
+// <-- Short Description -->
+// Uncontrolled message of grand/clcg4 should be displayed as warning
+
+warning("off");
+grand('setgen',"clcg4");
+grand('setsd',123456,123456,123456,123456);
+
+warning("on");
+grand('setgen',"clcg4");
+grand('setsd',123456,123456,123456,123456);
diff --git a/modules/randlib/tests/nonreg_tests/bug_9159.dia.ref b/modules/randlib/tests/nonreg_tests/bug_9159.dia.ref
new file mode 100755
index 000000000..3f3f36643
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_9159.dia.ref
@@ -0,0 +1,18 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009 - Digiteo - Allan CORNET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 9159 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=9159
+//
+// <-- Short Description -->
+// grand(1,3,'uin',1,0) crashed scilab.
+if execstr("grand(1,3,''uin'',1,0)", "errcatch") <> 999 then bugmes();quit;end
+if execstr("grand(1,3,''unf'',1,0)", "errcatch") <> 999 then bugmes();quit;end
diff --git a/modules/randlib/tests/nonreg_tests/bug_9159.tst b/modules/randlib/tests/nonreg_tests/bug_9159.tst
new file mode 100755
index 000000000..870bdee4e
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_9159.tst
@@ -0,0 +1,19 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009 - Digiteo - Allan CORNET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 9159 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=9159
+//
+// <-- Short Description -->
+// grand(1,3,'uin',1,0) crashed scilab.
+
+if execstr("grand(1,3,''uin'',1,0)", "errcatch") <> 999 then pause, end
+if execstr("grand(1,3,''unf'',1,0)", "errcatch") <> 999 then pause, end
diff --git a/modules/randlib/tests/nonreg_tests/bug_9584.dia.ref b/modules/randlib/tests/nonreg_tests/bug_9584.dia.ref
new file mode 100755
index 000000000..70bb60816
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_9584.dia.ref
@@ -0,0 +1,22 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 9584 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=9584
+//
+// <-- Short Description -->
+// grand returns non-empty matrix for negative size argument
+errMsg1 = msprintf(_("%s: Wrong value for input argument #%d: Positive scalar expected.\n"),"grand",1);
+errMsg2 = msprintf(_("%s: Wrong value for input argument #%d: Positive scalar expected.\n"),"grand",2);
+assert_checkequal(grand(1,0,'def'),[]);
+assert_checkerror("grand(-1,1,""def"");", errMsg1);
+assert_checkerror("grand(1,-1,""def"");", errMsg2);
+assert_checkequal(size(grand(-1,-1,'def')),[-1,-1]); //size(eye())==[-1 -1]
diff --git a/modules/randlib/tests/nonreg_tests/bug_9584.tst b/modules/randlib/tests/nonreg_tests/bug_9584.tst
new file mode 100755
index 000000000..51477659e
--- /dev/null
+++ b/modules/randlib/tests/nonreg_tests/bug_9584.tst
@@ -0,0 +1,23 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 9584 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=9584
+//
+// <-- Short Description -->
+// grand returns non-empty matrix for negative size argument
+
+errMsg1 = msprintf(_("%s: Wrong value for input argument #%d: Positive scalar expected.\n"),"grand",1);
+errMsg2 = msprintf(_("%s: Wrong value for input argument #%d: Positive scalar expected.\n"),"grand",2);
+assert_checkequal(grand(1,0,'def'),[]);
+assert_checkerror("grand(-1,1,""def"");", errMsg1);
+assert_checkerror("grand(1,-1,""def"");", errMsg2);
+assert_checkequal(size(grand(-1,-1,'def')),[-1,-1]); //size(eye())==[-1 -1]
diff --git a/modules/randlib/tests/unit_tests/grand_clcg4.dia.ref b/modules/randlib/tests/unit_tests/grand_clcg4.dia.ref
new file mode 100755
index 000000000..413af850e
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_clcg4.dia.ref
@@ -0,0 +1,115 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+grand('setgen','clcg4');
+// Low level test for grand
+//**********************************************************************
+//     A test program for the bottom level routines
+//     Scilab implementation of tstbot.f
+//**********************************************************************
+NB       = 10
+ NB  =
+ 
+    10.  
+NR       = 1000
+ NR  =
+ 
+    1000.  
+answer   = ones(NB,NR);
+genlst   = [1,5,10,20,32]
+ genlst  =
+ 
+    1.    5.    10.    20.    32.  
+nbad     = 0;
+str      = ['For five virtual generators of the 101';..
+			' This test generates'+string(NB)+' numbers then resets the block';..
+			'      and does it again';..
+			' Any disagreements are reported -- there should be none'];
+			
+write (%io(2),str);
+For five virtual generators of the 101
+ This test generates10 numbers then resets the block
+      and does it again
+ Any disagreements are reported -- there should be none
+//
+//     Set up Generators
+//
+grand('setall',12345,54321,6789,9876);
+//
+//     For a selected set of generators
+//
+for ixgen = 1:5
+	igen = genlst(ixgen)
+	grand('setcgn',igen);
+	write(%io(2),' Testing generator '+string(igen));
+	//
+	//     Use NB blocks
+	//
+	grand('initgn',-1);
+	SD=grand('getsd');iseed1=SD(1);iseed2=SD(2);
+	for iblock = 1:NB
+		//     Generate NR numbers
+		answer(iblock,1:NR)= grand(1,NR,'lgi');
+		grand('initgn',1);
+	end
+	grand('initgn',-1);
+	//
+	//     Do it again and compare answers
+	//
+	SD=grand('getsd');iseed1=SD(1);iseed2=SD(2);
+	//
+	//     Use NB blocks
+	//
+	for iblock = 1:NB
+		
+		//     Generate NR numbers
+		itmp = grand(1,NR,'lgi');
+		if itmp<>answer(iblock,:) then
+			str=[' Disagreement on regeneration of numbers';...
+			' Block '+string(iblock)+' N within Block '];
+			write(%io(2),str);
+		end
+		
+		if itmp<>answer(iblock,:) then bugmes();quit;end
+		
+		grand('initgn',1);
+		
+	end
+	
+	write (%io(2), ' Finished testing generator '+string(igen));
+	write (%io(2), ' Test completed successfully');
+	
+end
+ igen  =
+ 
+    1.  
+ Testing generator 1
+ Finished testing generator 1
+ Test completed successfully
+ igen  =
+ 
+    5.  
+ Testing generator 5
+ Finished testing generator 5
+ Test completed successfully
+ igen  =
+ 
+    10.  
+ Testing generator 10
+ Finished testing generator 10
+ Test completed successfully
+ igen  =
+ 
+    20.  
+ Testing generator 20
+ Finished testing generator 20
+ Test completed successfully
+ igen  =
+ 
+    32.  
+ Testing generator 32
+ Finished testing generator 32
+ Test completed successfully
diff --git a/modules/randlib/tests/unit_tests/grand_clcg4.tst b/modules/randlib/tests/unit_tests/grand_clcg4.tst
new file mode 100755
index 000000000..c7f3e76e7
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_clcg4.tst
@@ -0,0 +1,80 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+grand('setgen','clcg4');
+
+// Low level test for grand
+
+//**********************************************************************
+//     A test program for the bottom level routines
+//     Scilab implementation of tstbot.f
+//**********************************************************************
+NB       = 10
+NR       = 1000
+answer   = ones(NB,NR);
+genlst   = [1,5,10,20,32]
+nbad     = 0;
+
+str      = ['For five virtual generators of the 101';..
+			' This test generates'+string(NB)+' numbers then resets the block';..
+			'      and does it again';..
+			' Any disagreements are reported -- there should be none'];
+			
+write (%io(2),str);
+
+//
+//     Set up Generators
+//
+
+grand('setall',12345,54321,6789,9876);
+
+//
+//     For a selected set of generators
+//
+
+for ixgen = 1:5
+	igen = genlst(ixgen)
+	grand('setcgn',igen);
+	write(%io(2),' Testing generator '+string(igen));
+	//
+	//     Use NB blocks
+	//
+	grand('initgn',-1);
+	SD=grand('getsd');iseed1=SD(1);iseed2=SD(2);
+	for iblock = 1:NB
+		//     Generate NR numbers
+		answer(iblock,1:NR)= grand(1,NR,'lgi');
+		grand('initgn',1);
+	end
+	grand('initgn',-1);
+	//
+	//     Do it again and compare answers
+	//
+	SD=grand('getsd');iseed1=SD(1);iseed2=SD(2);
+	//
+	//     Use NB blocks
+	//
+	for iblock = 1:NB
+		
+		//     Generate NR numbers
+		itmp = grand(1,NR,'lgi');
+		if itmp<>answer(iblock,:) then
+			str=[' Disagreement on regeneration of numbers';...
+			' Block '+string(iblock)+' N within Block '];
+			write(%io(2),str);
+		end
+		
+		if itmp<>answer(iblock,:) then pause,end
+		
+		grand('initgn',1);
+		
+	end
+	
+	write (%io(2), ' Finished testing generator '+string(igen));
+	write (%io(2), ' Test completed successfully');
+	
+end
diff --git a/modules/randlib/tests/unit_tests/grand_generators.dia.ref b/modules/randlib/tests/unit_tests/grand_generators.dia.ref
new file mode 100755
index 000000000..5810b2751
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_generators.dia.ref
@@ -0,0 +1,499 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- ENGLISH IMPOSED -->
+function checkMeanVariance2arg ( m , n , name , A , B , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  R=grand(m,n,name,A,B);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+function checkMeanVariance0arg ( m , n , name , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  R=grand(m,n,name);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+function checkLaw0arg ( name , cdffun , N , NC , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  R = grand(1,N,name);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X)
+  TheoricPDF = diff(CDF);
+  assert_checktrue ( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+function checkPieceLaw0arg ( name , cdffun , N , NC , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+  R=grand(N,1,name);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+function checkPieceLaw2arg ( name , cdffun , N , NC , A , B , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+  R=grand(N,1,name,A,B);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X,A,B)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+function [x,y] = histcompute(n,data)
+   //
+   // Computes the histogram of a data.
+   //
+   // Parameters
+   // n : a 1-by-1 matrix of floating point integers, the number of classes.
+   // data : a matrix of doubles, the data
+   // x : 1-by-n+1 matrix of doubles, the classes of the histogram
+   // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
+   x = linspace(min(data), max(data), n+1);
+   [ind , y] = dsearch(data, x)
+   y = y./length(data)
+endfunction
+function p = mycdfdef (X)
+  p = X
+endfunction
+function p = mycdfuin (X,A,B)
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+function p = mycdflgi (X)
+  // Here, A and B are read from the environment
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+//
+// These tests makes checks that all uniform generators are correct.
+// We run the following actions : 'setgen', 'getgen', 'setsd','getsd'
+// We also check the first ten uniform numbers from each generator with a known seed.
+// seedsize : size of the state for each generator
+// MinInt : minimum of the uniform integer interval for random number generation
+// MaxInt : maximum of the uniform integer interval for random number generation
+//
+NGen = 6;
+generators = ["mt"   "kiss" "clcg2"    "clcg4" "fsultra" "urand"];
+seedsize =   [625    4      2          4       40        1];
+MinInt =     [0      0      0          0       0         0];
+MaxInt =     [2^32-1 2^32-1 2147483561 2^31-2  2^32-1    2^31-1];
+rtol = 1.e-2;
+// The number of classes in the histogram
+// NC must be even.
+NC = 2*12;
+N=10000;
+//
+// The default generator must be Mersenne-Twister
+S=grand('getgen');
+assert_checkequal ( S , "mt" );
+// The maximum integer generable with uin option
+UinMax = 2147483560;
+////////////////////////////////////////////////////////////////////
+//
+// "mt"
+//
+kgen = 1;
+gen = "mt";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',0);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.5488135    0.6027634    0.4236548    0.4375872    0.9636628    0.7917250  
+0.5928446    0.8579456    0.6235637    0.2975346    0.2726563    0.8121687  
+0.7151894    0.5448832    0.6458941    0.891773     0.3834415    0.5288949  
+0.8442657    0.8472517    0.3843817    0.0567130    0.4776651    0.4799772  
+];
+grand('setsd',0);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    2357136044.    2588848963.    1819583497.    1879422756.    4138900056.    3400433126.  
+    2546248239.    3684848379.    2678185683.    1277901399.    1171049868.    3488238119.  
+    3071714933.    2340255427.    2774094101.    3830135878.    1646868794.    2271586391.  
+    3626093760.    3638918503.    1650906866.    243580376.     2051556033.    2061486254.  
+];
+grand('setsd',0);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "kiss"
+//
+kgen = 2;
+gen = "kiss";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',0,0,0,0);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456,123456,123456,123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+    2.874450D-04    9.423555D-01    7.707249D-02    2.078324D-02    4.746445D-01    1.895302D-01  
+    8.538282D-01    5.493145D-01    3.200836D-01    4.775516D-01    2.245108D-01    6.637360D-01  
+    5.815227D-02    6.006782D-01    8.569004D-01    1.235649D-02    7.357421D-01    5.837571D-01  
+    5.196679D-01    2.448867D-01    2.568304D-01    4.503826D-01    9.680347D-01    5.214808D-01  
+];
+grand('setsd',0,0,0,0);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    1234567.       4047385867.    331023823.     89263315.      2038582807.    814026139.   
+    3667164066.    2359287638.    1374748746.    2051068542.    964266482.     2850724518.  
+    249762113.     2579893349.    3680359369.    53070701.      3159988049.    2507217781.  
+    2231956628.    1051780200.    1103078268.    1934378448.    4157677412.    2239743032.  
+];
+grand('setsd',0,0,0,0);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "clcg2"
+//
+kgen = 3;
+gen = "clcg2";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',1,1);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456,123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.9999997    0.0369445    0.2041364    0.9100817    0.6998243    0.9596867  
+0.9745196    0.1617119    0.1673069    0.1117162    0.9502824    0.9149753  
+0.6474839    0.6646450    0.6549574    0.2990212    0.0918107    0.4411791  
+0.3330856    0.0846729    0.1288161    0.2654475    0.9023415    0.0735483  
+];
+grand('setsd',1,1);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-5 );
+//
+// Check integers
+expected = [
+    2147482884.    79337801.      438379562.     1954385533.    1502861091.    2060911403.  
+    2092764894.    347273588.     359288887.     239908781.     2040715732.    1964894377.  
+    1390461064.    1427314282.    1406510214.    642143055.     197161966.     947424871.   
+    715295839.     181833602.     276630551.     570044126.     1937763493.    157943826.   
+];
+grand('setsd',1,1);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "clcg4"
+//
+kgen = 4;
+gen = "clcg4";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+warning("off");
+grand('setsd',1,1,1,1);
+warning("on");
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+warning("off");
+grand('setsd',123456,123456,123456,123456);
+warning("on");
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.9999661    0.0552914    0.6345306    0.0640227    0.2885048    0.2781458  
+0.6852419    0.1806991    0.8665501    0.0981421    0.2660715    0.4279616  
+0.4370514    0.4956021    0.6870544    0.8501209    0.1271038    0.4554926  
+0.4202952    0.2903676    0.5712601    0.4764120    0.1818799    0.3121748  
+];
+warning("off");
+grand('setsd',1,1,1,1);
+warning("on");
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    2147410798.    118737467.     1362644105.    137487708.     619559402.    597313629.  
+    1471545799.    388048305.     1860902184.    210758531.     571384155.    919040470.  
+    938560647.     1064297484.    1475437993.    1825620628.    272953383.    978162913.  
+    902576998.     623559632.     1226771622.    1023086907.    390584072.    670390361.  
+];
+warning("off");
+grand('setsd',1,1,1,1);
+warning("on");
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "fsultra"
+//
+kgen = 5;
+gen = "fsultra";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',1,1);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456,123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.3314877    0.3699260    0.4383216    0.99706      0.0577929    0.4836669  
+0.5826624    0.9600475    0.2037475    0.6774254    0.4450278    0.3082941  
+0.1630857    0.2033307    0.4214824    0.6372521    0.0782678    0.4409892  
+0.7211611    0.1833922    0.8065496    0.6375251    0.2572713    0.8039582  
+];
+grand('setsd',1,1);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    1423728774.    1588820113.    1882577034.    4282340079.    248218608.     2077333598.  
+    2502516061.    4123372468.    875088704.     2909519830.    1911379739.    1324113135.  
+    700447838.     873298853.     1810253313.    2736976953.    336157762.     1894034123.  
+    3097363227.    787663378.     3464104206.    2738149622.    1104971606.    3452974260.  
+];
+grand('setsd',1,1);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "urand"
+//
+kgen = 6;
+gen = "urand";
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',1);
+S=grand('getsd');
+assert_checkequal ( S , 1 );
+//
+grand('setsd',123456);
+S=grand('getsd');
+assert_checkequal ( S , 123456 );
+//
+// Check numbers
+expected = [
+0.6040239    0.5321420    0.2276811    0.8979351    0.8925854    0.6928366  
+0.0079647    0.4138784    0.6656067    0.8274020    0.0848163    0.6776849  
+0.6643966    0.5036204    0.9694369    0.0664231    0.2566682    0.4284010  
+0.9832111    0.6850569    0.0775390    0.1099766    0.6507795    0.3725794  
+];
+grand('setsd',1);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-5 );
+//
+// Check integers
+expected = [
+    1297131554.    1142766270.    488941338.     1928300854.    1916812562.    1487855278.  
+    17103983.      888797147.     1429379591.    1776832243.    182141599.     1455317259.  
+    1426780792.    1081516660.    2081849904.    142642604.     551190760.     919984100.   
+    2111429773.    1471148505.    166513637.     236172977.     1397538365.    800108169.   
+];
+grand('setsd',1);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
diff --git a/modules/randlib/tests/unit_tests/grand_generators.tst b/modules/randlib/tests/unit_tests/grand_generators.tst
new file mode 100755
index 000000000..0c73d3bdf
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_generators.tst
@@ -0,0 +1,522 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// <-- ENGLISH IMPOSED -->
+
+function checkMeanVariance2arg ( m , n , name , A , B , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  
+  R=grand(m,n,name,A,B);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+
+function checkMeanVariance0arg ( m , n , name , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  
+  R=grand(m,n,name);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+
+function checkLaw0arg ( name , cdffun , N , NC , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+
+  R = grand(1,N,name);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X)
+  TheoricPDF = diff(CDF);
+  assert_checktrue ( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+
+function checkPieceLaw0arg ( name , cdffun , N , NC , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+
+  R=grand(N,1,name);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+
+function checkPieceLaw2arg ( name , cdffun , N , NC , A , B , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+
+  R=grand(N,1,name,A,B);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X,A,B)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+
+function [x,y] = histcompute(n,data)
+   //
+   // Computes the histogram of a data.
+   //
+   // Parameters
+   // n : a 1-by-1 matrix of floating point integers, the number of classes.
+   // data : a matrix of doubles, the data
+   // x : 1-by-n+1 matrix of doubles, the classes of the histogram
+   // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
+   x = linspace(min(data), max(data), n+1);
+   [ind , y] = dsearch(data, x)
+   y = y./length(data)
+endfunction
+
+function p = mycdfdef (X)
+  p = X
+endfunction
+
+function p = mycdfuin (X,A,B)
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+
+function p = mycdflgi (X)
+  // Here, A and B are read from the environment
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+
+
+//
+// These tests makes checks that all uniform generators are correct.
+// We run the following actions : 'setgen', 'getgen', 'setsd','getsd'
+// We also check the first ten uniform numbers from each generator with a known seed.
+// seedsize : size of the state for each generator
+// MinInt : minimum of the uniform integer interval for random number generation
+// MaxInt : maximum of the uniform integer interval for random number generation
+//
+NGen = 6;
+generators = ["mt"   "kiss" "clcg2"    "clcg4" "fsultra" "urand"];
+seedsize =   [625    4      2          4       40        1];
+MinInt =     [0      0      0          0       0         0];
+MaxInt =     [2^32-1 2^32-1 2147483561 2^31-2  2^32-1    2^31-1];
+
+rtol = 1.e-2;
+
+// The number of classes in the histogram
+// NC must be even.
+NC = 2*12;
+N=10000;
+
+//
+// The default generator must be Mersenne-Twister
+S=grand('getgen');
+assert_checkequal ( S , "mt" );
+
+// The maximum integer generable with uin option
+UinMax = 2147483560;
+
+////////////////////////////////////////////////////////////////////
+//
+// "mt"
+//
+kgen = 1;
+gen = "mt";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',0);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.5488135    0.6027634    0.4236548    0.4375872    0.9636628    0.7917250  
+0.5928446    0.8579456    0.6235637    0.2975346    0.2726563    0.8121687  
+0.7151894    0.5448832    0.6458941    0.891773     0.3834415    0.5288949  
+0.8442657    0.8472517    0.3843817    0.0567130    0.4776651    0.4799772  
+];
+grand('setsd',0);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    2357136044.    2588848963.    1819583497.    1879422756.    4138900056.    3400433126.  
+    2546248239.    3684848379.    2678185683.    1277901399.    1171049868.    3488238119.  
+    3071714933.    2340255427.    2774094101.    3830135878.    1646868794.    2271586391.  
+    3626093760.    3638918503.    1650906866.    243580376.     2051556033.    2061486254.  
+];
+grand('setsd',0);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+
+////////////////////////////////////////////////////////////////////
+//
+// "kiss"
+//
+kgen = 2;
+gen = "kiss";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',0,0,0,0);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456,123456,123456,123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+    2.874450D-04    9.423555D-01    7.707249D-02    2.078324D-02    4.746445D-01    1.895302D-01  
+    8.538282D-01    5.493145D-01    3.200836D-01    4.775516D-01    2.245108D-01    6.637360D-01  
+    5.815227D-02    6.006782D-01    8.569004D-01    1.235649D-02    7.357421D-01    5.837571D-01  
+    5.196679D-01    2.448867D-01    2.568304D-01    4.503826D-01    9.680347D-01    5.214808D-01  
+];
+grand('setsd',0,0,0,0);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    1234567.       4047385867.    331023823.     89263315.      2038582807.    814026139.   
+    3667164066.    2359287638.    1374748746.    2051068542.    964266482.     2850724518.  
+    249762113.     2579893349.    3680359369.    53070701.      3159988049.    2507217781.  
+    2231956628.    1051780200.    1103078268.    1934378448.    4157677412.    2239743032.  
+];
+grand('setsd',0,0,0,0);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "clcg2"
+//
+kgen = 3;
+gen = "clcg2";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',1,1);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456,123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.9999997    0.0369445    0.2041364    0.9100817    0.6998243    0.9596867  
+0.9745196    0.1617119    0.1673069    0.1117162    0.9502824    0.9149753  
+0.6474839    0.6646450    0.6549574    0.2990212    0.0918107    0.4411791  
+0.3330856    0.0846729    0.1288161    0.2654475    0.9023415    0.0735483  
+];
+grand('setsd',1,1);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-5 );
+//
+// Check integers
+expected = [
+    2147482884.    79337801.      438379562.     1954385533.    1502861091.    2060911403.  
+    2092764894.    347273588.     359288887.     239908781.     2040715732.    1964894377.  
+    1390461064.    1427314282.    1406510214.    642143055.     197161966.     947424871.   
+    715295839.     181833602.     276630551.     570044126.     1937763493.    157943826.   
+];
+grand('setsd',1,1);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "clcg4"
+//
+kgen = 4;
+gen = "clcg4";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+warning("off");
+grand('setsd',1,1,1,1);
+warning("on");
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+warning("off");
+grand('setsd',123456,123456,123456,123456);
+warning("on");
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.9999661    0.0552914    0.6345306    0.0640227    0.2885048    0.2781458  
+0.6852419    0.1806991    0.8665501    0.0981421    0.2660715    0.4279616  
+0.4370514    0.4956021    0.6870544    0.8501209    0.1271038    0.4554926  
+0.4202952    0.2903676    0.5712601    0.4764120    0.1818799    0.3121748  
+];
+warning("off");
+grand('setsd',1,1,1,1);
+warning("on");
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    2147410798.    118737467.     1362644105.    137487708.     619559402.    597313629.  
+    1471545799.    388048305.     1860902184.    210758531.     571384155.    919040470.  
+    938560647.     1064297484.    1475437993.    1825620628.    272953383.    978162913.  
+    902576998.     623559632.     1226771622.    1023086907.    390584072.    670390361.  
+];
+warning("off");
+grand('setsd',1,1,1,1);
+warning("on");
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "fsultra"
+//
+kgen = 5;
+gen = "fsultra";
+sdsize = seedsize(kgen);
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',1,1);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+grand('setsd',123456,123456);
+S=grand('getsd');
+assert_checkequal ( typeof(S) , "constant" );
+assert_checkequal ( size(S) , [sdsize 1] );
+//
+// Check numbers
+expected = [
+0.3314877    0.3699260    0.4383216    0.99706      0.0577929    0.4836669  
+0.5826624    0.9600475    0.2037475    0.6774254    0.4450278    0.3082941  
+0.1630857    0.2033307    0.4214824    0.6372521    0.0782678    0.4409892  
+0.7211611    0.1833922    0.8065496    0.6375251    0.2572713    0.8039582  
+];
+grand('setsd',1,1);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-6 );
+//
+// Check integers
+expected = [
+    1423728774.    1588820113.    1882577034.    4282340079.    248218608.     2077333598.  
+    2502516061.    4123372468.    875088704.     2909519830.    1911379739.    1324113135.  
+    700447838.     873298853.     1810253313.    2736976953.    336157762.     1894034123.  
+    3097363227.    787663378.     3464104206.    2738149622.    1104971606.    3452974260.  
+];
+grand('setsd',1,1);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+////////////////////////////////////////////////////////////////////
+//
+// "urand"
+//
+kgen = 6;
+gen = "urand";
+grand('setgen',gen);
+S=grand('getgen');
+assert_checkequal ( S , gen );
+//
+grand('setsd',1);
+S=grand('getsd');
+assert_checkequal ( S , 1 );
+//
+grand('setsd',123456);
+S=grand('getsd');
+assert_checkequal ( S , 123456 );
+//
+// Check numbers
+expected = [
+0.6040239    0.5321420    0.2276811    0.8979351    0.8925854    0.6928366  
+0.0079647    0.4138784    0.6656067    0.8274020    0.0848163    0.6776849  
+0.6643966    0.5036204    0.9694369    0.0664231    0.2566682    0.4284010  
+0.9832111    0.6850569    0.0775390    0.1099766    0.6507795    0.3725794  
+];
+grand('setsd',1);
+computed = grand(4,6,"def");
+assert_checkalmostequal ( computed , expected, 1.e-5 );
+//
+// Check integers
+expected = [
+    1297131554.    1142766270.    488941338.     1928300854.    1916812562.    1487855278.  
+    17103983.      888797147.     1429379591.    1776832243.    182141599.     1455317259.  
+    1426780792.    1081516660.    2081849904.    142642604.     551190760.     919984100.   
+    2111429773.    1471148505.    166513637.     236172977.     1397538365.    800108169.   
+];
+grand('setsd',1);
+computed = grand(4,6,"lgi");
+assert_checkequal ( computed , expected );
+//
+// Check distribution of uniform numbers in [0,1[
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+//
+// Check distribution of uniform integers in [A,B]
+A = MinInt(kgen);
+B = MaxInt(kgen);
+checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
+checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
+
diff --git a/modules/randlib/tests/unit_tests/grand_hypermat.dia.ref b/modules/randlib/tests/unit_tests/grand_hypermat.dia.ref
new file mode 100755
index 000000000..2995e1005
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_hypermat.dia.ref
@@ -0,0 +1,98 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010-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('randlib', 'grand_hypermat', ['no_check_error_output'])
+///////////////////////////////////////////////////////////////////////////////
+//
+// Dimensions
+mat = grand(100, 101, 102, 'unf', 0, 1);
+assert_checktrue(size(mat) == [100 101 102]);
+///////////////////////////////////////////////////////////////////////////////
+//
+// Generators
+// The grand(i, j, ...) should be equal to the first element of grand(i, j, k, ...).
+// mt generator
+grand('setgen', 'mt');
+grand('setsd', 0);
+expected = grand(4, 6, 'def');
+grand('setsd', 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 0);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+// kiss generator
+grand('setgen', 'kiss');
+grand('setsd', 0, 0, 0, 0);
+expected = grand(4, 6, 'def');
+grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 0, 0, 0, 0);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+// clcg2 generator
+grand('setgen', 'clcg2');
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'def');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+// clcg4 generator
+grand('setgen', 'clcg4');
+warning('off');
+grand('setsd',1,1,1,1);
+warning('on');
+expected = grand(4, 6, 'def');
+warning('off');
+grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
+warning('on');
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+warning('off');
+grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
+warning('on');
+expected = grand(4, 6, 'lgi');
+warning('off');
+grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
+warning('on');
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+// fsultra generator
+grand('setgen', 'fsultra');
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'def');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+// urand generator
+grand('setgen', 'urand');
+grand('setsd', 1);
+expected = grand(4, 6, 'def');
+grand('setsd', 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 1);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
diff --git a/modules/randlib/tests/unit_tests/grand_hypermat.tst b/modules/randlib/tests/unit_tests/grand_hypermat.tst
new file mode 100755
index 000000000..08702c6a3
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_hypermat.tst
@@ -0,0 +1,109 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010-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('randlib', 'grand_hypermat', ['no_check_error_output'])
+
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// Dimensions
+mat = grand(100, 101, 102, 'unf', 0, 1);
+assert_checktrue(size(mat) == [100 101 102]);
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// Generators
+// The grand(i, j, ...) should be equal to the first element of grand(i, j, k, ...).
+
+// mt generator
+grand('setgen', 'mt');
+grand('setsd', 0);
+expected = grand(4, 6, 'def');
+grand('setsd', 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 0);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+
+// kiss generator
+grand('setgen', 'kiss');
+grand('setsd', 0, 0, 0, 0);
+expected = grand(4, 6, 'def');
+grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 0, 0, 0, 0);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+
+// clcg2 generator
+grand('setgen', 'clcg2');
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'def');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+
+// clcg4 generator
+grand('setgen', 'clcg4');
+warning('off');
+grand('setsd',1,1,1,1);
+warning('on');
+expected = grand(4, 6, 'def');
+warning('off');
+grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
+warning('on');
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+warning('off');
+grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
+warning('on');
+expected = grand(4, 6, 'lgi');
+warning('off');
+grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
+warning('on');
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+
+// fsultra generator
+grand('setgen', 'fsultra');
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'def');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 1, 1);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
+
+// urand generator
+grand('setgen', 'urand');
+grand('setsd', 1);
+expected = grand(4, 6, 'def');
+grand('setsd', 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'def');
+assert_checkequal(expected, computed(:, :, 1));
+grand('setsd', 1);
+expected = grand(4, 6, 'lgi');
+grand('setsd', 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, 'lgi');
+assert_checkequal(expected, computed(:, :, 1));
diff --git a/modules/randlib/tests/unit_tests/grand_laws.dia.ref b/modules/randlib/tests/unit_tests/grand_laws.dia.ref
new file mode 100755
index 000000000..bbbf1562c
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_laws.dia.ref
@@ -0,0 +1,594 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+function [x,y] = histcompute(n,data)
+   //
+   // Computes the histogram of a data.
+   //
+   // Parameters
+   // n : a 1-by-1 matrix of floating point integers, the number of classes.
+   // data : a matrix of doubles, the data
+   // x : 1-by-n+1 matrix of doubles, the classes of the histogram
+   // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
+   x = linspace(min(data), max(data), n+1);
+   [ind , y] = dsearch(data, x)
+   y = y./length(data)
+endfunction
+function checkLaw0arg ( name , cdffun , N , NC , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  R = grand(1,N,name);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+function checkLaw1arg ( name , cdffun , N , NC , A , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  R = grand(1,N,name,A);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X,A)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+function checkLaw2arg ( name , cdffun , N , NC , A , B , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  R = grand(1,N,name,A,B);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X,A,B)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+function checkLaw3arg ( name , cdffun , N , NC , A , B , C , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // C : a 1-by-1 matrix of doubles, the value of the 3d parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  R = grand(1,N,name,A,B,C);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X,A,B,C)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+function checkPieceLaw1arg ( name , cdffun , N , NC , A , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+  R=grand(N,1,name,A);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X,A)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+function checkPieceLaw2arg ( name , cdffun , N , NC , A , B , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+  R=grand(N,1,name,A,B);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X,A,B)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+function checkMeanVariance0arg ( m , n , name , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  R=grand(m,n,name);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+function checkMeanVariance1arg ( m , n , name , A , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  R=grand(m,n,name,A);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+function checkMeanVariance2arg ( m , n , name , A , B , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  R=grand(m,n,name,A,B);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+function checkMeanVariance3arg ( m , n , name , A , B , C , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // C : a 1-by-1 matrix of doubles, the value of the 3d parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  R=grand(m,n,name,A,B,C);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+//
+// Calling sequences
+//
+// Y=grand(m,n,'nor',Av,Sd)       : OK
+// Y=grand(m,n,'poi',mu)          : OK
+// Y=grand(m,n,'exp',Av)          : OK
+// Y=grand(m,n,'bet',A,B)         : OK
+// Y=grand(m,n,'bin',N,p)         : OK
+// Y=grand(m,n,'gam',shape,scale) : OK
+// Y=grand(m,n,'f',Dfn,Dfd)       : OK
+// Y=grand(m,n,'chi', Df)         : OK
+// Y=grand(m,n,'nbn',N,p)         : OK
+// Y=grand(m,n,'nch',Df,Xnon)     : OK
+// Y=grand(m,n,'nf',Dfn,Dfd,Xnon) : OK
+// Y=grand(n,'mn',Mean,Cov)       : OK
+// Y=grand(m,n,'geom', p)         : OK
+// Y=grand(m,n,'def')             : OK
+// Y=grand(m,n,'unf',Low,High)    : OK
+// Y=grand(m,n,'uin',Low,High)    : OK
+// Y=grand(n,'markov',P,x0)       : TODO
+// Y=grand(n,'mul',nb,P)          : TODO
+// Y=grand(n,'prm',vect)          : in grand_prm.tst
+// Y=grand(m,n,'lgi')  : in grand_generators.tst
+// S=grand('getgen')   : in grand_generators.tst
+// grand('setgen',gen) : in grand_generators.tst
+// S=grand('getsd')    : in grand_generators.tst
+// grand('setsd',S)    : in grand_generators.tst
+// For CLCG4 only:
+// grand('setcgn',G)           : in grand_clcg4.tst
+// S=grand('getcgn')           : TODO
+// grand('initgn',I)           : in grand_clcg4.tst
+// grand('setall',s1,s2,s3,s4) : in grand_clcg4.tst
+// grand('advnst',K)           : TODO
+//
+// Rationale for the unit tests of random number generators.
+//
+// There are several ways to check the output of random number generators.
+//
+// 1. Uniform random number generators are based on deterministic integer sequences.
+// These integers are then scaled into the [0,1[ interval.
+// When we set the seed to a constant (e.g. seed = 0, or seed=123456),
+// the doubles which are produced are always the same.
+// This does not check the quality of the generator.
+// But it does check that the generator under evaluation is, indeed, the 
+// expected generator.
+//
+// 2. Non-uniform random numbers have a given mean and a variance which 
+// can be predicted by there associated distribution function.
+// The problem is there to set the relative tolerance. 
+// The convergence of the sample mean to the expectation is driven by the
+// law of large numbers. 
+// The standard deviation of the sample mean is sigma/sqrt(n), 
+// meaning that to reduce the shift by a factor 2, we have to use 4x more points. 
+//
+// 3. The quality of a random number generator can be assessed by 
+// statistical tests. 
+// The Kolmogorov-Smirnov test is based on the comparison between the 
+// empirical CDF and the theoretical CDF. The difference between the two is 
+// a random variable, which must evolve according to the Normal distribution function. 
+// Then, the sum of the squares of the differences should evolve according to the 
+// Chi-square distribution function. By inverting this last distribution, we can 
+// compute the probability that the empirical distribution matches the theoretical distribution.
+//
+// 4. The empirical histogram must match the theoretical histogram. 
+// The number of classes is NC.
+// The classes XC are computed from NC and the actual output ; X has NC+1 entries.
+// The empirical histogram is the probability that X is the intervals defined by X.
+// The theoretical histogram is computed by differences 
+// of the cumulated probability distribution function.
+// The histograms can be graphically compared with the statements:
+// plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+// plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+//
+//
+// In this script, we check the non-uniform random numbers.
+// We use the strategies #2 and #4.
+// The strategy #3 should be used in a future version of the tests (it requires to 
+// develop the kstest function).
+//
+rtol = 1.e-2;
+// The number of classes in the histogram
+// NC must be even.
+NC = 2*12;
+N=10000;
+// Set the seed to always get the same random numbers
+grand("setsd",0);
+/////////////////////////////////////////////////////////////
+//
+// 'exp'
+//
+function p = mycdfexp(X,A)
+  p = 1 - exp(-(1/A)*X);
+endfunction
+A = 20;
+checkMeanVariance1arg ( 400 , 800 , 'exp' , A , A , A^2 , 2*rtol );
+//
+for A = linspace(0.1,50,10)
+  checkLaw1arg ( "exp" , mycdfexp , N , NC , A , 2*rtol );
+end
+/////////////////////////////////////////////////////////////
+//
+// 'gam'
+//
+function p = mycdfgam (X,A,B)
+  p = cdfgam("PQ",X,A*ones(X),B*ones(X));
+endfunction
+A = 5;
+B = 3;
+checkMeanVariance2arg ( 400 , 800 , 'gam' , A , B , A/B , A/(B^2) , rtol );
+//
+for A = linspace(1,15,4)
+  for B = linspace(1,15,4)
+    checkLaw2arg ( 'gam' , mycdfgam , N , NC , A , B , rtol );
+  end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'bet'
+//
+function p = mycdfbet (X,A,B)
+  p = cdfbet("PQ",X,1-X,A*ones(X),B*ones(X));
+endfunction
+A=3;
+B=10;
+checkMeanVariance2arg ( 400 , 800 , 'bet' , A , B , A/(A+B) , (A*B)/((A+B)^2*(A+B+1)) , rtol );
+//
+for A = linspace(1,20,4)
+  for B = linspace(1,20,4)
+    checkLaw2arg ( 'bet' , mycdfbet , N , NC , A , B , 2*rtol );
+  end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'poi'
+//
+// Caution : this is a piecewise distribution
+function p = mycdfpoi (X,A)
+  p = cdfpoi("PQ",X,A*ones(X));
+endfunction
+A = 10;
+checkMeanVariance1arg ( 400 , 800 , 'poi' , A , A , A , rtol );
+//
+for A = floor(linspace(50,70,10))
+  checkPieceLaw1arg ( 'poi' , mycdfpoi , N , NC , A , rtol );
+end
+/////////////////////////////////////////////////////////////
+//
+// 'bin'
+//
+// Caution : this is a piecewise distribution
+function p = mycdfbin (X,A,B)
+  p = cdfbin("PQ",X,A*ones(X),B*ones(X),(1-B)*ones(X));
+endfunction
+A = 10;
+B = 0.4;
+checkMeanVariance2arg ( 400 , 800 , 'bin' , A , B , A*B , A*B*(1-B) , rtol );
+//
+for A = floor(linspace(10,50,4))
+for B = linspace(0.1,0.9,4)
+  checkPieceLaw2arg ( 'bin' , mycdfbin , N , NC , A , B , 2*rtol );
+end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'nor'
+//
+function p = mycdfnor (X,A,B)
+  p = cdfnor("PQ",X,A*ones(X),B*ones(X));
+endfunction
+A = 7;
+B = 12;
+checkMeanVariance2arg ( 400 , 800 , 'nor' , A , B , A , B^2 , rtol );
+//
+//
+for A = linspace(1,20,4)
+  for B = linspace(1,20,4)
+    checkLaw2arg ( 'nor' , mycdfnor , N , NC , A , B , rtol );
+  end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'f'
+//
+function p = mycdff (X,A,B)
+  p = cdff("PQ",X,A*ones(X),B*ones(X));
+endfunction
+//
+// Increase the tolerance for this test.
+A = 7;
+B = 12;
+checkMeanVariance2arg ( 400 , 800 , 'f' , A , B , B/(B-2) , (2*B^2*(A+B-2))/(A*(B-2)^2*(B-4)) , 4*rtol );
+//
+for A = floor(linspace(1,20,4))
+  for B = floor(linspace(5,20,4))
+    checkLaw2arg ( 'f' , mycdff , N , NC , A , B , 2*rtol );
+  end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'chi'
+//
+function p = mycdfchi (X,A)
+  p = cdfchi("PQ",X,A*ones(X));
+endfunction
+A = 7;
+checkMeanVariance1arg ( 400 , 800 , 'chi' , A , A , 2*A , rtol );
+// Increase the tolerance for this test.
+for A = floor(linspace(1,20,10))
+    checkLaw1arg ( 'chi' , mycdfchi , N , NC , A , 2*rtol );
+end
+/////////////////////////////////////////////////////////////
+//
+// 'nbn'
+//
+// This is a piecewise function.
+function p = mycdfnbn (X,A,B)
+  p = cdfnbn("PQ",X,A*ones(X),B*ones(X),(1-B)*ones(X));
+endfunction
+A = 7;
+B = 0.1;
+checkMeanVariance2arg ( 400 , 800 , 'nbn' , A , B , A*(1-B)/B , A*(1-B)/B^2 , rtol );
+// Increase the tolerance for this test.
+for A = floor(linspace(1,20,4))
+for B = linspace(0.1,0.9,4)
+    checkPieceLaw2arg ( 'nbn' , mycdfnbn , N , NC , A , B , 2*rtol );
+end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'nch'
+//
+function p = mycdfnch (X,A,B)
+  p = cdfchn("PQ",X,A*ones(X),B*ones(X));
+endfunction
+A = 4;
+B = 3;
+checkMeanVariance2arg ( 400 , 800 , 'nch' , A , B , A+B , 2*(A+2*B) , rtol );
+for A = floor(linspace(1,20,4))
+for B = floor(linspace(1,20,4))
+    checkLaw2arg ( 'nch' , mycdfnch , N , NC , A , B , 2*rtol );
+end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'nf'
+// http://mathworld.wolfram.com/NoncentralF-Distribution.html
+//
+function p = mycdfnf (X,A,B,C)
+  p = cdffnc("PQ",X,A*ones(X),B*ones(X),C*ones(X));
+endfunction
+A = 4;
+B = 3;
+C = 10;
+checkMeanVariance3arg ( 400 , 800 , 'nf' , A , B , C , ((C+A)*B)/(A*(B-2)) , [] , 4*rtol );
+for A = floor(linspace(1,20,4))
+for B = floor(linspace(1,20,4))
+for C = floor(linspace(1,20,4))
+    checkLaw3arg ( 'nf' , mycdfnf , N , NC , A , B , C , 2*rtol );
+end
+end
+end
+/////////////////////////////////////////////////////////////
+//
+// "mn"
+//
+n = 2^16;
+A = [1;2;3];
+// The covariance matrix B is symetric, positive definite.
+// Its diagonal entries are [4;6;5], the variances.
+B = [4,2,3;2,6,4;3,4,5];
+  R=grand(n,"mn",A,B);
+  assert_checkequal ( size(R) , [3,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R,"c") , A , 2*rtol );
+  assert_checkalmostequal ( variance(R,"c") , [4;6;5] , rtol );
+//
+// No CDF for this function => no histogram test. 
+/////////////////////////////////////////////////////////////
+//
+// "geom"
+// http://en.wikipedia.org/wiki/Geometric_distribution
+//
+A = 0.1;
+checkMeanVariance1arg ( 400 , 800 , "geom" , A , 1/A , (1-A)/A^2 , rtol );
+function p = mycdfgeom (X,A)
+  p = 1 -(1-A)^X
+endfunction
+//
+for A = linspace(0.1,0.9,10)
+   checkPieceLaw1arg ( "geom" , mycdfgeom , N , NC , A , 2*rtol );
+end
+/////////////////////////////////////////////////////////////
+//
+// "unf"
+// http://en.wikipedia.org/wiki/Uniform_distribution_(continuous)
+//
+A = 0.1;
+B = 2.3;
+checkMeanVariance2arg ( 400 , 800 , "unf" , A , B , (A+B)/2 , (B-A)^2/12 , rtol );
+function p = mycdfunf (X,A,B)
+  p = (X-A)/(B-A)
+endfunction
+//
+for A = linspace(0.1,0.9,4)
+for B = linspace(2.5,4.2,4)
+   checkLaw2arg ( "unf" , mycdfunf , N , NC , A , B , rtol );
+end
+end
+/////////////////////////////////////////////////////////////
+//
+// "uin"
+// http://en.wikipedia.org/wiki/Uniform_distribution_(discrete)
+//
+// Piecewise.
+A = 10;
+B = 20;
+checkMeanVariance2arg ( 400 , 800 , "uin" , A , B , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+function p = mycdfuin (X,A,B)
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+//
+for A = floor(linspace(10,20,4))
+for B = floor(linspace(30,40,4))
+   checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , A , B , rtol );
+end
+end
+/////////////////////////////////////////////////////////////
+//
+// "def"
+// http://en.wikipedia.org/wiki/Uniform_distribution_(continuous)
+//
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+function p = mycdfdef (X)
+  p = X
+endfunction
+//
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
diff --git a/modules/randlib/tests/unit_tests/grand_laws.tst b/modules/randlib/tests/unit_tests/grand_laws.tst
new file mode 100755
index 000000000..722e0ef88
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_laws.tst
@@ -0,0 +1,677 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+
+function [x,y] = histcompute(n,data)
+   //
+   // Computes the histogram of a data.
+   //
+   // Parameters
+   // n : a 1-by-1 matrix of floating point integers, the number of classes.
+   // data : a matrix of doubles, the data
+   // x : 1-by-n+1 matrix of doubles, the classes of the histogram
+   // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
+   x = linspace(min(data), max(data), n+1);
+   [ind , y] = dsearch(data, x)
+   y = y./length(data)
+endfunction
+
+function checkLaw0arg ( name , cdffun , N , NC , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+
+  R = grand(1,N,name);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+
+function checkLaw1arg ( name , cdffun , N , NC , A , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+
+  R = grand(1,N,name,A);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X,A)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+
+function checkLaw2arg ( name , cdffun , N , NC , A , B , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+
+  R = grand(1,N,name,A,B);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X,A,B)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+
+function checkLaw3arg ( name , cdffun , N , NC , A , B , C , rtol )
+  //
+  // Check the random number generator for a continuous distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // C : a 1-by-1 matrix of doubles, the value of the 3d parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+
+  R = grand(1,N,name,A,B,C);
+  [X,EmpiricalPDF] = histcompute(NC,R);
+  CDF = cdffun(X,A,B,C)
+  TheoricPDF = diff(CDF);
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+    end
+endfunction
+
+function checkPieceLaw1arg ( name , cdffun , N , NC , A , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+
+  R=grand(N,1,name,A);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X,A)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+
+function checkPieceLaw2arg ( name , cdffun , N , NC , A , B , rtol )
+  //
+  // Check the random number generator for a piecewise distribution function.
+  //
+  // Parameters
+  // name: a 1-by-1 string, the name of the distribution function
+  // cdffun : a function, the Cumulated Distribution Function
+  // NC : a 1-by-1 matrix of floating point integers, the number of classes
+  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  //
+  // Description
+  //  Compare the empirical histogram with the theoretical histogram.
+  //  The classes of the histogram are computed from the random samples, 
+  //  from which the unique entries are extracted.
+
+  R=grand(N,1,name,A,B);
+  X = unique(R);
+  for k = 1 : size(X,"*")
+    EmpiricalPDF(k) = length(find(R==X(k)));
+  end
+  EmpiricalPDF = EmpiricalPDF./N;
+  CDF = cdffun(X,A,B)
+  TheoricPDF=[CDF(1);diff(CDF)];
+  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+  end
+endfunction
+
+function checkMeanVariance0arg ( m , n , name , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  
+  R=grand(m,n,name);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+
+
+function checkMeanVariance1arg ( m , n , name , A , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  
+  R=grand(m,n,name,A);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+
+function checkMeanVariance2arg ( m , n , name , A , B , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  
+  R=grand(m,n,name,A,B);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+  
+
+function checkMeanVariance3arg ( m , n , name , A , B , C , mu , va , rtol )
+  // Check the mean and variance of random numbers.
+  //
+  // Parameters
+  // m : a 1-by-1 matrix of floating point integers, the number of rows
+  // n : a 1-by-1 matrix of floating point integers, the number of columns
+  // name: a 1-by-1 string, the name of the distribution function
+  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+  // C : a 1-by-1 matrix of doubles, the value of the 3d parameter
+  // mu : a 1-by-1 matrix of doubles, the expected mean
+  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+  
+  R=grand(m,n,name,A,B,C);
+  assert_checkequal ( size(R) , [m,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R) , mu , rtol );
+  if ( va<>[] ) then
+    assert_checkalmostequal ( variance(R) , va , rtol );
+  end
+endfunction
+
+//
+// Calling sequences
+//
+// Y=grand(m,n,'nor',Av,Sd)       : OK
+// Y=grand(m,n,'poi',mu)          : OK
+// Y=grand(m,n,'exp',Av)          : OK
+// Y=grand(m,n,'bet',A,B)         : OK
+// Y=grand(m,n,'bin',N,p)         : OK
+// Y=grand(m,n,'gam',shape,scale) : OK
+// Y=grand(m,n,'f',Dfn,Dfd)       : OK
+// Y=grand(m,n,'chi', Df)         : OK
+// Y=grand(m,n,'nbn',N,p)         : OK
+// Y=grand(m,n,'nch',Df,Xnon)     : OK
+// Y=grand(m,n,'nf',Dfn,Dfd,Xnon) : OK
+// Y=grand(n,'mn',Mean,Cov)       : OK
+// Y=grand(m,n,'geom', p)         : OK
+// Y=grand(m,n,'def')             : OK
+// Y=grand(m,n,'unf',Low,High)    : OK
+// Y=grand(m,n,'uin',Low,High)    : OK
+// Y=grand(n,'markov',P,x0)       : TODO
+// Y=grand(n,'mul',nb,P)          : TODO
+// Y=grand(n,'prm',vect)          : in grand_prm.tst
+// Y=grand(m,n,'lgi')  : in grand_generators.tst
+// S=grand('getgen')   : in grand_generators.tst
+// grand('setgen',gen) : in grand_generators.tst
+// S=grand('getsd')    : in grand_generators.tst
+// grand('setsd',S)    : in grand_generators.tst
+
+// For CLCG4 only:
+// grand('setcgn',G)           : in grand_clcg4.tst
+// S=grand('getcgn')           : TODO
+// grand('initgn',I)           : in grand_clcg4.tst
+// grand('setall',s1,s2,s3,s4) : in grand_clcg4.tst
+// grand('advnst',K)           : TODO
+
+//
+// Rationale for the unit tests of random number generators.
+//
+// There are several ways to check the output of random number generators.
+//
+// 1. Uniform random number generators are based on deterministic integer sequences.
+// These integers are then scaled into the [0,1[ interval.
+// When we set the seed to a constant (e.g. seed = 0, or seed=123456),
+// the doubles which are produced are always the same.
+// This does not check the quality of the generator.
+// But it does check that the generator under evaluation is, indeed, the 
+// expected generator.
+//
+// 2. Non-uniform random numbers have a given mean and a variance which 
+// can be predicted by there associated distribution function.
+// The problem is there to set the relative tolerance. 
+// The convergence of the sample mean to the expectation is driven by the
+// law of large numbers. 
+// The standard deviation of the sample mean is sigma/sqrt(n), 
+// meaning that to reduce the shift by a factor 2, we have to use 4x more points. 
+//
+// 3. The quality of a random number generator can be assessed by 
+// statistical tests. 
+// The Kolmogorov-Smirnov test is based on the comparison between the 
+// empirical CDF and the theoretical CDF. The difference between the two is 
+// a random variable, which must evolve according to the Normal distribution function. 
+// Then, the sum of the squares of the differences should evolve according to the 
+// Chi-square distribution function. By inverting this last distribution, we can 
+// compute the probability that the empirical distribution matches the theoretical distribution.
+//
+// 4. The empirical histogram must match the theoretical histogram. 
+// The number of classes is NC.
+// The classes XC are computed from NC and the actual output ; X has NC+1 entries.
+// The empirical histogram is the probability that X is the intervals defined by X.
+// The theoretical histogram is computed by differences 
+// of the cumulated probability distribution function.
+// The histograms can be graphically compared with the statements:
+// plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+// plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+//
+
+//
+// In this script, we check the non-uniform random numbers.
+// We use the strategies #2 and #4.
+// The strategy #3 should be used in a future version of the tests (it requires to 
+// develop the kstest function).
+//
+
+rtol = 1.e-2;
+
+// The number of classes in the histogram
+// NC must be even.
+NC = 2*12;
+N=10000;
+
+// Set the seed to always get the same random numbers
+grand("setsd",0);
+
+/////////////////////////////////////////////////////////////
+//
+// 'exp'
+//
+
+function p = mycdfexp(X,A)
+  p = 1 - exp(-(1/A)*X);
+endfunction
+
+A = 20;
+checkMeanVariance1arg ( 400 , 800 , 'exp' , A , A , A^2 , 2*rtol );
+//
+for A = linspace(0.1,50,10)
+  checkLaw1arg ( "exp" , mycdfexp , N , NC , A , 2*rtol );
+end
+
+
+/////////////////////////////////////////////////////////////
+//
+// 'gam'
+//
+
+function p = mycdfgam (X,A,B)
+  p = cdfgam("PQ",X,A*ones(X),B*ones(X));
+endfunction
+
+A = 5;
+B = 3;
+checkMeanVariance2arg ( 400 , 800 , 'gam' , A , B , A/B , A/(B^2) , rtol );
+//
+for A = linspace(1,15,4)
+  for B = linspace(1,15,4)
+    checkLaw2arg ( 'gam' , mycdfgam , N , NC , A , B , rtol );
+  end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'bet'
+//
+
+function p = mycdfbet (X,A,B)
+  p = cdfbet("PQ",X,1-X,A*ones(X),B*ones(X));
+endfunction
+
+A=3;
+B=10;
+checkMeanVariance2arg ( 400 , 800 , 'bet' , A , B , A/(A+B) , (A*B)/((A+B)^2*(A+B+1)) , rtol );
+//
+for A = linspace(1,20,4)
+  for B = linspace(1,20,4)
+    checkLaw2arg ( 'bet' , mycdfbet , N , NC , A , B , 2*rtol );
+  end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// 'poi'
+//
+// Caution : this is a piecewise distribution
+
+function p = mycdfpoi (X,A)
+  p = cdfpoi("PQ",X,A*ones(X));
+endfunction
+
+A = 10;
+checkMeanVariance1arg ( 400 , 800 , 'poi' , A , A , A , rtol );
+//
+for A = floor(linspace(50,70,10))
+  checkPieceLaw1arg ( 'poi' , mycdfpoi , N , NC , A , rtol );
+end
+
+/////////////////////////////////////////////////////////////
+//
+// 'bin'
+//
+// Caution : this is a piecewise distribution
+
+function p = mycdfbin (X,A,B)
+  p = cdfbin("PQ",X,A*ones(X),B*ones(X),(1-B)*ones(X));
+endfunction
+
+A = 10;
+B = 0.4;
+checkMeanVariance2arg ( 400 , 800 , 'bin' , A , B , A*B , A*B*(1-B) , rtol );
+//
+for A = floor(linspace(10,50,4))
+for B = linspace(0.1,0.9,4)
+  checkPieceLaw2arg ( 'bin' , mycdfbin , N , NC , A , B , 2*rtol );
+end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// 'nor'
+//
+
+function p = mycdfnor (X,A,B)
+  p = cdfnor("PQ",X,A*ones(X),B*ones(X));
+endfunction
+
+A = 7;
+B = 12;
+checkMeanVariance2arg ( 400 , 800 , 'nor' , A , B , A , B^2 , rtol );
+//
+//
+for A = linspace(1,20,4)
+  for B = linspace(1,20,4)
+    checkLaw2arg ( 'nor' , mycdfnor , N , NC , A , B , rtol );
+  end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// 'f'
+//
+
+function p = mycdff (X,A,B)
+  p = cdff("PQ",X,A*ones(X),B*ones(X));
+endfunction
+
+//
+// Increase the tolerance for this test.
+A = 7;
+B = 12;
+checkMeanVariance2arg ( 400 , 800 , 'f' , A , B , B/(B-2) , (2*B^2*(A+B-2))/(A*(B-2)^2*(B-4)) , 4*rtol );
+//
+for A = floor(linspace(1,20,4))
+  for B = floor(linspace(5,20,4))
+    checkLaw2arg ( 'f' , mycdff , N , NC , A , B , 2*rtol );
+  end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// 'chi'
+//
+
+function p = mycdfchi (X,A)
+  p = cdfchi("PQ",X,A*ones(X));
+endfunction
+
+A = 7;
+checkMeanVariance1arg ( 400 , 800 , 'chi' , A , A , 2*A , rtol );
+
+// Increase the tolerance for this test.
+for A = floor(linspace(1,20,10))
+    checkLaw1arg ( 'chi' , mycdfchi , N , NC , A , 2*rtol );
+end
+
+/////////////////////////////////////////////////////////////
+//
+// 'nbn'
+//
+// This is a piecewise function.
+
+function p = mycdfnbn (X,A,B)
+  p = cdfnbn("PQ",X,A*ones(X),B*ones(X),(1-B)*ones(X));
+endfunction
+
+A = 7;
+B = 0.1;
+checkMeanVariance2arg ( 400 , 800 , 'nbn' , A , B , A*(1-B)/B , A*(1-B)/B^2 , rtol );
+
+// Increase the tolerance for this test.
+for A = floor(linspace(1,20,4))
+for B = linspace(0.1,0.9,4)
+    checkPieceLaw2arg ( 'nbn' , mycdfnbn , N , NC , A , B , 2*rtol );
+end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// 'nch'
+//
+
+function p = mycdfnch (X,A,B)
+  p = cdfchn("PQ",X,A*ones(X),B*ones(X));
+endfunction
+
+A = 4;
+B = 3;
+checkMeanVariance2arg ( 400 , 800 , 'nch' , A , B , A+B , 2*(A+2*B) , rtol );
+
+for A = floor(linspace(1,20,4))
+for B = floor(linspace(1,20,4))
+    checkLaw2arg ( 'nch' , mycdfnch , N , NC , A , B , 2*rtol );
+end
+end
+/////////////////////////////////////////////////////////////
+//
+// 'nf'
+// http://mathworld.wolfram.com/NoncentralF-Distribution.html
+//
+
+function p = mycdfnf (X,A,B,C)
+  p = cdffnc("PQ",X,A*ones(X),B*ones(X),C*ones(X));
+endfunction
+
+A = 4;
+B = 3;
+C = 10;
+checkMeanVariance3arg ( 400 , 800 , 'nf' , A , B , C , ((C+A)*B)/(A*(B-2)) , [] , 4*rtol );
+
+for A = floor(linspace(1,20,4))
+for B = floor(linspace(1,20,4))
+for C = floor(linspace(1,20,4))
+    checkLaw3arg ( 'nf' , mycdfnf , N , NC , A , B , C , 2*rtol );
+end
+end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// "mn"
+//
+
+n = 2^16;
+A = [1;2;3];
+// The covariance matrix B is symetric, positive definite.
+// Its diagonal entries are [4;6;5], the variances.
+B = [4,2,3;2,6,4;3,4,5];
+  R=grand(n,"mn",A,B);
+  assert_checkequal ( size(R) , [3,n] );
+  assert_checkequal ( typeof(R) , "constant" );
+  assert_checkalmostequal ( mean(R,"c") , A , 2*rtol );
+  assert_checkalmostequal ( variance(R,"c") , [4;6;5] , rtol );
+//
+// No CDF for this function => no histogram test. 
+
+/////////////////////////////////////////////////////////////
+//
+// "geom"
+// http://en.wikipedia.org/wiki/Geometric_distribution
+//
+
+A = 0.1;
+checkMeanVariance1arg ( 400 , 800 , "geom" , A , 1/A , (1-A)/A^2 , rtol );
+
+function p = mycdfgeom (X,A)
+  p = 1 -(1-A)^X
+endfunction
+//
+for A = linspace(0.1,0.9,10)
+   checkPieceLaw1arg ( "geom" , mycdfgeom , N , NC , A , 2*rtol );
+end
+
+/////////////////////////////////////////////////////////////
+//
+// "unf"
+// http://en.wikipedia.org/wiki/Uniform_distribution_(continuous)
+//
+
+A = 0.1;
+B = 2.3;
+checkMeanVariance2arg ( 400 , 800 , "unf" , A , B , (A+B)/2 , (B-A)^2/12 , rtol );
+
+function p = mycdfunf (X,A,B)
+  p = (X-A)/(B-A)
+endfunction
+//
+for A = linspace(0.1,0.9,4)
+for B = linspace(2.5,4.2,4)
+   checkLaw2arg ( "unf" , mycdfunf , N , NC , A , B , rtol );
+end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// "uin"
+// http://en.wikipedia.org/wiki/Uniform_distribution_(discrete)
+//
+// Piecewise.
+
+A = 10;
+B = 20;
+checkMeanVariance2arg ( 400 , 800 , "uin" , A , B , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
+
+function p = mycdfuin (X,A,B)
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+//
+for A = floor(linspace(10,20,4))
+for B = floor(linspace(30,40,4))
+   checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , A , B , rtol );
+end
+end
+
+/////////////////////////////////////////////////////////////
+//
+// "def"
+// http://en.wikipedia.org/wiki/Uniform_distribution_(continuous)
+//
+
+checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
+
+function p = mycdfdef (X)
+  p = X
+endfunction
+//
+checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
+
+
+
diff --git a/modules/randlib/tests/unit_tests/grand_laws2.dia.ref b/modules/randlib/tests/unit_tests/grand_laws2.dia.ref
new file mode 100755
index 000000000..cf8cfcf1c
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_laws2.dia.ref
@@ -0,0 +1,154 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+//
+// These tests makes comparisons between the empirical cumulated 
+// distribution function and the theoretical distribution function.
+// They do not make use of the Chi-square distribution function,
+// and, therefore, are not Kolmogorov-Smirnov tests.
+//
+// For each test, we can compare the two plots with the statements :
+// plot(RdevS,PS,"b-"); // Empirical distribution
+// plot(RdevS,P,"ro-"); // Theoretical distribution
+// 
+// Not all distribution functions are tested.
+// Moreover, the comparison is not done for integer-based (piecewise) 
+// distribution functions.
+// This work is to be updated.
+//
+// Set the seed to always get the same random numbers
+grand("setsd",0);
+prec = 1;
+// test for exp 
+for i=linspace(0.1,50,10)
+  N=10000;A=i;
+  Rdev=grand(1,N,'exp',A);
+  RdevS=gsort(Rdev,"g","i")';
+  PS=(1:N)'/N;
+  P=1-exp(-RdevS/A);
+  if norm(P-PS) > 2*prec then bugmes();quit;end
+end
+// test for gamma 
+for i=linspace(1,15,4)
+  for j=linspace(1,15,4)
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'gam',A,B); 
+    RdevS=gsort(Rdev,"g","i")';
+    PS=(1:N)'/N;
+    [P]=cdfgam("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+    if norm(P-PS)> 2*prec then bugmes();quit;end
+  end
+end
+// test for beta random deviate 
+for i=linspace(1,20,4)
+  for j=linspace(1,20,4)
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'bet',A,B); 
+    RdevS=gsort(Rdev,"g","i")';
+    PS=(1:N)'/N;
+    [P]=cdfbet("PQ",RdevS,1-RdevS,A*ones(RdevS),B*ones(RdevS));
+    if norm(P-PS)> 2*prec then bugmes();quit;end
+  end
+end
+// test for poi 
+for i=floor(linspace(50,70,10))
+  N=10000;A=i;
+  Rdev=grand(1,N,'poi',A);
+  RdevS=gsort(Rdev,"g","i")';
+  PS=(1:N)'/N;
+  [P]=cdfpoi("PQ",RdevS,A*ones(RdevS));
+  // Need an other test P is piecewize linear and PS
+  // linear
+  //if norm(P-PS) > prec then bugmes();quit;end
+end
+N=100;A=5;B=0.7;
+Rdev=grand(N,N,'bin',A,B); 
+Rdev=matrix(Rdev,1,N^2);
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfbin("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),(1-B)*ones(RdevS));
+//if norm(P-PS) > prec then bugmes();quit;end
+// test for f 
+N=10000;A=1;B=3;
+Rdev=grand(1,N,'f',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdff("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then bugmes();quit;end
+// test for mul
+// TODO 
+// test for nor 
+N=10000;A=1;B=2;
+Rdev=grand(1,N,'nor',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfnor("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then bugmes();quit;end
+// test for unf 
+N=10000;A=1;B=2;
+Rdev=grand(1,N,'unf',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS-A;
+if norm(P-PS) > prec then bugmes();quit;end
+// test for uin ( a finir ) 
+N=10000;A=1;B=10;
+Rdev=grand(1,N,'uin',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS-A;
+//TODO need an other test
+//if norm(P-PS) > prec then bugmes();quit;end
+// test for lgi 
+// This is a completely wrong test:
+// The output depends on the random number generator...
+N=10000;
+Rdev=grand(1,N,'lgi');
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS-A;
+//TODO need an other test 
+//if norm(P-PS) > prec then bugmes();quit;end
+// test for nbn 
+N=10000;A=5;B=0.7;
+Rdev=grand(1,N,'nbn',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfnbn("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),(1-B)*ones(RdevS));
+//TODO need an other test 
+//if norm(P-PS) > prec then bugmes();quit;end
+// test for mn 
+// TODO
+// test for 'def'
+N=10000;
+Rdev=grand(1,N,'def');
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS;
+if norm(P-PS) > prec then bugmes();quit;end
+// test for nch or chn 
+N=10000;A=5;B=4;
+Rdev=grand(1,N,'nch',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfchn("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then bugmes();quit;end
+// test for nf or fnc
+N=10000;A=5;B=4;C=10;
+Rdev=grand(1,N,'nf',A,B,C); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdffnc("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),C*ones(RdevS));
+if norm(P-PS) > prec then bugmes();quit;end
+// test for chi 
+N=10000;A=5;
+Rdev=grand(1,N,'chi',A);
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfchi("PQ",RdevS,A*ones(RdevS));
+if norm(P-PS) > prec then bugmes();quit;end
diff --git a/modules/randlib/tests/unit_tests/grand_laws2.tst b/modules/randlib/tests/unit_tests/grand_laws2.tst
new file mode 100755
index 000000000..7a3eb234b
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_laws2.tst
@@ -0,0 +1,194 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+
+//
+// These tests makes comparisons between the empirical cumulated 
+// distribution function and the theoretical distribution function.
+// They do not make use of the Chi-square distribution function,
+// and, therefore, are not Kolmogorov-Smirnov tests.
+//
+// For each test, we can compare the two plots with the statements :
+// plot(RdevS,PS,"b-"); // Empirical distribution
+// plot(RdevS,P,"ro-"); // Theoretical distribution
+// 
+// Not all distribution functions are tested.
+// Moreover, the comparison is not done for integer-based (piecewise) 
+// distribution functions.
+// This work is to be updated.
+//
+
+
+// Set the seed to always get the same random numbers
+grand("setsd",0);
+
+prec = 1;
+
+// test for exp 
+
+for i=linspace(0.1,50,10)
+  N=10000;A=i;
+  Rdev=grand(1,N,'exp',A);
+  RdevS=gsort(Rdev,"g","i")';
+  PS=(1:N)'/N;
+  P=1-exp(-RdevS/A);
+  if norm(P-PS) > 2*prec then pause,end
+end
+
+
+// test for gamma 
+for i=linspace(1,15,4)
+  for j=linspace(1,15,4)
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'gam',A,B); 
+    RdevS=gsort(Rdev,"g","i")';
+    PS=(1:N)'/N;
+    [P]=cdfgam("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+    if norm(P-PS)> 2*prec then pause,end
+  end
+end
+
+
+// test for beta random deviate 
+for i=linspace(1,20,4)
+  for j=linspace(1,20,4)
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'bet',A,B); 
+    RdevS=gsort(Rdev,"g","i")';
+    PS=(1:N)'/N;
+    [P]=cdfbet("PQ",RdevS,1-RdevS,A*ones(RdevS),B*ones(RdevS));
+    if norm(P-PS)> 2*prec then pause,end
+  end
+end
+
+// test for poi 
+for i=floor(linspace(50,70,10))
+  N=10000;A=i;
+  Rdev=grand(1,N,'poi',A);
+  RdevS=gsort(Rdev,"g","i")';
+  PS=(1:N)'/N;
+  [P]=cdfpoi("PQ",RdevS,A*ones(RdevS));
+  // Need an other test P is piecewize linear and PS
+  // linear
+  //if norm(P-PS) > prec then pause,end
+end
+
+
+N=100;A=5;B=0.7;
+Rdev=grand(N,N,'bin',A,B); 
+Rdev=matrix(Rdev,1,N^2);
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfbin("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),(1-B)*ones(RdevS));
+//if norm(P-PS) > prec then pause,end
+
+// test for f 
+
+N=10000;A=1;B=3;
+Rdev=grand(1,N,'f',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdff("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then pause,end
+
+// test for mul
+// TODO 
+
+// test for nor 
+
+N=10000;A=1;B=2;
+Rdev=grand(1,N,'nor',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfnor("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then pause,end
+
+// test for unf 
+
+N=10000;A=1;B=2;
+Rdev=grand(1,N,'unf',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS-A;
+if norm(P-PS) > prec then pause,end
+
+// test for uin ( a finir ) 
+
+N=10000;A=1;B=10;
+Rdev=grand(1,N,'uin',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS-A;
+//TODO need an other test
+//if norm(P-PS) > prec then pause,end
+
+// test for lgi 
+// This is a completely wrong test:
+// The output depends on the random number generator...
+N=10000;
+Rdev=grand(1,N,'lgi');
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS-A;
+//TODO need an other test 
+//if norm(P-PS) > prec then pause,end
+
+
+
+// test for nbn 
+
+N=10000;A=5;B=0.7;
+Rdev=grand(1,N,'nbn',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfnbn("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),(1-B)*ones(RdevS));
+//TODO need an other test 
+//if norm(P-PS) > prec then pause,end
+
+
+
+// test for mn 
+// TODO
+
+// test for 'def'
+
+N=10000;
+Rdev=grand(1,N,'def');
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=RdevS;
+if norm(P-PS) > prec then pause,end
+
+// test for nch or chn 
+
+N=10000;A=5;B=4;
+Rdev=grand(1,N,'nch',A,B); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfchn("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then pause,end
+
+// test for nf or fnc
+
+N=10000;A=5;B=4;C=10;
+Rdev=grand(1,N,'nf',A,B,C); 
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdffnc("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),C*ones(RdevS));
+if norm(P-PS) > prec then pause,end
+
+// test for chi 
+
+N=10000;A=5;
+Rdev=grand(1,N,'chi',A);
+RdevS=gsort(Rdev,"g","i")';
+PS=(1:N)'/N;
+[P]=cdfchi("PQ",RdevS,A*ones(RdevS));
+if norm(P-PS) > prec then pause,end
+
+
diff --git a/modules/randlib/tests/unit_tests/grand_plot.dia.ref b/modules/randlib/tests/unit_tests/grand_plot.dia.ref
new file mode 100755
index 000000000..5a1577b79
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_plot.dia.ref
@@ -0,0 +1,111 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- TEST WITH GRAPHIC -->
+//
+// These tests makes comparisons between the empirical cumulated 
+// distribution function and the theoretical distribution function.
+// They do not make use of the Chi-square distribution function,
+// and, therefore, are not Kolmogorov-Smirnov tests.
+// The tester is asked to visually compare the two plots, which 
+// cannot be automated.
+//
+//Comparison of pseudo-random numbers following an exponential distribution 
+//and the density of this distribution
+//Parameter of the distribution which can be modified
+lambda=1.6;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following an exponential distribution
+X = grand(1,N,"exp",lambda);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,12,25);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density 
+x=linspace(0,12,25);
+y = (1/lambda)*exp(-(1/lambda)*x);
+plot2d(x,y,3);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a beta distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=1;B=3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a beta distribution
+X = grand(1,N,"bet",A,B);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,1,50);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density 
+x=linspace(0,1,50);
+y = (1/(beta(A,B))).*(x.^(A-1)).*((1-x).^(B-1)) ;
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a gamma distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=2;B=1;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a gamma distribution
+X = grand(1,N,"gam",A,B);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,2,50);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density 
+x=linspace(0,2,50);
+y = (B/(gamma(A))).*exp(-B*x).*(B*x).^(A-1);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a binomial distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+n=50;p=0.3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a binomial distribution
+X = grand(1,N,"bin",n,p);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,n,n+1);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density
+x=linspace(0,n,n+1);
+y = binomial(p,n);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a poisson distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+mu=50;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a poisson distribution
+X = grand(1,N,"poi",mu);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,2*mu,101);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density
+[x]=linspace(0,2*mu,101);
+y = exp(-mu).*(mu.^x)./factorial(x);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
diff --git a/modules/randlib/tests/unit_tests/grand_plot.tst b/modules/randlib/tests/unit_tests/grand_plot.tst
new file mode 100755
index 000000000..3c90e48a7
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_plot.tst
@@ -0,0 +1,120 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA - Sabine Gaüzere
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- TEST WITH GRAPHIC -->
+
+//
+// These tests makes comparisons between the empirical cumulated
+// distribution function and the theoretical distribution function.
+// They do not make use of the Chi-square distribution function,
+// and, therefore, are not Kolmogorov-Smirnov tests.
+// The tester is asked to visually compare the two plots, which
+// cannot be automated.
+//
+
+//Comparison of pseudo-random numbers following an exponential distribution
+//and the density of this distribution
+//Parameter of the distribution which can be modified
+lambda=1.6;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following an exponential distribution
+X = grand(1,N,"exp",lambda);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,12,25);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density
+x=linspace(0,12,25);
+y = (1/lambda)*exp(-(1/lambda)*x);
+plot2d(x,y,3);
+f=gcf();
+delete(f);
+
+//Comparison of pseudo-random numbers following a beta distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=1;B=3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a beta distribution
+X = grand(1,N,"bet",A,B);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,1,50);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density
+x=linspace(0,1,50);
+y = (1/(beta(A,B))).*(x.^(A-1)).*((1-x).^(B-1)) ;
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+//Comparison of pseudo-random numbers following a gamma distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=2;B=1;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a gamma distribution
+X = grand(1,N,"gam",A,B);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,2,50);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density
+x=linspace(0,2,50);
+y = (B/(gamma(A))).*exp(-B*x).*(B*x).^(A-1);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+
+//Comparison of pseudo-random numbers following a binomial distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+n=50;p=0.3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a binomial distribution
+X = grand(1,N,"bin",n,p);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,n,n+1);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density
+x=linspace(0,n,n+1);
+y = binomial(p,n);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+//Comparison of pseudo-random numbers following a poisson distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+mu=50;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a poisson distribution
+X = grand(1,N,"poi",mu);
+clf();
+//Discretization of the abscisses in classes
+classes = linspace(0,2*mu,101);
+//Draw in histogram
+histplot(classes,X);
+//Draw the density
+[x]=linspace(0,2*mu,101);
+y = exp(-mu).*(mu.^x)./factorial(x);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+
diff --git a/modules/randlib/tests/unit_tests/grand_prm.dia.ref b/modules/randlib/tests/unit_tests/grand_prm.dia.ref
new file mode 100755
index 000000000..10e3e323b
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_prm.dia.ref
@@ -0,0 +1,80 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+rtol = 1.e-2;
+// Set the seed to always get the same random numbers
+grand("setsd",0);
+//
+// 1. Permute some vectors, and check that the output is basically correct.
+//
+X = (2:10)';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+//
+X = (10:100)';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+//
+X = [-12 4 9 365]';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+//
+X = (2:11)';
+P = grand(5,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , [10 5] );
+assert_checkequal ( gsort(P,"r","i") , X * ones(1,5) );
+//
+// 2. Check that the produced permutations are uniformly chosen in the 
+// set of all possible permutations.
+//
+X = [-12 4 9 365]';
+S = perms(X);
+// Total number of permutations : F = 4! = 24
+F = size(S,"r");
+// Number of random permutations : N
+N = 10000;
+P = grand(N,"prm",X);
+// Set in R(k), k=1,2,...,N the index i of the permutation, with i=1,2,...,F.
+// R must be a random variable uniformly distributed in the interval [1,2,...,F].
+R = zeros(N,1);
+for k = 1:N
+  permk = P(:,k);
+  // Search the index i of the permutation permk
+  for i = 1 : F
+    permi = S(i,:);
+    if ( permi'==permk ) then
+      R(k) = i;
+      break
+    end
+  end
+end
+assert_checkalmostequal ( mean(R) , (1+F)/2 , 0.01 );
+assert_checkalmostequal ( variance(R) , (F^2-1)/12 , 0.1 );
+//
+// Check the distribution of R
+function p = mycdfuin (X,A,B)
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+X = (1:F)';
+for k = 1 : size(X,"*")
+    EmpicicalPDF(k) = length(find(R==X(k)));
+end
+EmpicicalPDF = EmpicicalPDF./N;
+CDF = mycdfuin(X,1,F);
+TheoricPDF=[CDF(1);diff(CDF)];
+assert_checktrue( abs(EmpicicalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpicicalPDF,"bo-"); // Empirical PDF
+    plot(X,TheoricPDF,"rox-"); // Theoretical PDF
+  end
diff --git a/modules/randlib/tests/unit_tests/grand_prm.tst b/modules/randlib/tests/unit_tests/grand_prm.tst
new file mode 100755
index 000000000..0d33d454c
--- /dev/null
+++ b/modules/randlib/tests/unit_tests/grand_prm.tst
@@ -0,0 +1,85 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+
+rtol = 1.e-2;
+
+// Set the seed to always get the same random numbers
+grand("setsd",0);
+
+//
+// 1. Permute some vectors, and check that the output is basically correct.
+//
+X = (2:10)';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+//
+X = (10:100)';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+//
+X = [-12 4 9 365]';
+P = grand(1,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , size(X) );
+assert_checkequal ( gsort(P,"g","i") , X );
+//
+X = (2:11)';
+P = grand(5,"prm",X);
+assert_checkequal ( typeof(P) , "constant" );
+assert_checkequal ( size(P) , [10 5] );
+assert_checkequal ( gsort(P,"r","i") , X * ones(1,5) );
+//
+// 2. Check that the produced permutations are uniformly chosen in the 
+// set of all possible permutations.
+//
+X = [-12 4 9 365]';
+S = perms(X);
+// Total number of permutations : F = 4! = 24
+F = size(S,"r");
+// Number of random permutations : N
+N = 10000;
+P = grand(N,"prm",X);
+// Set in R(k), k=1,2,...,N the index i of the permutation, with i=1,2,...,F.
+// R must be a random variable uniformly distributed in the interval [1,2,...,F].
+R = zeros(N,1);
+for k = 1:N
+  permk = P(:,k);
+  // Search the index i of the permutation permk
+  for i = 1 : F
+    permi = S(i,:);
+    if ( permi'==permk ) then
+      R(k) = i;
+      break
+    end
+  end
+end
+assert_checkalmostequal ( mean(R) , (1+F)/2 , 0.01 );
+assert_checkalmostequal ( variance(R) , (F^2-1)/12 , 0.1 );
+//
+// Check the distribution of R
+function p = mycdfuin (X,A,B)
+  p = (floor(X)-A+1)/(B-A+1)
+endfunction
+X = (1:F)';
+for k = 1 : size(X,"*")
+    EmpicicalPDF(k) = length(find(R==X(k)));
+end
+EmpicicalPDF = EmpicicalPDF./N;
+CDF = mycdfuin(X,1,F);
+TheoricPDF=[CDF(1);diff(CDF)];
+assert_checktrue( abs(EmpicicalPDF-TheoricPDF) < rtol );
+  if ( %f ) then
+    plot(X,EmpicicalPDF,"bo-"); // Empirical PDF
+    plot(X,TheoricPDF,"rox-"); // Theoretical PDF
+  end
+
+