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