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