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