fix spgrand for :
[scilab.git] / scilab / modules / randlib / sci_gateway / cpp / sci_grand.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Cedric DELAMARRE
4 * Copyright (C) 2014 - Scilab Enterprises - Sylvain GENIN
5 *
6 * This file must be used under the terms of the CeCILL.
7 * This source file is licensed as described in the file COPYING, which
8 * you should have received as part of this distribution.  The terms
9 * are also available at
10 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 *
12 */
13 /*--------------------------------------------------------------------------*/
14 #include "randlib_gw.hxx"
15 #include "function.hxx"
16 #include "double.hxx"
17 #include "string.hxx"
18 #include "configvariable.hxx"
19 #include "int.hxx"
20 #include "polynom.hxx"
21 #include "sparse.hxx"
22 #include "overload.hxx"
23 #include "execvisitor.hxx"
24
25 extern "C"
26 {
27 #include "sci_malloc.h"
28 #include "localization.h"
29 #include "Scierror.h"
30 #include "sciprint.h"
31 #include "grand.h"
32 #include "clcg4.h"
33 #include "others_generators.h"
34 }
35 /*--------------------------------------------------------------------------*/
36
37 template<class U>
38 void sci_grand_prm(int iNumIter, U *pIn, types::InternalType** pOut);
39
40 types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, types::typed_list &out)
41 {
42     enum
43     {
44         MT, KISS, CLCG4, CLCG2, URAND, FSULTRA
45     };
46
47     //  names at the scilab level
48     const wchar_t* names_gen[6] = { L"mt", L"kiss", L"clcg4", L"clcg2", L"urand", L"fsultra" };
49
50     types::String* pStrMethod = NULL;
51     types::String* pStrGenOrPhr = NULL;
52
53     std::vector<types::Double*> vectpDblInput;
54
55     int iStrPos = 0;
56     int iPos = 0;
57     int meth = -1;
58     int iRows = -1;
59     int iCols = -1;
60     int iNumIter = -1;
61
62
63     int current_base_gen = ConfigVariable::getCurrentBaseGen();
64
65     // *** check the maximal number of input args. ***
66     /* if (in.size() > 6)
67     {
68     Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "grand", 1, 6);
69     return types::Function::Error;
70     }*/
71
72     // *** check number of output args. ***
73     if (_iRetCount > 1)
74     {
75         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), "grand", 1);
76         return types::Function::Error;
77     }
78
79     // *** find the mothod string. ***
80     for (int i = 0; i < in.size(); i++)
81     {
82         if (in[i]->isString())
83         {
84             pStrMethod = in[i]->getAs<types::String>();
85             iStrPos = i;
86
87             break;
88         }
89     }
90
91     if ((iStrPos == 0) && (in[0]->isString() == false))
92     {
93         std::wstring wstFuncName = L"%" + in[0]->getShortTypeStr() + L"_grand";
94         return Overload::call(wstFuncName, in, _iRetCount, out, new ast::ExecVisitor());
95     }
96
97     int iDims = iStrPos > 1 ? iStrPos : 2;
98     int *itab = new int[iStrPos > 1 ? iStrPos : 2];
99
100     //all variables are matrix
101     itab[0] = 1;
102     itab[1] = 1;
103
104     if (pStrMethod == NULL)
105     {
106         Scierror(78, _("%s: Wrong number of output argument(s): At least %d string expected.\n"), "grand", 1);
107         return types::Function::Error;
108     }
109
110     wchar_t* wcsMeth = pStrMethod->get(0);
111     int iNumInputArg = 5;
112     if (wcscmp(wcsMeth, L"bet") == 0) // beta
113     {
114         meth = 1;
115     }
116     else if (wcscmp(wcsMeth, L"bin") == 0) // binomial
117     {
118         meth = 2;
119     }
120     else if (wcscmp(wcsMeth, L"nbn") == 0) // negative binomial
121     {
122         meth = 3;
123     }
124     else if (wcscmp(wcsMeth, L"chi") == 0) // chisquare
125     {
126         iNumInputArg = 4;
127         meth = 4;
128     }
129     else if (wcscmp(wcsMeth, L"nch") == 0) // non central chisquare
130     {
131         meth = 5;
132     }
133     else if (wcscmp(wcsMeth, L"exp") == 0) // exponential
134     {
135         iNumInputArg = 4;
136         meth = 6;
137     }
138     else if (wcscmp(wcsMeth, L"f") == 0) // F variance ratio
139     {
140         meth = 7;
141     }
142     else if (wcscmp(wcsMeth, L"nf") == 0) // non central F variance ratio
143     {
144         iNumInputArg = 6;
145         meth = 8;
146     }
147     else if (wcscmp(wcsMeth, L"gam") == 0) // gamma
148     {
149         meth = 9;
150     }
151     else if (wcscmp(wcsMeth, L"nor") == 0) // Gauss Laplace (normal)
152     {
153         meth = 10;
154     }
155     else if (wcscmp(wcsMeth, L"mn") == 0) // multivariate gaussian (multivariate normal)
156     {
157         iNumInputArg = 4;
158         meth = 11;
159     }
160     else if (wcscmp(wcsMeth, L"geom") == 0) // geometric
161     {
162         iNumInputArg = 4;
163         meth = 12;
164     }
165     else if (wcscmp(wcsMeth, L"markov") == 0) // markov
166     {
167         iNumInputArg = 4;
168         meth = 13;
169     }
170     else if (wcscmp(wcsMeth, L"mul") == 0) // multinomial
171     {
172         iNumInputArg = 4;
173         meth = 14;
174     }
175     else if (wcscmp(wcsMeth, L"poi") == 0) // Poisson
176     {
177         iNumInputArg = 4;
178         meth = 15;
179     }
180     else if (wcscmp(wcsMeth, L"prm") == 0) // random permutations
181     {
182         iNumInputArg = 3;
183         meth = 16;
184     }
185     else if (wcscmp(wcsMeth, L"def") == 0) // uniform (def)
186     {
187         iNumInputArg = 3;
188         meth = 17;
189     }
190     else if (wcscmp(wcsMeth, L"unf") == 0) // uniform (unf)
191     {
192         meth = 18;
193     }
194     else if (wcscmp(wcsMeth, L"uin") == 0) // uniform (uin)
195     {
196         meth = 19;
197     }
198     else if (wcscmp(wcsMeth, L"lgi") == 0) // uniform (lgi)
199     {
200         iNumInputArg = 3;
201         meth = 20;
202     }
203     else if (wcscmp(wcsMeth, L"getgen") == 0) // getgen
204     {
205         iNumInputArg = 1;
206         meth = 21;
207     }
208     else if (wcscmp(wcsMeth, L"setgen") == 0) // setgen
209     {
210         iNumInputArg = 2;
211         meth = 22;
212     }
213     else if (wcscmp(wcsMeth, L"getsd") == 0) // getsd
214     {
215         iNumInputArg = 1;
216         meth = 23;
217     }
218     else if (wcscmp(wcsMeth, L"setsd") == 0) // setsd
219     {
220         switch (current_base_gen)
221         {
222             case MT:
223                 iNumInputArg = 2;
224                 break;
225             case KISS:
226                 iNumInputArg = 5;
227                 break;
228             case CLCG2:
229                 iNumInputArg = 3;
230                 break;
231             case CLCG4:
232                 iNumInputArg = 5;
233                 break;
234             case URAND:
235                 iNumInputArg = 2;
236                 break;
237             case FSULTRA:
238             {
239                 if (in.size() != 2 && in.size() != 3)
240                 {
241                     char* pstMeth = wide_string_to_UTF8(wcsMeth);
242                     Scierror(77, _("%s: Wrong number of input argument(s) for method %s: %d or %d expected.\n"), "grand", pstMeth, 2, 3);
243                     FREE(pstMeth);
244                     return types::Function::Error;
245                 }
246
247                 iNumInputArg = (int)in.size();
248                 break;
249             }
250         }
251         meth = 24;
252     }
253     else if (wcscmp(wcsMeth, L"phr2sd") == 0) // phr2sd
254     {
255         iNumInputArg = 2;
256         meth = 25;
257     }
258     else if (wcscmp(wcsMeth, L"setcgn") == 0) // setcgn
259     {
260         iNumInputArg = 2;
261         meth = 26;
262     }
263     else if (wcscmp(wcsMeth, L"getcgn") == 0) // getcgn
264     {
265         iNumInputArg = 1;
266         meth = 27;
267     }
268     else if (wcscmp(wcsMeth, L"initgn") == 0) // initgn
269     {
270         iNumInputArg = 2;
271         meth = 28;
272     }
273     else if (wcscmp(wcsMeth, L"setall") == 0) // setall
274     {
275         iNumInputArg = 5;
276         meth = 29;
277     }
278     else if (wcscmp(wcsMeth, L"advnst") == 0) // advnst
279     {
280         iNumInputArg = 2;
281         meth = 30;
282     }
283     else
284     {
285         char* pstMeth = wide_string_to_UTF8(wcsMeth);
286         Scierror(78, _("%s: Wrong method for input argument #%d: %s unknown.\n"), "grand", iStrPos + 1, pstMeth);
287         FREE(pstMeth);
288         return types::Function::Error;
289     }
290
291     // *** get arguments before methode string. ***
292     if (meth == 11 || meth == 13 || meth == 14 || meth == 16) // grand(n, "...", ...
293     {
294         if (iStrPos != 1)
295         {
296             Scierror(999, _("%s: Wrong position for input argument #%d : Must be in position %d.\n"), "grand", iStrPos + 1, 2);
297             return types::Function::Error;
298         }
299
300         if (in[iPos]->isDouble() == false)
301         {
302             Scierror(999, _("%s: Wrong type for input argument #%d : A matrix expected.\n"), "grand", iPos + 1);
303             return types::Function::Error;
304         }
305
306         types::Double* pDblTemp = in[iPos]->getAs<types::Double>();
307
308         if (pDblTemp->isScalar() == false)
309         {
310             Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos + 1);
311             return types::Function::Error;
312         }
313
314         iNumIter = (int)pDblTemp->get(0);
315         iPos++;
316     }
317     else if (meth < 21) // grand(m, n, "...", ... || grand(matrix, "...", ...
318     {
319
320         /*if (iStrPos != 1 && iStrPos != 2)
321         {
322         Scierror(999, _("%s: Wrong position for input argument #%d : Must be in position %d or %d.\n"), "grand", iStrPos + 1, 2, 3);
323         return types::Function::Error;
324         }*/
325
326         std::vector<types::Double*> vectpDblTemp;
327         for (int i = 0; i < iStrPos; i++)
328         {
329             if (in[iPos]->isDouble() == false)
330             {
331                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos + 1);
332                 return types::Function::Error;
333             }
334
335             vectpDblTemp.push_back(in[iPos]->getAs<types::Double>());
336
337             if (iStrPos == 3 && vectpDblTemp[i]->isScalar() == false)
338             {
339                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos + 1);
340                 return types::Function::Error;
341             }
342             iPos++;
343         }
344
345         //get number of dimensions to output
346         for (int i = 0; i < iStrPos; i++)
347         {
348             itab[i] = static_cast<int>(vectpDblTemp[i]->get(0));
349         }
350     }
351
352     if (in[iPos]->isString() == false)
353     {
354         Scierror(999, _("%s: Wrong type for input argument #%d : A string expected.\n"), "grand", iPos + 1);
355         return types::Function::Error;
356     }
357     iPos++; // method string has been already got.
358
359     // *** get arguments after methode string. ***
360     if (meth == 22 || meth == 25) // setgen || phr2sd
361     {
362         if (in[iPos]->isString() == false)
363         {
364             Scierror(999, _("%s: Wrong type for input argument #%d : A string expected.\n"), "grand", iPos + 1);
365             return types::Function::Error;
366         }
367
368         pStrGenOrPhr = in[iPos]->getAs<types::String>();
369
370         if (pStrGenOrPhr->isScalar() == false)
371         {
372             Scierror(999, _("%s: Wrong size for input argument #%d : Only one string expected.\n"), "grand", iPos + 1);
373             return types::Function::Error;
374         }
375     }
376     else if (meth == 16)
377     {
378         if (in.size() > 3)
379         {
380             Scierror(999, _("Missing vect for random permutation\n"));
381             return types::Function::Error;
382         }
383     }
384     else
385     {
386
387         for (int i = iPos; i < in.size(); i++)
388         {
389             if (in[i]->isDouble() == false)
390             {
391                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", iPos + 1);
392                 return types::Function::Error;
393             }
394
395             vectpDblInput.push_back(in[i]->getAs<types::Double>());
396         }
397     }
398
399     // *** perform operation in according method string and return result. ***
400
401     types::Double* pDblOut = NULL;
402
403     if (iStrPos >= 2)
404     {
405         int iProdiTab = 1;
406         for (int i = 0; i < iStrPos; i++)
407         {
408             iProdiTab = iProdiTab * itab[i];
409         }
410
411         if (iProdiTab == 0)
412         {
413             pDblOut = types::Double::Empty();
414         }
415         else if ((itab[1] != itab[0]) && (itab[0] <= -1))
416         {
417             Scierror(999, _("%s: Wrong value for input argument #%d: Positive scalar expected.\n"), "grand", 1);
418             return types::Function::Error;
419         }
420         else if ((itab[1] != itab[0]) && (itab[1] <= -1))
421         {
422             Scierror(999, _("%s: Wrong value for input argument #%d: Positive scalar expected.\n"), "grand", 2);
423             return types::Function::Error;
424         }
425         else
426         {
427             pDblOut = new types::Double(iDims, itab);
428         }
429     }
430     else
431     {
432
433         types::Double* pDblIn = in[0]->getAs<types::Double>();
434         pDblOut = new types::Double(pDblIn->getDims(), pDblIn->getDimsArray());
435
436     }
437
438     delete[] itab;
439
440     switch (meth)
441     {
442         case 1: // beta
443         {
444             double minlog   = 1.e-37;
445
446             for (int i = 0; i < 2; i++)
447             {
448                 if (vectpDblInput[i]->isScalar() == false)
449                 {
450                     delete pDblOut;
451                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
452                     return types::Function::Error;
453                 }
454
455                 if (vectpDblInput[i]->get(0) < minlog)
456                 {
457                     delete pDblOut;
458                     Scierror(999, _("%s: Wrong value for input argument #%d : At least %lf expected.\n"), "grand", iPos + 1, minlog);
459                     return types::Function::Error;
460                 }
461             }
462
463             for (int i = 0; i < pDblOut->getSize(); i++)
464             {
465                 pDblOut->set(i, C2F(genbet)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
466             }
467
468             out.push_back(pDblOut);
469             break;
470         }
471         case 2: // binomial
472         case 3: // negative binomial
473         {
474             for (int i = 0; i < 2; i++)
475             {
476                 if (vectpDblInput[i]->isScalar() == false)
477                 {
478                     delete pDblOut;
479                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
480                     return types::Function::Error;
481                 }
482             }
483
484             if (vectpDblInput[0]->get(0) < 0.0) // N
485             {
486                 delete pDblOut;
487                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive integer expected.\n"), "grand", 4);
488                 return types::Function::Error;
489             }
490
491             if (vectpDblInput[1]->get(0) < 0.0 || vectpDblInput[1]->get(0) > 1.0) // p
492             {
493                 delete pDblOut;
494                 Scierror(999, _("%s: Wrong value for input argument #%d : A value  expected.\n"), "grand", 5);
495                 return types::Function::Error;
496             }
497
498             int N = static_cast<int>(vectpDblInput[0]->get(0));
499             if (meth == 2)
500             {
501                 for (int i = 0; i < pDblOut->getSize(); i++)
502                 {
503                     pDblOut->set(i, static_cast<double>(C2F(ignbin)(&N, vectpDblInput[1]->get())));
504                 }
505             }
506             else // meth == 3
507             {
508                 for (int i = 0; i < pDblOut->getSize(); i++)
509                 {
510                     pDblOut->set(i, static_cast<double>(C2F(ignnbn)(&N, vectpDblInput[1]->get())));
511                 }
512             }
513
514             out.push_back(pDblOut);
515             break;
516         }
517         case 4: // chisquare
518         {
519             if (vectpDblInput[0]->isScalar() == false)
520             {
521                 delete pDblOut;
522                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
523                 return types::Function::Error;
524             }
525
526             if (vectpDblInput[0]->get(0) <= 0.0) // Df
527             {
528                 delete pDblOut;
529                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive no null value expected.\n"), "grand", 4);
530                 return types::Function::Error;
531             }
532
533             for (int i = 0; i < pDblOut->getSize(); i++)
534             {
535                 pDblOut->set(i, C2F(genchi)(vectpDblInput[0]->get()));
536             }
537
538             out.push_back(pDblOut);
539             break;
540         }
541         case 5: // non central chisquare
542         {
543             for (int i = 0; i < 2; i++)
544             {
545                 if (vectpDblInput[i]->isScalar() == false)
546                 {
547                     delete pDblOut;
548                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
549                     return types::Function::Error;
550                 }
551             }
552
553             if (vectpDblInput[0]->get(0) < 1.0) // Df
554             {
555                 delete pDblOut;
556                 Scierror(999, _("%s: Wrong value for input argument #%d : value >1.0 expected.\n"), "grand", 4);
557                 return types::Function::Error;
558             }
559
560             if (vectpDblInput[1]->get(0) < 0.0) // Xnon
561             {
562                 delete pDblOut;
563                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 5);
564                 return types::Function::Error;
565             }
566
567             for (int i = 0; i < pDblOut->getSize(); i++)
568             {
569                 pDblOut->set(i, C2F(gennch)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
570             }
571
572             out.push_back(pDblOut);
573             break;
574         }
575         case 6: // exponential
576         {
577             if (vectpDblInput[0]->isScalar() == false)
578             {
579                 delete pDblOut;
580                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
581                 return types::Function::Error;
582             }
583
584             if (vectpDblInput[0]->get(0) < 0.0) // Av
585             {
586                 delete pDblOut;
587                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 4);
588                 return types::Function::Error;
589             }
590
591             for (int i = 0; i < pDblOut->getSize(); i++)
592             {
593                 pDblOut->set(i, C2F(genexp)(vectpDblInput[0]->get()));
594             }
595
596             out.push_back(pDblOut);
597             break;
598         }
599         case 7: // F variance ratio
600         {
601             for (int i = 0; i < 2; i++)
602             {
603                 if (vectpDblInput[i]->isScalar() == false)
604                 {
605                     delete pDblOut;
606                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
607                     return types::Function::Error;
608                 }
609
610                 if (vectpDblInput[i]->get(0) <= 0.0) // Dfn Dfd
611                 {
612                     delete pDblOut;
613                     Scierror(999, _("%s: Wrong value for input argument #%d : Positive no null value expected.\n"), "grand", i + 4);
614                     return types::Function::Error;
615                 }
616             }
617
618             for (int i = 0; i < pDblOut->getSize(); i++)
619             {
620                 pDblOut->set(i, C2F(genf)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
621             }
622
623             out.push_back(pDblOut);
624             break;
625         }
626         case 8: // non F variance ratio
627         {
628             for (int i = 0; i < 3; i++)
629             {
630                 if (vectpDblInput[i]->isScalar() == false)
631                 {
632                     delete pDblOut;
633                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
634                     return types::Function::Error;
635                 }
636             }
637
638             if (vectpDblInput[0]->get(0) < 1.0) // Dfn
639             {
640                 delete pDblOut;
641                 Scierror(999, _("%s: Wrong value for input argument #%d : value > 1.0 expected.\n"), "grand", 4);
642                 return types::Function::Error;
643             }
644
645             if (vectpDblInput[1]->get(0) <= 0.0) // Dfd
646             {
647                 delete pDblOut;
648                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive non null value expected.\n"), "grand", 5);
649                 return types::Function::Error;
650             }
651
652             if (vectpDblInput[2]->get(0) < 0.0) // Xnon
653             {
654                 delete pDblOut;
655                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 6);
656                 return types::Function::Error;
657             }
658
659             for (int i = 0; i < pDblOut->getSize(); i++)
660             {
661                 pDblOut->set(i, C2F(gennf)(vectpDblInput[0]->get(), vectpDblInput[1]->get(), vectpDblInput[2]->get()));
662             }
663
664             out.push_back(pDblOut);
665             break;
666         }
667         case 9: // gamma
668         {
669             for (int i = 0; i < 2; i++)
670             {
671                 if (vectpDblInput[i]->isScalar() == false)
672                 {
673                     delete pDblOut;
674                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
675                     return types::Function::Error;
676                 }
677
678                 if (vectpDblInput[i]->get(0) <= 0.0)
679                 {
680                     delete pDblOut;
681                     Scierror(999, _("%s: Wrong value for input argument #%d : Positive non null value expected.\n"), "grand", i + 4);
682                     return types::Function::Error;
683                 }
684             }
685
686             for (int i = 0; i < pDblOut->getSize(); i++)
687             {
688                 //** WARNING ** : order is changed in parameters for
689                 //compatibility between Rand(...'gam',..) and cdfgam
690
691                 pDblOut->set(i, C2F(gengam)(vectpDblInput[1]->get(), vectpDblInput[0]->get()));
692             }
693
694             out.push_back(pDblOut);
695             break;
696         }
697         case 10: // Gauss Laplace (normal)
698         {
699             for (int i = 0; i < 2; i++)
700             {
701                 if (vectpDblInput[i]->isScalar() == false)
702                 {
703                     delete pDblOut;
704                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
705                     return types::Function::Error;
706                 }
707             }
708
709             if (vectpDblInput[1]->get(0) < 0.0)
710             {
711                 delete pDblOut;
712                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 5);
713                 return types::Function::Error;
714             }
715
716             for (int i = 0; i < pDblOut->getSize(); i++)
717             {
718                 pDblOut->set(i, C2F(gennor)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
719             }
720
721             out.push_back(pDblOut);
722             break;
723         }
724         case 11: // multivariate gaussian (multivariate normal)
725         {
726             if (vectpDblInput[0]->getCols() != 1 || vectpDblInput[0]->getSize() == 0)
727             {
728                 delete pDblOut;
729                 Scierror(999, _("%s: Wrong size for input argument #%d : A matrix of size m x 1 expected.(m > 0)\n"), "grand", 4);
730                 return types::Function::Error;
731             }
732
733             if (vectpDblInput[0]->getRows() != vectpDblInput[1]->getRows())
734             {
735                 delete pDblOut;
736                 Scierror(999, _("%s: Wrong size for input argument #%d and #%d: Mean and Cov have incompatible dimensions.\n"), "grand", 4, 5);
737                 return types::Function::Error;
738             }
739
740             if (vectpDblInput[1]->getRows() != vectpDblInput[1]->getCols())
741             {
742                 delete pDblOut;
743                 Scierror(999, _("%s: Wrong size for input argument #%d : A square symmetric positive definite matrix expected.\n"), "grand", 5);
744                 return types::Function::Error;
745             }
746
747             int ierr = 0;
748             int size = vectpDblInput[0]->getRows();
749             int mp   = (size * (size + 3)) / 2 + 1;
750
751             delete pDblOut;
752             types::Double* pDblOut = new types::Double(size, iNumIter);
753             double* work = new double[size];
754             double* param = new double[mp];
755
756             types::Double* pDblMean = vectpDblInput[0]->clone()->getAs<types::Double>();
757             types::Double* pDblCov  = vectpDblInput[1]->clone()->getAs<types::Double>();
758
759             C2F(setgmn)(pDblMean->get(), pDblCov->get(), &size, &size, param, &ierr);
760
761             delete pDblMean;
762             delete pDblCov;
763
764             if (ierr == 1)
765             {
766                 delete work;
767                 delete param;
768                 delete pDblOut;
769                 Scierror(999, _("%s: setgmn return with state %d.\n"), "grand", ierr);
770                 return types::Function::Error;
771             }
772
773             for (int i = 0; i < iNumIter; i++)
774             {
775                 C2F(genmn)(param, pDblOut->get() + (size * i), work);
776             }
777
778             delete work;
779             delete param;
780
781             out.push_back(pDblOut);
782             break;
783         }
784         case 12: // geometric
785         {
786             double pmin = 1.3e-307;
787
788             if (vectpDblInput[0]->isScalar() == false)
789             {
790                 delete pDblOut;
791                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
792                 return types::Function::Error;
793             }
794
795             if (vectpDblInput[0]->get(0) < pmin || vectpDblInput[0]->get(0) > 1.0)
796             {
797                 delete pDblOut;
798                 Scierror(999, _("%s: Wrong value for input argument #%d : Must be between %lf and %d.\n"), "grand", 4, pmin, 1);
799                 return types::Function::Error;
800             }
801
802             for (int i = 0; i < pDblOut->getSize(); i++)
803             {
804                 pDblOut->set(i, igngeom(vectpDblInput[0]->get(0)));
805             }
806
807             out.push_back(pDblOut);
808             break;
809         }
810         case 13: // markov
811         {
812             if ( vectpDblInput[0]->getRows() != vectpDblInput[0]->getCols() &&
813                     vectpDblInput[0]->getRows() != 1)
814             {
815                 delete pDblOut;
816                 Scierror(999, _("%s: Wrong size for input argument #%d : A square matrix or a row vector expected.\n"), "grand", 4);
817                 return types::Function::Error;
818             }
819
820             if (vectpDblInput[1]->getSize() == 0)
821             {
822                 delete pDblOut;
823                 Scierror(999, _("%s: Wrong size for input argument #%d: No empty matrix expected.\n"), "grand", 5);
824                 return types::Function::Error;
825             }
826
827             int sizeOfX0 = vectpDblInput[1]->getSize();
828
829             for ( int i = 0; i < sizeOfX0; i++)
830             {
831                 if (vectpDblInput[1]->get(i) < 1 || (vectpDblInput[1]->get(i) - 1) >= vectpDblInput[0]->getCols())
832                 {
833                     Scierror(999, _("%s: X0(%d) must be in the range [1,%d[.\n"), "grand", i + 1, vectpDblInput[0]->getCols() + 1);
834                     return types::Function::Error;
835                 }
836             }
837
838             int size = vectpDblInput[0]->getSize() + vectpDblInput[0]->getCols();
839             double* work = (double*)malloc(size * sizeof(double));
840
841             for (int i = 0; i < vectpDblInput[0]->getRows(); i++)
842             {
843                 double ptot = 0.0;
844                 for (int j = 0; j < vectpDblInput[0]->getCols(); j++)
845                 {
846                     int position = vectpDblInput[0]->getRows() * j + i;
847
848                     if (vectpDblInput[0]->get(position) < 0 || vectpDblInput[0]->get(position) > 1)
849                     {
850                         delete pDblOut;
851                         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);
852                         return types::Function::Error;
853                     }
854
855                     ptot += vectpDblInput[0]->get(position);
856                 }
857
858                 if (fabs(ptot - 1.0) > 1e-8)
859                 {
860                     delete pDblOut;
861                     Scierror(999, _("%s: Sum of P(%d,1:%d)=%lf ~= 1.\n"), "grand", i + 1, vectpDblInput[0]->getCols(), ptot);
862                     return types::Function::Error;
863                 }
864             }
865             // ** Computing the cumulative sum of the P matrix **
866
867             for (int i = 0; i < vectpDblInput[0]->getRows(); i++)
868             {
869                 double cumsum = 0.0;
870                 work[i] = 0.0;
871                 for (int j = 1; j < vectpDblInput[0]->getCols() + 1; j++)
872                 {
873                     cumsum += vectpDblInput[0]->get(vectpDblInput[0]->getRows() * (j - 1) + i);
874                     work[vectpDblInput[0]->getRows() * j + i] = cumsum;
875                 }
876             }
877
878             for (int i = 0; i < vectpDblInput[1]->getSize(); i++)
879             {
880                 int iCur = static_cast<int>(vectpDblInput[1]->get(i)) - 1;
881                 for (int j = 0; j < iNumIter; j++)
882                 {
883                     int niv = 0;
884                     double rr = C2F(ranf)();
885
886                     if (vectpDblInput[0]->getRows() == 1)
887                     {
888                         iCur = 0;
889                     }
890
891                     while (  rr >= work[iCur + vectpDblInput[0]->getRows() * niv] &&
892                              niv < (vectpDblInput[0]->getRows() + 1))
893                     {
894                         niv++;
895                     }
896
897                     // ** projection to avoid boundaries **
898                     niv = std::max(std::min(niv, vectpDblInput[0]->getCols()), 1);
899                     pDblOut->set(vectpDblInput[1]->getSize() * j + i, static_cast<double>(niv));
900                     iCur = niv - 1;
901                 }
902             }
903
904             out.push_back(pDblOut);
905             break;
906         }
907         case 14: // multinomial
908         {
909             if (vectpDblInput[0]->isScalar() == false)
910             {
911                 delete pDblOut;
912                 Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), "grand", 4);
913                 return types::Function::Error;
914             }
915
916             if (vectpDblInput[0]->get(0) < 0)
917             {
918                 delete pDblOut;
919                 Scierror(999, _("%s: Wrong value for input argument #%d: A positive scalar expected.\n"), "grand", 4);
920                 return types::Function::Error;
921             }
922
923             if (vectpDblInput[1]->getCols() != 1 || vectpDblInput[1]->getRows() <= 0)
924             {
925                 delete pDblOut;
926                 Scierror(999, _("%s: Wrong size for input argument #%d: A colomn vector expected.\n"), "grand", 5);
927                 return types::Function::Error;
928             }
929
930             double ptot = 0.0;
931             int ncat = vectpDblInput[1]->getRows() + 1;
932
933             for (int i = 0; i < vectpDblInput[1]->getCols(); i++)
934             {
935                 if (vectpDblInput[1]->get(i) < 0.0 || vectpDblInput[1]->get(i) > 1.0)
936                 {
937                     delete pDblOut;
938                     Scierror(999, _("%s: Wrong value for input argument #%d: P(%d) must be in the range [0 1].\n"), "grand", 5, i + 1);
939                     return types::Function::Error;
940                 }
941                 ptot += vectpDblInput[1]->get(i);
942             }
943
944             if (ptot > 1.0)
945             {
946                 delete pDblOut;
947                 Scierror(999, _("%s: Wrong value for input argument #%d: Sum of P(i) > 1.\n"), "grand", 5);
948                 return types::Function::Error;
949             }
950
951             int* piP = new int[vectpDblInput[0]->getSize()];
952             int* piOut = new int[pDblOut->getSize()];
953
954             for (int i = 0; i < vectpDblInput[0]->getSize(); i++)
955             {
956                 piP[i] = static_cast<int>(vectpDblInput[0]->get(i));
957             }
958
959             for (int i = 0; i < iNumIter; i++)
960             {
961                 C2F(genmul)(piP, vectpDblInput[1]->get(), &ncat,  + &piOut[ncat * i]);
962             }
963
964             for (int i = 0; i < pDblOut->getSize(); i++)
965             {
966                 pDblOut->set(i, static_cast<double>(piOut[i]));
967             }
968
969             delete piP;
970             delete piOut;
971
972             out.push_back(pDblOut);
973             break;
974         }
975         case 15: // Poisson
976         {
977             if (vectpDblInput[0]->isScalar() == false)
978             {
979                 delete pDblOut;
980                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 4);
981                 return types::Function::Error;
982             }
983
984             if (vectpDblInput[0]->get(0) < 0.0)
985             {
986                 delete pDblOut;
987                 Scierror(999, _("%s: Wrong value for input argument #%d : Positive value expected.\n"), "grand", 4);
988                 return types::Function::Error;
989             }
990
991             for (int i = 0; i < pDblOut->getSize(); i++)
992             {
993                 pDblOut->set(i, static_cast<double>(C2F(ignpoi)(vectpDblInput[0]->get())));
994             }
995
996             out.push_back(pDblOut);
997             break;
998         }
999         case 16: // random permutations
1000         {
1001             delete pDblOut;
1002             types::InternalType* pITOut = NULL;
1003
1004             // get dims and isComplex
1005             int iDims = 0;
1006             int* piDimsArray = NULL;
1007             bool bIsComplex = false;
1008             if (in[2]->isGenericType())
1009             {
1010                 types::GenericType* pGT = in[2]->getAs<types::GenericType>();
1011                 iDims = pGT->getDims();
1012                 piDimsArray = pGT->getDimsArray();
1013                 bIsComplex = pGT->isComplex();
1014             }
1015
1016             switch (in[2]->getType())
1017             {
1018                 case types::InternalType::ScilabInt8:
1019                     pITOut = new types::Int8(iDims, piDimsArray);
1020                     sci_grand_prm(iNumIter, in[2]->getAs<types::Int8>(), &pITOut);
1021                     break;
1022                 case types::InternalType::ScilabUInt8:
1023                     pITOut = new types::UInt8(iDims, piDimsArray);
1024                     sci_grand_prm(iNumIter, in[2]->getAs<types::UInt8>(), &pITOut);
1025                     break;
1026                 case types::InternalType::ScilabInt16:
1027                     pITOut = new types::Int16(iDims, piDimsArray);
1028                     sci_grand_prm(iNumIter, in[2]->getAs<types::Int16>(), &pITOut);
1029                     break;
1030                 case types::InternalType::ScilabUInt16:
1031                     pITOut = new types::UInt16(iDims, piDimsArray);
1032                     sci_grand_prm(iNumIter, in[2]->getAs<types::UInt16>(), &pITOut);
1033                     break;
1034                 case types::InternalType::ScilabInt32:
1035                     pITOut = new types::Int32(iDims, piDimsArray);
1036                     sci_grand_prm(iNumIter, in[2]->getAs<types::Int32>(), &pITOut);
1037                     break;
1038                 case types::InternalType::ScilabUInt32:
1039                     pITOut = new types::UInt32(iDims, piDimsArray);
1040                     sci_grand_prm(iNumIter, in[2]->getAs<types::UInt32>(), &pITOut);
1041                     break;
1042                 case types::InternalType::ScilabInt64:
1043                     pITOut = new types::Int64(iDims, piDimsArray);
1044                     sci_grand_prm(iNumIter, in[2]->getAs<types::Int64>(), &pITOut);
1045                     break;
1046                 case types::InternalType::ScilabUInt64:
1047                     pITOut = new types::UInt64(iDims, piDimsArray);
1048                     sci_grand_prm(iNumIter, in[2]->getAs<types::UInt64>(), &pITOut);
1049                     break;
1050                 case types::InternalType::ScilabDouble:
1051                     pITOut = new types::Double(iDims, piDimsArray, bIsComplex);
1052                     sci_grand_prm(iNumIter, in[2]->getAs<types::Double>(), &pITOut);
1053                     break;
1054                 case types::InternalType::ScilabBool:
1055                     pITOut = new types::Bool(iDims, piDimsArray);
1056                     sci_grand_prm(iNumIter, in[2]->getAs<types::Bool>(), &pITOut);
1057                     break;
1058                 case types::InternalType::ScilabString:
1059                     pITOut = new types::String(iDims, piDimsArray);
1060                     sci_grand_prm(iNumIter, in[2]->getAs<types::String>(), &pITOut);
1061                     break;
1062                 case types::InternalType::ScilabPolynom:
1063                     pITOut = new types::Polynom(in[2]->getAs<types::Polynom>()->getVariableName(), iDims, piDimsArray);
1064                     sci_grand_prm(iNumIter, in[2]->getAs<types::Polynom>(), &pITOut);
1065                     break;
1066                 case types::InternalType::ScilabSparse:
1067                 {
1068                     std::complex<double> cplxDbl;
1069                     types::InternalType* pITOutTempo = NULL;
1070                     types::Double* pDblTempo = NULL;
1071                     types::Sparse* pSP = in[2]->getAs<types::Sparse>();
1072                     int isize = pSP->getSize();
1073                     pITOut = new types::Sparse(pSP->getRows(), pSP->getCols(), pSP->isComplex());
1074                     pDblTempo = new types::Double(isize, 1, pSP->isComplex());
1075                     pITOutTempo = new types::Double(isize, iNumIter, pSP->isComplex());
1076
1077                     if (pDblTempo->isComplex())
1078                     {
1079                         for (int i = 0; i < isize; i++)
1080                         {
1081                             cplxDbl = in[2]->getAs<types::Sparse>()->getImg(i);
1082                             pDblTempo->set(i, cplxDbl.real());
1083                             pDblTempo->set(i, cplxDbl.imag());
1084                         }
1085                     }
1086                     else
1087                     {
1088                         for (int i = 0; i < isize; i++)
1089                         {
1090                             pDblTempo->set(i, in[2]->getAs<types::Sparse>()->get(i));
1091                         }
1092                     }
1093                     sci_grand_prm(iNumIter, pDblTempo, &pITOutTempo);
1094
1095                     if (pDblTempo->isComplex())
1096                     {
1097                         for (int i = 0; i < isize; i++)
1098                         {
1099                             cplxDbl.real(pITOutTempo->getAs<types::Double>()->get(i));
1100                             cplxDbl.imag(pITOutTempo->getAs<types::Double>()->getImg(i));
1101                             pITOutTempo->getAs<types::Sparse>()->set(i, cplxDbl);
1102                         }
1103                     }
1104                     else
1105                     {
1106                         for (int i = 0; i < isize; i++)
1107                         {
1108                             pITOut->getAs<types::Sparse>()->set(i, pITOutTempo->getAs<types::Double>()->get(i));
1109                         }
1110                     }
1111
1112                     delete pITOutTempo;
1113                     delete pDblTempo;
1114                     break;
1115                 }
1116                 default:
1117                 {
1118                     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");
1119                     return types::Function::Error;
1120                 }
1121             }
1122
1123             out.push_back(pITOut);
1124             break;
1125         }
1126         case 17: // uniform (def)
1127         {
1128             for (int i = 0; i < pDblOut->getSize(); i++)
1129             {
1130                 pDblOut->set(i, C2F(ranf)());
1131             }
1132
1133             out.push_back(pDblOut);
1134             break;
1135         }
1136         case 18: // uniform (unf)
1137         {
1138             for (int i = 0; i < 2; i++)
1139             {
1140                 if (vectpDblInput[i]->isScalar() == false)
1141                 {
1142                     delete pDblOut;
1143                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
1144                     return types::Function::Error;
1145                 }
1146             }
1147
1148             double low  = vectpDblInput[0]->get(0);
1149             double high = vectpDblInput[1]->get(0);
1150
1151             if (low > high)
1152             {
1153                 delete pDblOut;
1154                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low < High expected.\n"), "grand", 4, 5);
1155                 return types::Function::Error;
1156             }
1157
1158             for (int i = 0; i < pDblOut->getSize(); i++)
1159             {
1160                 pDblOut->set(i, low + (high - low) * C2F(ranf)());
1161             }
1162
1163             out.push_back(pDblOut);
1164             break;
1165         }
1166         case 19: // uniform (uin)
1167         {
1168             for (int i = 0; i < 2; i++)
1169             {
1170                 if (vectpDblInput[i]->isScalar() == false)
1171                 {
1172                     delete pDblOut;
1173                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 4);
1174                     return types::Function::Error;
1175                 }
1176             }
1177
1178             int low  = static_cast<int>(vectpDblInput[0]->get(0));
1179             int high = static_cast<int>(vectpDblInput[1]->get(0));
1180
1181             if (low > high)
1182             {
1183                 delete pDblOut;
1184                 Scierror(999, _("%s: Wrong value for input argument #%d and #%d: Low < High expected.\n"), "grand", 4, 5);
1185                 return types::Function::Error;
1186             }
1187
1188             if ( low  != vectpDblInput[0]->get(0) ||
1189                     high != vectpDblInput[1]->get(0) ||
1190                     (high - low + 1) > 2147483561)
1191             {
1192                 delete pDblOut;
1193                 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);
1194                 return types::Function::Error;
1195             }
1196
1197             for (int i = 0; i < pDblOut->getSize(); i++)
1198             {
1199                 pDblOut->set(i, C2F(ignuin)(vectpDblInput[0]->get(), vectpDblInput[1]->get()));
1200             }
1201
1202             out.push_back(pDblOut);
1203             break;
1204         }
1205         case 20: // uniform (lgi)
1206         {
1207             for (int i = 0; i < pDblOut->getSize(); i++)
1208             {
1209                 pDblOut->set(i, static_cast<double>(ignlgi()));
1210             }
1211
1212             out.push_back(pDblOut);
1213             break;
1214         }
1215         case 21: // getgen
1216         {
1217             delete pDblOut;
1218             types::String* pStrOut = new types::String(names_gen[current_base_gen]);
1219             out.push_back(pStrOut);
1220             break;
1221         }
1222         case 22: // setgen
1223         {
1224             delete pDblOut;
1225             wchar_t* wcsGen = pStrGenOrPhr->get(0);
1226
1227             if (wcscmp(wcsGen, L"mt") == 0)
1228             {
1229                 ConfigVariable::setCurrentBaseGen(MT);
1230             }
1231             else if (wcscmp(wcsGen, L"kiss") == 0)
1232             {
1233                 ConfigVariable::setCurrentBaseGen(KISS);
1234             }
1235             else if (wcscmp(wcsGen, L"clcg4") == 0)
1236             {
1237                 ConfigVariable::setCurrentBaseGen(CLCG4);
1238             }
1239             else if (wcscmp(wcsGen, L"clcg2") == 0)
1240             {
1241                 ConfigVariable::setCurrentBaseGen(CLCG2);
1242             }
1243             else if (wcscmp(wcsGen, L"urand") == 0)
1244             {
1245                 ConfigVariable::setCurrentBaseGen(URAND);
1246             }
1247             else if (wcscmp(wcsGen, L"fsultra") == 0)
1248             {
1249                 ConfigVariable::setCurrentBaseGen(FSULTRA);
1250             }
1251             else
1252             {
1253                 Scierror(999, _("%s: Wrong value for input argument #%d: '%s', '%s', '%s', '%s', '%s' or '%s' expected.\n"), "grand", 2, "mt", "kiss", "clcg4", "clcg2", "urand", "fsultra");
1254                 return types::Function::Error;
1255             }
1256
1257             int current_gen = ConfigVariable::getCurrentBaseGen();
1258             types::String* pStrOut = new types::String(names_gen[current_gen]);
1259             out.push_back(pStrOut);
1260             break;
1261         }
1262         case 23: // getsd
1263         {
1264             delete pDblOut;
1265             switch (current_base_gen)
1266             {
1267                 case MT:
1268                 {
1269                     pDblOut = new types::Double(625, 1);
1270                     get_state_mt(pDblOut->get());
1271                     break;
1272                 }
1273                 case KISS:
1274                 {
1275                     pDblOut = new types::Double(4, 1);
1276                     get_state_kiss(pDblOut->get());
1277                     break;
1278                 }
1279                 case CLCG4:
1280                 {
1281                     pDblOut = new types::Double(4, 1);
1282                     int current_clcg4 = ConfigVariable::getCurrentClcg4();
1283                     get_state_clcg4(current_clcg4, pDblOut->get());
1284                     break;
1285                 }
1286                 case CLCG2:
1287                 {
1288                     pDblOut = new types::Double(2, 1);
1289                     get_state_clcg2(pDblOut->get());
1290                     break;
1291                 }
1292                 case URAND:
1293                 {
1294                     pDblOut = new types::Double(1, 1);
1295                     get_state_urand(pDblOut->get());
1296                     break;
1297                 }
1298                 case FSULTRA:
1299                 {
1300                     pDblOut = new types::Double(40, 1);
1301                     get_state_fsultra(pDblOut->get());
1302                     break;
1303                 }
1304             }
1305
1306             out.push_back(pDblOut);
1307             break;
1308         }
1309         case 24: // setsd
1310         {
1311             delete pDblOut;
1312             int ierr = 0;
1313             switch (current_base_gen)
1314             {
1315                 case MT:
1316                 {
1317                     if (vectpDblInput[0]->isScalar())
1318                     {
1319                         ierr = set_state_mt_simple(vectpDblInput[0]->get(0));
1320                     }
1321                     else if (vectpDblInput[0]->getSize() != 625)
1322                     {
1323                         ierr = set_state_mt(vectpDblInput[0]->get());
1324                     }
1325                     else
1326                     {
1327                         Scierror(999, _("%s: Wrong size for input argument #%d : A scalar or a vector of size %d expected.\n"), "grand", 4, 625);
1328                         return types::Function::Error;
1329                     }
1330
1331                     break;
1332                 }
1333                 case KISS:
1334                 case CLCG4:
1335                 {
1336                     for (int i = 0; i < 4; i++)
1337                     {
1338                         if (vectpDblInput[i]->isScalar() == false)
1339                         {
1340                             Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 2);
1341                             return types::Function::Error;
1342                         }
1343                     }
1344
1345                     if (current_base_gen == KISS)
1346                     {
1347                         ierr = set_state_kiss(  vectpDblInput[0]->get(0),
1348                                                 vectpDblInput[1]->get(0),
1349                                                 vectpDblInput[2]->get(0),
1350                                                 vectpDblInput[3]->get(0));
1351                     }
1352                     else // CLCG4
1353                     {
1354                         ierr = set_seed_clcg4(  ConfigVariable::getCurrentClcg4(),
1355                                                 vectpDblInput[0]->get(0),
1356                                                 vectpDblInput[1]->get(0),
1357                                                 vectpDblInput[2]->get(0),
1358                                                 vectpDblInput[3]->get(0));
1359                     }
1360
1361                     break;
1362                 }
1363                 case CLCG2:
1364                 {
1365                     for (int i = 0; i < 2; i++)
1366                     {
1367                         if (vectpDblInput[i]->isScalar() == false)
1368                         {
1369                             Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 2);
1370                             return types::Function::Error;
1371                         }
1372                     }
1373
1374                     ierr = set_state_clcg2(vectpDblInput[0]->get(0), vectpDblInput[1]->get(0));
1375                     break;
1376                 }
1377                 case URAND:
1378                 {
1379                     if (vectpDblInput[0]->isScalar() == false)
1380                     {
1381                         Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1382                         return types::Function::Error;
1383                     }
1384
1385                     ierr = set_state_urand(vectpDblInput[0]->get(0));
1386                     break;
1387                 }
1388                 case FSULTRA:
1389                 {
1390                     if (in.size() == 2)
1391                     {
1392                         if (vectpDblInput[0]->getRows() != 40 || vectpDblInput[0]->getCols() != 1)
1393                         {
1394                             Scierror(999, _("%s: Wrong size for input argument #%d : A vector of size %d x %d expected.\n"), "grand", 2, 40, 1);
1395                             return types::Function::Error;
1396                         }
1397
1398                         ierr = set_state_fsultra(vectpDblInput[0]->get());
1399                     }
1400                     else // in.size() == 3
1401                     {
1402                         for (int i = 0; i < 2; i++)
1403                         {
1404                             if (vectpDblInput[i]->isScalar() == false)
1405                             {
1406                                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 2);
1407                                 return types::Function::Error;
1408                             }
1409                         }
1410
1411                         ierr = set_state_fsultra_simple(vectpDblInput[0]->get(0), vectpDblInput[1]->get(0));
1412                     }
1413
1414                     break;
1415                 }
1416             }
1417
1418             if (ierr == 0)
1419             {
1420                 Scierror(999, _("%s: Wrong value for the last %d input argument(s).\n"), "grand", in.size() - 1);
1421                 return types::Function::Error;
1422             }
1423
1424             break;
1425         }
1426         case 25: // phr2sd
1427         {
1428             delete pDblOut;
1429             if (pStrGenOrPhr->isScalar() == false)
1430             {
1431                 Scierror(999, _("%s: Wrong type for input argument #%d : One string expected.\n"), "grand", 2);
1432                 return types::Function::Error;
1433             }
1434
1435             types::Double* pDblOut = new types::Double(1, 2);
1436             int size = (int)wcslen(pStrGenOrPhr->get(0));
1437             int piOut[2];
1438             char* strPhr = wide_string_to_UTF8(pStrGenOrPhr->get(0));
1439
1440             C2F(phrtsd)(strPhr, &size, piOut, piOut + 1, size);
1441
1442             pDblOut->set(0, static_cast<double>(piOut[0]));
1443             pDblOut->set(1, static_cast<double>(piOut[1]));
1444
1445             FREE(strPhr);
1446             out.push_back(pDblOut);
1447             break;
1448         }
1449         case 26: // setcgn
1450         {
1451             delete pDblOut;
1452             if (current_base_gen != CLCG4)
1453             {
1454                 sciprint(_("The %s option affects only the %s generator\n"), "setcgn", "clcg4");
1455             }
1456
1457             if (vectpDblInput[0]->isScalar() == false)
1458             {
1459                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1460                 return types::Function::Error;
1461             }
1462
1463             if (vectpDblInput[0]->get(0) < 0 || vectpDblInput[0]->get(0) > Maxgen)
1464             {
1465                 Scierror(999, _("%s: Wrong value for input argument #%d : Must be between %d and %d.\n"), "grand", 0, Maxgen);
1466                 return types::Function::Error;
1467             }
1468
1469             ConfigVariable::setCurrentClcg4(static_cast<int>(vectpDblInput[0]->get(0)));
1470             double dOut = static_cast<double>(ConfigVariable::getCurrentClcg4());
1471             out.push_back(new types::Double(dOut));
1472             break;
1473         }
1474         case 27: // getcgn
1475         {
1476             delete pDblOut;
1477             double dOut = static_cast<double>(ConfigVariable::getCurrentClcg4());
1478             out.push_back(new types::Double(dOut));
1479             break;
1480         }
1481         case 28: // initgn
1482         {
1483             delete pDblOut;
1484             SeedType where;
1485             if (current_base_gen != CLCG4)
1486             {
1487                 sciprint(_("The %s option affects only the %s generator\n"), "initgn", "clcg4");
1488             }
1489
1490             if (vectpDblInput[0]->isScalar() == false)
1491             {
1492                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1493                 return types::Function::Error;
1494             }
1495
1496             if ( vectpDblInput[0]->get(0) != 0 &&
1497                     vectpDblInput[0]->get(0) != -1 &&
1498                     vectpDblInput[0]->get(0) != 1)
1499             {
1500                 Scierror(999, _("%s: Wrong value for input argument #%d : Must be between %d, %d or %d.\n"), "grand", 2, -1, 0, 1);
1501                 return types::Function::Error;
1502             }
1503
1504             where = (SeedType)(int)(vectpDblInput[0]->get(0) + 1);
1505             init_generator_clcg4(ConfigVariable::getCurrentClcg4(), where);
1506             out.push_back(vectpDblInput[0]);
1507             break;
1508         }
1509         case 29: // setall
1510         {
1511             delete pDblOut;
1512             if (current_base_gen != CLCG4)
1513             {
1514                 sciprint(_("The %s option affects only the %s generator\n"), "setall", "clcg4");
1515             }
1516
1517             for (int i = 0; i < 4; i++)
1518             {
1519                 if (vectpDblInput[i]->isScalar() == false)
1520                 {
1521                     Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 2);
1522                     return types::Function::Error;
1523                 }
1524             }
1525
1526             int ierr = set_initial_seed_clcg4(  vectpDblInput[0]->get(0),
1527                                                 vectpDblInput[1]->get(0),
1528                                                 vectpDblInput[2]->get(0),
1529                                                 vectpDblInput[3]->get(0));
1530             if (ierr == 0)
1531             {
1532                 Scierror(999, _("%s: Wrong value for the last %d input argument(s).\n"), "grand", 4);
1533                 return types::Function::Error;
1534             }
1535
1536             out.push_back(pStrMethod);
1537             break;
1538         }
1539         case 30: // advnst
1540         {
1541             delete pDblOut;
1542             if (current_base_gen != CLCG4)
1543             {
1544                 sciprint(_("The %s option affects only the %s generator\n"), "advnst", "clcg4");
1545             }
1546
1547             if (vectpDblInput[0]->isScalar() == false)
1548             {
1549                 Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", 2);
1550                 return types::Function::Error;
1551             }
1552
1553             int k = static_cast<int>(vectpDblInput[0]->get(0));
1554
1555             if (k < 1)
1556             {
1557                 Scierror(999, _("%s: Wrong value for input argument #%d : Must be > %d.\n"), "grand", 2, 0);
1558                 return types::Function::Error;
1559             }
1560
1561             advance_state_clcg4(ConfigVariable::getCurrentClcg4(), k);
1562             out.push_back(new types::Double(static_cast<double>(k)));
1563             break;
1564         }
1565     }
1566
1567     return types::Function::OK;
1568 }
1569 /*--------------------------------------------------------------------------*/
1570
1571 template<class U>
1572 void sci_grand_prm(int iNumIter, U *pIn, types::InternalType** pOut)
1573 {
1574     U* pUTempo = NULL;
1575     types::InternalType* pITTempo = NULL;
1576     int* piDimsArray = NULL;
1577     int Dims = 0;
1578
1579     if ((pIn->getCols() == 1) && (pIn->getDims() == 2))
1580     {
1581         pOut[0]->getAs<U>()->resize(pIn->getRows(), iNumIter);
1582         pUTempo = pIn;
1583     }
1584     else if ((pIn->getRows() == 1) && (pIn->getDims() == 2))
1585     {
1586         pIn->transpose(pITTempo);
1587         pOut[0]->getAs<U>()->resize(iNumIter, pIn->getCols());
1588         pUTempo = pITTempo->getAs<U>();
1589     }
1590     else
1591     {
1592         piDimsArray = pOut[0]->getAs<U>()->getDimsArray();
1593         Dims = pOut[0]->getAs<U>()->getDims();
1594         piDimsArray[Dims] = iNumIter;
1595         pOut[0]->getAs<U>()->resize(piDimsArray, Dims + 1);
1596         pUTempo = pIn;
1597     }
1598
1599     int isize = pUTempo->getSize();
1600
1601     types::Double* pDblOut = new types::Double(isize, iNumIter, pUTempo->isComplex());
1602
1603     for (int i = 0; i < iNumIter; i++)
1604     {
1605         for (int j = 0; j < isize; j++)
1606         {
1607             pDblOut->set(isize * i + j, j);
1608         }
1609         C2F(genprm)(pDblOut->get() + (isize * i), &isize);
1610     }
1611
1612     if ((pIn->getCols() != 1) && (pIn->getRows() == 1) && (pIn->getDims() == 2))
1613     {
1614         pDblOut->transpose(pITTempo);
1615         delete pDblOut;
1616         pDblOut = pITTempo->getAs<types::Double>();
1617     }
1618
1619     if (pUTempo->isComplex() && pUTempo->isPoly() == false)
1620     {
1621         for (int i = 0; i < pOut[0]->getAs<U>()->getSize(); i++)
1622         {
1623             pOut[0]->getAs<U>()->set(i , pIn->get(pDblOut->get(i)));
1624             pOut[0]->getAs<U>()->setImg(i , pIn->getImg(pDblOut->get(i)));
1625         }
1626     }
1627     else
1628     {
1629         for (int i = 0; i < pOut[0]->getAs<U>()->getSize(); i++)
1630         {
1631             pOut[0]->getAs<U>()->set(i, pIn->get(pDblOut->get(i)));
1632         }
1633     }
1634
1635     if ((pIn->getCols() != 1) && (pIn->getRows() == 1) && (pIn->getDims() == 2))
1636     {
1637         delete pUTempo;
1638     }
1639
1640     delete pDblOut;
1641 }
1642
1643 /*--------------------------------------------------------------------------*/