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