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