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