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