* Bug #14057 fixed - grand(m,n) returned a wrong error and grand(m,n,p) called %s_gra...
[scilab.git] / scilab / modules / randlib / sci_gateway / cpp / sci_grand.cpp
index 08d9fe9..a999b0b 100644 (file)
@@ -1,6 +1,7 @@
 /*
 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
+* Copyright (C) 2014 - Scilab Enterprises - Sylvain GENIN
 *
 * This file must be used under the terms of the CeCILL.
 * This source file is licensed as described in the file COPYING, which
 #include "double.hxx"
 #include "string.hxx"
 #include "configvariable.hxx"
+#include "int.hxx"
+#include "polynom.hxx"
+#include "sparse.hxx"
+#include "overload.hxx"
+#include "execvisitor.hxx"
 
 extern "C"
 {
-#include "MALLOC.h"
+#include "sci_malloc.h"
 #include "localization.h"
 #include "Scierror.h"
 #include "sciprint.h"
@@ -28,35 +34,40 @@ extern "C"
 }
 /*--------------------------------------------------------------------------*/
 
+template<class U>
+void sci_grand_prm(int iNumIter, U *pIn, types::InternalType** pOut);
+
 types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, types::typed_list &out)
 {
-    enum {
+    enum
+    {
         MT, KISS, CLCG4, CLCG2, URAND, FSULTRA
     };
 
     //  names at the scilab level
-    wchar_t* names_gen[6] = {L"mt", L"kiss", L"clcg4", L"clcg2", L"urand", L"fsultra"};
+    const wchar_t* names_gen[6] = { L"mt", L"kiss", L"clcg4", L"clcg2", L"urand", L"fsultra" };
 
-    types::String* pStrMethod   = NULL;
+    types::String* pStrMethod = NULL;
     types::String* pStrGenOrPhr = NULL;
 
     std::vector<types::Double*> vectpDblInput;
 
-    int iStrPos     = 0;
-    int iPos        = 0;
-    int meth        = -1;
-    int iRows       = -1;
-    int iCols       = -1;
-    int iNumIter    = -1;
+    int iStrPos = 0;
+    int iPos = 0;
+    int meth = -1;
+    int iRows = -1;
+    int iCols = -1;
+    int iNumIter = -1;
+
 
     int current_base_gen = ConfigVariable::getCurrentBaseGen();
 
     // *** check the maximal number of input args. ***
-    if (in.size() > 6)
+    /* if (in.size() > 6)
     {
-        Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "grand", 1, 6);
-        return types::Function::Error;
-    }
+    Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "grand", 1, 6);
+    return types::Function::Error;
+    }*/
 
     // *** check number of output args. ***
     if (_iRetCount > 1)
@@ -65,7 +76,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         return types::Function::Error;
     }
 
-    // *** find the mothod string. ***
+    // *** find the method string. ***
     for (int i = 0; i < in.size(); i++)
     {
         if (in[i]->isString())
@@ -78,10 +89,27 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
     if (pStrMethod == NULL)
     {
-        Scierror(78, _("%s: Wrong number of output argument(s): At least %d string expected.\n"), "grand", 1);
+        for (int i = 0; i < in.size(); i++)
+        {
+            if (in[i]->isDouble() == false)
+            {
+                ast::ExecVisitor exec;
+                std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_grand";
+                return Overload::call(wstFuncName, in, _iRetCount, out, &exec);
+            }
+        }
+
+        Scierror(999, _("%s: Wrong type for input argument #%d: A string expected.\n"), "grand", in.size() + 1);
         return types::Function::Error;
     }
 
+    int iDims = iStrPos > 1 ? iStrPos : 2;
+    int *itab = new int[iStrPos > 1 ? iStrPos : 2];
+
+    //all variables are matrix
+    itab[0] = 1;
+    itab[1] = 1;
+
     wchar_t* wcsMeth = pStrMethod->get(0);
     int iNumInputArg = 5;
     if (wcscmp(wcsMeth, L"bet") == 0) // beta
@@ -219,7 +247,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                     return types::Function::Error;
                 }
 
-                iNumInputArg = in.size();
+                iNumInputArg = (int)in.size();
                 break;
             }
         }
@@ -286,16 +314,17 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             return types::Function::Error;
         }
 
-        iNumIter = pDblTemp->get(0);
+        iNumIter = (int)pDblTemp->get(0);
         iPos++;
     }
     else if (meth < 21) // grand(m, n, "...", ... || grand(matrix, "...", ...
     {
-        if (iStrPos != 1 && iStrPos != 2)
+
+        /*if (iStrPos != 1 && iStrPos != 2)
         {
-            Scierror(999, _("%s: Wrong position for input argument #%d : Must be in position %d or %d.\n"), "grand", iStrPos + 1, 2, 3);
-            return types::Function::Error;
-        }
+        Scierror(999, _("%s: Wrong position for input argument #%d : Must be in position %d or %d.\n"), "grand", iStrPos + 1, 2, 3);
+        return types::Function::Error;
+        }*/
 
         std::vector<types::Double*> vectpDblTemp;
         for (int i = 0; i < iStrPos; i++)
@@ -316,16 +345,10 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             iPos++;
         }
 
-        if (iStrPos == 2)
-        {
-            iRows = static_cast<int>(vectpDblTemp[0]->get(0));
-            iCols = static_cast<int>(vectpDblTemp[1]->get(0));;
-        }
-        else // iStrPos == 2
+        //get number of dimensions to output
+        for (int i = 0; i < iStrPos; i++)
         {
-            iRows = vectpDblTemp[0]->getRows();
-            iCols = vectpDblTemp[0]->getCols();
-            iNumInputArg--;
+            itab[i] = static_cast<int>(vectpDblTemp[i]->get(0));
         }
     }
 
@@ -336,15 +359,6 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
     }
     iPos++; // method string has been already got.
 
-    // *** check the number of input args according the methode. ***
-    if (in.size() != iNumInputArg)
-    {
-        char* pstMeth = wide_string_to_UTF8(wcsMeth);
-        Scierror(77, _("%s: Wrong number of input argument(s) for method %s: %d expected.\n"), "grand", pstMeth, iNumInputArg);
-        FREE(pstMeth);
-        return types::Function::Error;
-    }
-
     // *** get arguments after methode string. ***
     if (meth == 22 || meth == 25) // setgen || phr2sd
     {
@@ -362,8 +376,17 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             return types::Function::Error;
         }
     }
+    else if (meth == 16)
+    {
+        if (in.size() > 3)
+        {
+            Scierror(999, _("Missing vect for random permutation\n"));
+            return types::Function::Error;
+        }
+    }
     else
     {
+
         for (int i = iPos; i < in.size(); i++)
         {
             if (in[i]->isDouble() == false)
@@ -378,28 +401,73 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
     // *** perform operation in according method string and return result. ***
 
-    if (iCols * iRows == 0)
+    types::Double* pDblOut = NULL;
+
+    if (iStrPos >= 2)
+    {
+        int iProdiTab = 1;
+        for (int i = 0; i < iStrPos; i++)
+        {
+            iProdiTab = iProdiTab * itab[i];
+        }
+
+        if (iProdiTab == 0)
+        {
+            pDblOut = types::Double::Empty();
+        }
+        else if ((itab[1] != itab[0]) && (itab[0] <= -1))
+        {
+            Scierror(999, _("%s: Wrong value for input argument #%d: Positive scalar expected.\n"), "grand", 1);
+            return types::Function::Error;
+        }
+        else if ((itab[1] != itab[0]) && (itab[1] <= -1))
+        {
+            Scierror(999, _("%s: Wrong value for input argument #%d: Positive scalar expected.\n"), "grand", 2);
+            return types::Function::Error;
+        }
+        else
+        {
+            pDblOut = new types::Double(iDims, itab);
+        }
+    }
+    else
     {
-        out.push_back(types::Double::Empty());
+
+        types::Double* pDblIn = in[0]->getAs<types::Double>();
+        if (meth == 14)//'mul'
+        {
+            int* iDimsArraytempo = new int[2];
+            iDimsArraytempo[0] = in[3]->getAs<types::Double>()->getSize() + 1;
+            iDimsArraytempo[1] = iNumIter;
+            pDblOut = new types::Double(pDblIn->getDims(), iDimsArraytempo);
+        }
+        else
+        {
+            pDblOut = new types::Double(pDblIn->getDims(), pDblIn->getDimsArray());
+        }
+
     }
 
+    delete[] itab;
+
     switch (meth)
     {
         case 1: // beta
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
             double minlog   = 1.e-37;
 
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
 
                 if (vectpDblInput[i]->get(0) < minlog)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong value for input argument #%d : At least %lf expected.\n"), "grand", iPos + 1, minlog);
                     return types::Function::Error;
                 }
@@ -416,12 +484,11 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         case 2: // binomial
         case 3: // negative binomial
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -429,12 +496,14 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (vectpDblInput[0]->get(0) < 0.0) // N
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive integer expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[1]->get(0) < 0.0 || vectpDblInput[1]->get(0) > 1.0) // p
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : A value  expected.\n"), "grand", 5);
                 return types::Function::Error;
             }
@@ -460,16 +529,16 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 4: // chisquare
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             if (vectpDblInput[0]->isScalar() == false)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[0]->get(0) <= 0.0) // Df
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive no null value expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
@@ -484,12 +553,11 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 5: // non central chisquare
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -497,12 +565,14 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (vectpDblInput[0]->get(0) < 1.0) // Df
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : value >1.0 expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[1]->get(0) < 0.0) // Xnon
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 5);
                 return types::Function::Error;
             }
@@ -517,16 +587,16 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 6: // exponential
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             if (vectpDblInput[0]->isScalar() == false)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[0]->get(0) < 0.0) // Av
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
@@ -541,18 +611,18 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 7: // F variance ratio
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
 
                 if (vectpDblInput[i]->get(0) <= 0.0) // Dfn Dfd
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong value for input argument #%d : Positive no null value expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -568,12 +638,11 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 8: // non F variance ratio
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 3; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -581,18 +650,21 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (vectpDblInput[0]->get(0) < 1.0) // Dfn
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : value > 1.0 expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[1]->get(0) <= 0.0) // Dfd
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive non null value expected.\n"), "grand", 5);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[2]->get(0) < 0.0) // Xnon
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 6);
                 return types::Function::Error;
             }
@@ -607,18 +679,18 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 9: // gamma
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
 
                 if (vectpDblInput[i]->get(0) <= 0.0)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong value for input argument #%d : Positive non null value expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -637,12 +709,11 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 10: // Gauss Laplace (normal)
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -650,6 +721,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (vectpDblInput[1]->get(0) < 0.0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 5);
                 return types::Function::Error;
             }
@@ -666,18 +738,21 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         {
             if (vectpDblInput[0]->getCols() != 1 || vectpDblInput[0]->getSize() == 0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong size for input argument #%d : A matrix of size m x 1 expected.(m > 0)\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[0]->getRows() != vectpDblInput[1]->getRows())
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong size for input argument #%d and #%d: Mean and Cov have incompatible dimensions.\n"), "grand", 4, 5);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[1]->getRows() != vectpDblInput[1]->getCols())
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong size for input argument #%d : A square symmetric positive definite matrix expected.\n"), "grand", 5);
                 return types::Function::Error;
             }
@@ -686,9 +761,10 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             int size = vectpDblInput[0]->getRows();
             int mp   = (size * (size + 3)) / 2 + 1;
 
+            delete pDblOut;
             types::Double* pDblOut = new types::Double(size, iNumIter);
-            double* work  = (double*)malloc(size * sizeof(double));
-            double* param = (double*)malloc(mp   * sizeof(double));
+            double* work = new double[size];
+            double* param = new double[mp];
 
             types::Double* pDblMean = vectpDblInput[0]->clone()->getAs<types::Double>();
             types::Double* pDblCov  = vectpDblInput[1]->clone()->getAs<types::Double>();
@@ -700,8 +776,9 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (ierr == 1)
             {
-                free(work);
-                free(param);
+                delete work;
+                delete param;
+                delete pDblOut;
                 Scierror(999, _("%s: setgmn return with state %d.\n"), "grand", ierr);
                 return types::Function::Error;
             }
@@ -711,25 +788,26 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                 C2F(genmn)(param, pDblOut->get() + (size * i), work);
             }
 
-            free(work);
-            free(param);
+            delete work;
+            delete param;
 
             out.push_back(pDblOut);
             break;
         }
         case 12: // geometric
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
             double pmin = 1.3e-307;
 
             if (vectpDblInput[0]->isScalar() == false)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[0]->get(0) < pmin || vectpDblInput[0]->get(0) > 1.0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Must be between %lf and %d.\n"), "grand", 4, pmin, 1);
                 return types::Function::Error;
             }
@@ -747,12 +825,14 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             if ( vectpDblInput[0]->getRows() != vectpDblInput[0]->getCols() &&
                     vectpDblInput[0]->getRows() != 1)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong size for input argument #%d : A square matrix or a row vector expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[1]->getSize() == 0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong size for input argument #%d: No empty matrix expected.\n"), "grand", 5);
                 return types::Function::Error;
             }
@@ -768,7 +848,6 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                 }
             }
 
-            types::Double* pDblOut = new types::Double(sizeOfX0, iNumIter);
             int size = vectpDblInput[0]->getSize() + vectpDblInput[0]->getCols();
             double* work = (double*)malloc(size * sizeof(double));
 
@@ -781,6 +860,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
                     if (vectpDblInput[0]->get(position) < 0 || vectpDblInput[0]->get(position) > 1)
                     {
+                        delete pDblOut;
                         Scierror(999, _("%s: Wrong value for input argument #%d: P(%d,%d) must be in the range [0 1].\n"), "grand", i + 1, j + 1);
                         return types::Function::Error;
                     }
@@ -790,6 +870,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
                 if (fabs(ptot - 1.0) > 1e-8)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Sum of P(%d,1:%d)=%lf ~= 1.\n"), "grand", i + 1, vectpDblInput[0]->getCols(), ptot);
                     return types::Function::Error;
                 }
@@ -827,7 +908,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                     }
 
                     // ** projection to avoid boundaries **
-                    niv = Max(Min(niv, vectpDblInput[0]->getCols()), 1);
+                    niv = std::max(std::min(niv, vectpDblInput[0]->getCols()), 1);
                     pDblOut->set(vectpDblInput[1]->getSize() * j + i, static_cast<double>(niv));
                     iCur = niv - 1;
                 }
@@ -840,30 +921,33 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         {
             if (vectpDblInput[0]->isScalar() == false)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[0]->get(0) < 0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d: A positive scalar expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[1]->getCols() != 1 || vectpDblInput[1]->getRows() <= 0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong size for input argument #%d: A colomn vector expected.\n"), "grand", 5);
                 return types::Function::Error;
             }
 
             double ptot = 0.0;
             int ncat = vectpDblInput[1]->getRows() + 1;
-            types::Double* pDblOut = new types::Double(ncat, iNumIter);
 
             for (int i = 0; i < vectpDblInput[1]->getCols(); i++)
             {
                 if (vectpDblInput[1]->get(i) < 0.0 || vectpDblInput[1]->get(i) > 1.0)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong value for input argument #%d: P(%d) must be in the range [0 1].\n"), "grand", 5, i + 1);
                     return types::Function::Error;
                 }
@@ -872,12 +956,13 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (ptot > 1.0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d: Sum of P(i) > 1.\n"), "grand", 5);
                 return types::Function::Error;
             }
 
-            int* piP    = (int*)malloc(vectpDblInput[0]->getSize() * sizeof(int));
-            int* piOut  = (int*)malloc(pDblOut->getSize() * sizeof(int));
+            int* piP = new int[vectpDblInput[0]->getSize()];
+            int* piOut = new int[pDblOut->getSize()];
 
             for (int i = 0; i < vectpDblInput[0]->getSize(); i++)
             {
@@ -894,24 +979,24 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                 pDblOut->set(i, static_cast<double>(piOut[i]));
             }
 
-            free(piP);
-            free(piOut);
+            delete piP;
+            delete piOut;
 
             out.push_back(pDblOut);
             break;
         }
         case 15: // Poisson
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             if (vectpDblInput[0]->isScalar() == false)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
 
             if (vectpDblInput[0]->get(0) < 0.0)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 4);
                 return types::Function::Error;
             }
@@ -926,31 +1011,133 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 16: // random permutations
         {
-            if (vectpDblInput[0]->getCols() != 1)
+            delete pDblOut;
+            types::InternalType* pITOut = NULL;
+
+            // get dims and isComplex
+            int iDims = 0;
+            int* piDimsArray = NULL;
+            bool bIsComplex = false;
+            if (in[2]->isGenericType())
             {
-                Scierror(999, _("%s: Wrong type for input argument #%d : A colomn vector expected.\n"), "grand", 4);
-                return types::Function::Error;
+                types::GenericType* pGT = in[2]->getAs<types::GenericType>();
+                iDims = pGT->getDims();
+                piDimsArray = pGT->getDimsArray();
+                bIsComplex = pGT->isComplex();
             }
 
-            int iVectRows = vectpDblInput[0]->getRows();
-            types::Double* pDblOut = new types::Double(iVectRows, iNumIter);
-
-            for (int i = 0; i < iNumIter; i++)
+            switch (in[2]->getType())
             {
-                for (int j = 0; j < iVectRows; j++)
+                case types::InternalType::ScilabInt8:
+                    pITOut = new types::Int8(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::Int8>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabUInt8:
+                    pITOut = new types::UInt8(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::UInt8>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabInt16:
+                    pITOut = new types::Int16(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::Int16>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabUInt16:
+                    pITOut = new types::UInt16(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::UInt16>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabInt32:
+                    pITOut = new types::Int32(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::Int32>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabUInt32:
+                    pITOut = new types::UInt32(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::UInt32>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabInt64:
+                    pITOut = new types::Int64(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::Int64>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabUInt64:
+                    pITOut = new types::UInt64(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::UInt64>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabDouble:
+                    pITOut = new types::Double(iDims, piDimsArray, bIsComplex);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::Double>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabBool:
+                    pITOut = new types::Bool(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::Bool>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabString:
+                    pITOut = new types::String(iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::String>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabPolynom:
+                    pITOut = new types::Polynom(in[2]->getAs<types::Polynom>()->getVariableName(), iDims, piDimsArray);
+                    sci_grand_prm(iNumIter, in[2]->getAs<types::Polynom>(), &pITOut);
+                    break;
+                case types::InternalType::ScilabSparse:
                 {
-                    pDblOut->set(iVectRows * i + j, vectpDblInput[0]->get(j));
+                    std::complex<double> cplxDbl;
+                    types::InternalType* pITOutTempo = NULL;
+                    types::Double* pDblTempo = NULL;
+                    types::Sparse* pSP = in[2]->getAs<types::Sparse>();
+                    int isize = pSP->getSize();
+                    pITOut = new types::Sparse(pSP->getRows(), pSP->getCols(), pSP->isComplex());
+                    pDblTempo = new types::Double(isize, 1, pSP->isComplex());
+                    pITOutTempo = new types::Double(isize, iNumIter, pSP->isComplex());
+
+                    if (pDblTempo->isComplex())
+                    {
+                        for (int i = 0; i < isize; i++)
+                        {
+                            cplxDbl = in[2]->getAs<types::Sparse>()->getImg(i);
+                            pDblTempo->set(i, cplxDbl.real());
+                            pDblTempo->set(i, cplxDbl.imag());
+                        }
+                    }
+                    else
+                    {
+                        for (int i = 0; i < isize; i++)
+                        {
+                            pDblTempo->set(i, in[2]->getAs<types::Sparse>()->get(i));
+                        }
+                    }
+                    sci_grand_prm(iNumIter, pDblTempo, &pITOutTempo);
+
+                    if (pDblTempo->isComplex())
+                    {
+                        for (int i = 0; i < isize; i++)
+                        {
+                            cplxDbl.real(pITOutTempo->getAs<types::Double>()->get(i));
+                            cplxDbl.imag(pITOutTempo->getAs<types::Double>()->getImg(i));
+                            pITOutTempo->getAs<types::Sparse>()->set(i, cplxDbl);
+                        }
+                    }
+                    else
+                    {
+                        for (int i = 0; i < isize; i++)
+                        {
+                            pITOut->getAs<types::Sparse>()->set(i, pITOutTempo->getAs<types::Double>()->get(i));
+                        }
+                    }
+
+                    delete pITOutTempo;
+                    delete pDblTempo;
+                    break;
+                }
+                default:
+                {
+                    Scierror(999, _("%s: Wrong type for input argument: Matrix (full or sparse) or Hypermatrix of Reals, Complexes, Integers, Booleans, Strings or Polynomials expected.\n"), "grand");
+                    return types::Function::Error;
                 }
-                C2F(genprm)(pDblOut->get() + (iVectRows * i), &iVectRows);
             }
 
-            out.push_back(pDblOut);
+            out.push_back(pITOut);
             break;
         }
         case 17: // uniform (def)
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < pDblOut->getSize(); i++)
             {
                 pDblOut->set(i, C2F(ranf)());
@@ -961,12 +1148,11 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 18: // uniform (unf)
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -977,6 +1163,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (low > high)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low < High expected.\n"), "grand", 4, 5);
                 return types::Function::Error;
             }
@@ -991,12 +1178,11 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 19: // uniform (uin)
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < 2; i++)
             {
                 if (vectpDblInput[i]->isScalar() == false)
                 {
+                    delete pDblOut;
                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
                     return types::Function::Error;
                 }
@@ -1007,6 +1193,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 
             if (low > high)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low < High expected.\n"), "grand", 4, 5);
                 return types::Function::Error;
             }
@@ -1015,6 +1202,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                     high != vectpDblInput[1]->get(0) ||
                     (high - low + 1) > 2147483561)
             {
+                delete pDblOut;
                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low and High must be integers and (high - low + 1) <=  2147483561.\n"), "grand", 4, 5);
                 return types::Function::Error;
             }
@@ -1029,8 +1217,6 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 20: // uniform (lgi)
         {
-            types::Double* pDblOut = new types::Double(iRows, iCols);
-
             for (int i = 0; i < pDblOut->getSize(); i++)
             {
                 pDblOut->set(i, static_cast<double>(ignlgi()));
@@ -1041,12 +1227,14 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 21: // getgen
         {
+            delete pDblOut;
             types::String* pStrOut = new types::String(names_gen[current_base_gen]);
             out.push_back(pStrOut);
             break;
         }
         case 22: // setgen
         {
+            delete pDblOut;
             wchar_t* wcsGen = pStrGenOrPhr->get(0);
 
             if (wcscmp(wcsGen, L"mt") == 0)
@@ -1086,8 +1274,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 23: // getsd
         {
-            types::Double* pDblOut = NULL;
-
+            delete pDblOut;
             switch (current_base_gen)
             {
                 case MT:
@@ -1134,6 +1321,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 24: // setsd
         {
+            delete pDblOut;
             int ierr = 0;
             switch (current_base_gen)
             {
@@ -1143,7 +1331,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                     {
                         ierr = set_state_mt_simple(vectpDblInput[0]->get(0));
                     }
-                    else if (vectpDblInput[0]->getSize() != 625)
+                    else if (vectpDblInput[0]->getSize() == 625)
                     {
                         ierr = set_state_mt(vectpDblInput[0]->get());
                     }
@@ -1250,6 +1438,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 25: // phr2sd
         {
+            delete pDblOut;
             if (pStrGenOrPhr->isScalar() == false)
             {
                 Scierror(999, _("%s: Wrong type for input argument #%d : One string expected.\n"), "grand", 2);
@@ -1257,7 +1446,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             }
 
             types::Double* pDblOut = new types::Double(1, 2);
-            int size = wcslen(pStrGenOrPhr->get(0));
+            int size = (int)wcslen(pStrGenOrPhr->get(0));
             int piOut[2];
             char* strPhr = wide_string_to_UTF8(pStrGenOrPhr->get(0));
 
@@ -1272,6 +1461,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 26: // setcgn
         {
+            delete pDblOut;
             if (current_base_gen != CLCG4)
             {
                 sciprint(_("The %s option affects only the %s generator\n"), "setcgn", "clcg4");
@@ -1296,12 +1486,14 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 27: // getcgn
         {
+            delete pDblOut;
             double dOut = static_cast<double>(ConfigVariable::getCurrentClcg4());
             out.push_back(new types::Double(dOut));
             break;
         }
         case 28: // initgn
         {
+            delete pDblOut;
             SeedType where;
             if (current_base_gen != CLCG4)
             {
@@ -1329,6 +1521,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 29: // setall
         {
+            delete pDblOut;
             if (current_base_gen != CLCG4)
             {
                 sciprint(_("The %s option affects only the %s generator\n"), "setall", "clcg4");
@@ -1358,6 +1551,7 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
         }
         case 30: // advnst
         {
+            delete pDblOut;
             if (current_base_gen != CLCG4)
             {
                 sciprint(_("The %s option affects only the %s generator\n"), "advnst", "clcg4");
@@ -1387,3 +1581,76 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 }
 /*--------------------------------------------------------------------------*/
 
+template<class U>
+void sci_grand_prm(int iNumIter, U *pIn, types::InternalType** pOut)
+{
+    U* pUTempo = NULL;
+    types::InternalType* pITTempo = NULL;
+    int* piDimsArray = NULL;
+    int Dims = 0;
+
+    if ((pIn->getCols() == 1) && (pIn->getDims() == 2))
+    {
+        pOut[0]->getAs<U>()->resize(pIn->getRows(), iNumIter);
+        pUTempo = pIn;
+    }
+    else if ((pIn->getRows() == 1) && (pIn->getDims() == 2))
+    {
+        pIn->transpose(pITTempo);
+        pOut[0]->getAs<U>()->resize(iNumIter, pIn->getCols());
+        pUTempo = pITTempo->getAs<U>();
+    }
+    else
+    {
+        piDimsArray = pOut[0]->getAs<U>()->getDimsArray();
+        Dims = pOut[0]->getAs<U>()->getDims();
+        piDimsArray[Dims] = iNumIter;
+        pOut[0]->getAs<U>()->resize(piDimsArray, Dims + 1);
+        pUTempo = pIn;
+    }
+
+    int isize = pUTempo->getSize();
+
+    types::Double* pDblOut = new types::Double(isize, iNumIter, pUTempo->isComplex());
+
+    for (int i = 0; i < iNumIter; i++)
+    {
+        for (int j = 0; j < isize; j++)
+        {
+            pDblOut->set(isize * i + j, j);
+        }
+        C2F(genprm)(pDblOut->get() + (isize * i), &isize);
+    }
+
+    if ((pIn->getCols() != 1) && (pIn->getRows() == 1) && (pIn->getDims() == 2))
+    {
+        pDblOut->transpose(pITTempo);
+        delete pDblOut;
+        pDblOut = pITTempo->getAs<types::Double>();
+    }
+
+    if (pUTempo->isComplex() && pUTempo->isPoly() == false)
+    {
+        for (int i = 0; i < pOut[0]->getAs<U>()->getSize(); i++)
+        {
+            pOut[0]->getAs<U>()->set(i , pIn->get(static_cast<int>(pDblOut->get(i))));
+            pOut[0]->getAs<U>()->setImg(i, pIn->getImg(static_cast<int>(pDblOut->get(i))));
+        }
+    }
+    else
+    {
+        for (int i = 0; i < pOut[0]->getAs<U>()->getSize(); i++)
+        {
+            pOut[0]->getAs<U>()->set(i, pIn->get(static_cast<int>(pDblOut->get(i))));
+        }
+    }
+
+    if ((pIn->getCols() != 1) && (pIn->getRows() == 1) && (pIn->getDims() == 2))
+    {
+        delete pUTempo;
+    }
+
+    delete pDblOut;
+}
+
+/*--------------------------------------------------------------------------*/