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