1 /* -*- mode: C; mode: fold -*-
2 Copyright (C) 2007-2017,2018 John E. Davis
3
4 This file is part of the S-Lang Library.
5
6 The S-Lang Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 2 of the
9 License, or (at your option) any later version.
10
11 The S-Lang Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this library; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
19 USA.
20 */
21 #include "config.h"
22
23 #include <stdio.h>
24 #include <math.h>
25 #include <slang.h>
26
27 #ifdef HAVE_PROCESS_H
28 # include <process.h> /* for getpid */
29 #endif
30
31 #include <sys/types.h>
32
33 #ifdef HAVE_UNISTD_H
34 # include <unistd.h>
35 #endif
36
37 #include <time.h>
38
39 /* The random number generator used here uses a combination of 3 generators.
40 * One of the generators was derived from mzran13 published in
41 *
42 * G. Marsaglia and A. Zaman, ``Some portable very-long-period
43 * random number generators'', Computers in Physics, Vol. 8,
44 * No. 1, Jan/Feb 1994
45 *
46 * However, the algorithm (mzran13) published there contained a few
47 * typos and as written was not very good. The subtract-with-carry
48 * component of mzran13 adopted for the purposes of the generator in
49 * this file has a period of 4x10^28.
50 *
51 * The other two generators are described at
52 * <http://www.helsbreth.org/random/rng_combo.html>. Briefly, one is
53 * a simple product generator, r_{k+2} = r_{k+1} r_k, and the other
54 * is is a multiply with carry generator. The period of this
55 * combination exceeds 10^18.
56 *
57 * The total period of the combined generators exceeds 10^46.
58 *
59 * The generator was tested using dieharder-2.24.4, and the testing
60 * indicates that it may perform a bit better than the Mersenne
61 * twister mt19937_1999, which is thought to be a high-quality
62 * generator. The speed is about the same.
63 *
64 * Note: mt19937_1999 is available to S-Lang users via the GSL module.
65 */
66
67 SLANG_MODULE(rand);
68
69 #ifdef PI
70 # undef PI
71 #endif
72 #define PI 3.14159265358979323846264338327950288
73 #define LOG_SQRT_2PI 0.9189385332046727417803297364
74
75 /*{{{ The basic generator */
76
77 #define SEED_X 521288629U
78 #define SEED_Y 362436069U
79 #define SEED_Z 16163801U
80 #define SEED_C 1 /* (SEED_Y > SEED_Z) */
81
82 #if (SIZEOF_INT == 4)
83 typedef unsigned int uint32;
84 # define UINT32_TYPE SLANG_UINT_TYPE
85 # define PUSH_UINT32 SLang_push_uint
86 #else
87 typedef unsigned long uint32
88 # define UINT32_TYPE SLANG_ULONG_TYPE
89 # define PUSH_UINT32 SLang_push_ulong
90 #endif
91
92 #define CACHE_SIZE 4 /* do not change this */
93
94 typedef struct
95 {
96 int cache_index;
97 uint32 cache[CACHE_SIZE];
98 uint32 x, y, z;
99 uint32 cx, cy, cz;
100
101 /* Gaussian Random fields */
102 int one_available;
103 double g2;
104 }
105 Rand_Type;
106
107 static uint32 generate_uint32_random (Rand_Type *);
108
109 #define NUM_SEEDS 3
seed_random(Rand_Type * rt,unsigned long seeds[NUM_SEEDS])110 static void seed_random (Rand_Type *rt, unsigned long seeds[NUM_SEEDS])
111 {
112 int count;
113 unsigned long s1 = seeds[0];
114 unsigned long s2 = seeds[1];
115 unsigned long s3 = seeds[2];
116
117 rt->x = (uint32) (s1 + SEED_X);
118 rt->y = (uint32) (s1/2 + SEED_Y);
119 rt->z = (uint32) (2*s1+SEED_Z);
120 rt->x += (rt->y > rt->z);
121
122 rt->cache_index = CACHE_SIZE;
123
124 rt->cx = s2 * 8 + 3;
125 rt->cy = s2 * 2 + 1;
126 rt->cz = s3 | 1;
127 count = 32;
128 while (count)
129 {
130 count--;
131 (void) generate_uint32_random (rt);
132 }
133
134 rt->one_available = 0;
135 rt->g2 = 0.0;
136 }
137
generate_uint32_random(Rand_Type * rt)138 static uint32 generate_uint32_random (Rand_Type *rt)
139 {
140 uint32 s0, s1, s2, s3, r;
141 uint32 *cache;
142
143 cache = rt->cache;
144 if (rt->cache_index < CACHE_SIZE)
145 return cache[rt->cache_index++];
146
147 s1 = rt->x;
148 s2 = rt->y;
149 s3 = rt->z;
150
151 r = s0 = (s2-s1) - 18*(s2<=s1); s2 += (s2 <= s1);
152 rt->x = s1 = (s3-s2) - 18*(s3<=s2); s3 += (s3 <= s2);
153 rt->y = s2 = (s0-s3) - 18*(s0<=s3); s0 += (s0 <= s3);
154 rt->z = s3 = (s1-s0) - 18*(s1<=s0); s1 += (s1 <= s0);
155
156 #define COMBINE(y,x,a,b) y = ((x)+(a)+(b))
157 s1 = rt->cx;
158 s2 = rt->cy;
159 s3 = rt->cz;
160
161 s0 = s1*s2;
162 s3 = 30903U*(s3 & 0xFFFFU) + (s3 >> 16);
163 COMBINE(r,r,s0,s3);
164
165 s1 = s2*s0;
166 s3 = 30903U*(s3 & 0xFFFFU) + (s3 >> 16);
167 COMBINE(cache[1],rt->x,s1,s3);
168
169 rt->cx = s2 = s0*s1;
170 s3 = 30903U*(s3 & 0xFFFFU) + (s3 >> 16);
171 COMBINE(cache[2],rt->y,s2,s3);
172
173 rt->cy = s0 = s1*s2;
174 rt->cz = s3 = 30903U*(s3 & 0xFFFFU) + (s3 >> 16);
175 COMBINE(cache[3],rt->z,s0,s3);
176
177 rt->cache_index = 1;
178 return r;
179 }
180
free_random(Rand_Type * r)181 static void free_random (Rand_Type *r)
182 {
183 SLfree ((char *) r);
184 }
185
create_random(unsigned long seeds[NUM_SEEDS])186 static Rand_Type *create_random (unsigned long seeds[NUM_SEEDS])
187 {
188 Rand_Type *rt;
189
190 if (NULL != (rt = (Rand_Type *) SLmalloc (sizeof (Rand_Type))))
191 seed_random (rt, seeds);
192 return rt;
193 }
194
195 /*}}}*/
196
197 /*{{{ Flat distribution */
198
generate_random_doubles(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)199 static void generate_random_doubles (Rand_Type *rt, VOID_STAR ap,
200 SLuindex_Type num, VOID_STAR parms)
201 {
202 double *a = (double *)ap;
203 double *amax = a + num;
204 uint32 *cache = rt->cache;
205
206 (void) parms;
207
208 while (a < amax)
209 {
210 uint32 u;
211
212 if (rt->cache_index < CACHE_SIZE)
213 u = cache[rt->cache_index++];
214 else
215 u = generate_uint32_random (rt);
216
217 *a++ = u / 4294967296.0;
218 }
219 }
220
generate_random_open_doubles(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)221 static void generate_random_open_doubles (Rand_Type *rt, VOID_STAR ap,
222 SLuindex_Type num, VOID_STAR parms)
223 {
224 double *a = (double *)ap;
225 double *amax = a + num;
226 uint32 *cache = rt->cache;
227
228 (void) parms;
229
230 while (a < amax)
231 {
232 uint32 u;
233
234 if (rt->cache_index < CACHE_SIZE)
235 u = cache[rt->cache_index++];
236 else
237 u = generate_uint32_random (rt);
238
239 if (u == 0)
240 continue;
241
242 *a++ = u / 4294967296.0;
243 }
244 }
245
generate_random_uints(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)246 static void generate_random_uints (Rand_Type *rt, VOID_STAR ap,
247 SLuindex_Type num, VOID_STAR parms)
248 {
249 uint32 *a = (uint32 *)ap;
250 uint32 *amax = a + num;
251 uint32 *cache = rt->cache;
252
253 (void) parms;
254
255 while (a < amax)
256 {
257 if (rt->cache_index < CACHE_SIZE)
258 *a++ = cache[rt->cache_index++];
259 else
260 *a++ = generate_uint32_random (rt);
261 }
262 }
263
264 /* produces a random number in the open interval (0,1) */
open_interval_random(Rand_Type * rt)265 static double open_interval_random (Rand_Type *rt)
266 {
267 uint32 u;
268
269 do
270 {
271 if (rt->cache_index < CACHE_SIZE)
272 u = rt->cache[rt->cache_index++];
273 else
274 u = generate_uint32_random (rt);
275 }
276 while (u == 0);
277
278 return u / 4294967296.0;
279 }
280
uniform_random(Rand_Type * rt)281 static double uniform_random (Rand_Type *rt)
282 {
283 uint32 u;
284
285 if (rt->cache_index < CACHE_SIZE)
286 u = rt->cache[rt->cache_index++];
287 else
288 u = generate_uint32_random (rt);
289
290 return u / 4294967296.0;
291 }
292
293 /*}}}*/
294
295 /*{{{ Gaussian Distribution */
296
gaussian_box_muller(Rand_Type * rt)297 static double gaussian_box_muller (Rand_Type *rt)
298 {
299 double s, g1, g2, g;
300
301 if (rt->one_available)
302 {
303 rt->one_available = 0;
304 return rt->g2;
305 }
306
307 do
308 {
309 uint32 u;
310 if (rt->cache_index < CACHE_SIZE)
311 u = rt->cache[rt->cache_index++];
312 else
313 u = generate_uint32_random (rt);
314
315 g1 = 2.0*(u/4294967296.0) - 1.0;
316
317 if (rt->cache_index < CACHE_SIZE)
318 u = rt->cache[rt->cache_index++];
319 else
320 u = generate_uint32_random (rt);
321 g2 = 2.0*(u/4294967296.0) - 1.0;
322 g = g1*g1 + g2*g2;
323 }
324 while ((g >= 1.0) || (g == 0.0));
325
326 s = sqrt (-2.0 * log (g) / g);
327 rt->g2 = g2 * s;
328 rt->one_available = 1;
329 return g1*s;
330 }
331
generate_gaussian_randoms(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)332 static void generate_gaussian_randoms (Rand_Type *rt, VOID_STAR ap,
333 SLuindex_Type num, VOID_STAR parms)
334 {
335 double *a = (double *)ap;
336 double *amax = a + num;
337 double sigma = *(double *)parms;
338
339 if ((a < amax) && (rt->one_available))
340 {
341 *a++ = sigma*rt->g2;
342 rt->one_available = 0;
343 }
344
345 while (a < amax)
346 {
347 *a++ = sigma*gaussian_box_muller (rt);
348 if (a == amax)
349 break;
350
351 *a++ = sigma*rt->g2;
352 rt->one_available = 0;
353 }
354 }
355
356 /*}}}*/
357
358 /*{{{ Poisson Distribution */
359
360 /*
361 * The algorithm used is described in:
362 * W. Hörmann "The Transformed Rejection Method for Generating
363 * Poisson Random Variables", which is available from
364 * <http://statistik.wu-wien.ac.at/papers/92-04-13.wh.ps.gz>.
365 */
366 #define FACTORIAL_TABLE_SIZE 10
367 static double Log_Factorial_Table [FACTORIAL_TABLE_SIZE+1];
368
init_poisson(void)369 static void init_poisson (void)
370 {
371 unsigned int i;
372 double x = 1.0;
373 Log_Factorial_Table[0] = 0.0;
374 for (i = 1; i <= FACTORIAL_TABLE_SIZE; i++)
375 {
376 x *= i;
377 Log_Factorial_Table[i] = log(x);
378 }
379 }
380
381 /* Assumes x >= 0 */
log_factorial(double x)382 static double log_factorial (double x)
383 {
384 double x2;
385
386 if (x <= FACTORIAL_TABLE_SIZE)
387 return Log_Factorial_Table[(unsigned int) x];
388
389 x2 = x*x;
390
391 return LOG_SQRT_2PI + (x+0.5)*log(x) - x
392 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0/x2)/x2)/x2)/x2)/x/166320.0;
393 }
394
hoermann_ptrd_poisson(Rand_Type * rt,double mu,double a,double b,double vr,double alphainv,double lnmu,double smu)395 static unsigned int hoermann_ptrd_poisson
396 (Rand_Type *rt, double mu, double a, double b, double vr,
397 double alphainv, double lnmu, double smu)
398 {
399 while (1)
400 {
401 double u, v, fk, us;
402 unsigned int k;
403
404 v = open_interval_random (rt);
405 if (v <= 0.86*vr)
406 {
407 u = v/vr - 0.43;
408 fk = floor((2.0*a/(0.5-fabs(u))+b)*u + mu + 0.445);
409 return (unsigned int)fk;
410 }
411 if (v >= vr)
412 u = open_interval_random (rt) - 0.5;
413 else
414 {
415 u = v/vr - 0.93; u = ((u < 0.0) ? -0.5 : 0.5) - u;
416 v = vr*open_interval_random (rt);
417 }
418 us = 0.5 - fabs(u);
419 if ((us < 0.013) && (v > us))
420 continue;
421
422 fk = floor((2.0*a/us + b)*u + mu + 0.445);
423 if (fk < 0.0)
424 continue;
425
426 k = (unsigned int) fk;
427
428 v = v * alphainv/(a/(us*us) + b);
429 if ((k >= 10)
430 && (log(v*smu) <= (fk+0.5)*log(mu/fk)-mu-LOG_SQRT_2PI+fk
431 - (1.0/12-1.0/(360.0*fk*fk))/fk))
432 return (unsigned int)fk;
433
434 if ((k <= 9)
435 && (log(v)<=fk*lnmu-mu-Log_Factorial_Table[k]))
436 return k;
437 }
438 }
439
knuth_poisson(Rand_Type * rt,double cutoff)440 static unsigned int knuth_poisson (Rand_Type *rt, double cutoff)
441 {
442 double x = 1.0;
443 unsigned int n = 0;
444
445 do
446 {
447 if (rt->cache_index < CACHE_SIZE)
448 x *= rt->cache[rt->cache_index++] / 4294967296.0;
449 else
450 x *= generate_uint32_random (rt) / 4294967296.0;
451 n++;
452 }
453 while (x >= cutoff);
454
455 return n - 1;
456 }
457
generate_poisson_randoms(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)458 static void generate_poisson_randoms (Rand_Type *rt, VOID_STAR ap,
459 SLuindex_Type num, VOID_STAR parms)
460 {
461 unsigned int *x = (unsigned int *)ap;
462 unsigned int *xmax = x + num;
463 double mu = *(double *)parms;
464 double cutoff;
465
466 if (mu > 10.0)
467 {
468 double smu = sqrt(mu);
469 double b = 0.931 + 2.53*smu;
470 double a = -0.059 + 0.02483*b;
471 double vr = 0.9277 - 3.6224/(b-2.0);
472 double alphainv = 1.1239 + 1.1328/(b-3.4);
473 double lnmu = log(mu);
474
475 while (x < xmax)
476 {
477 *x++ = hoermann_ptrd_poisson (rt, mu, a, b, vr, alphainv, lnmu, smu);
478 }
479 return;
480 }
481
482 cutoff = exp(-mu);
483 while (x < xmax)
484 {
485 *x++ = knuth_poisson (rt, cutoff);
486 }
487 }
488
489 /*}}}*/
490
491 /*{{{ Gamma Distribution */
492
493 /* The gamma distribution is:
494 *
495 * f(x; k,theta)dx = x^{k-1} \frac{e^{-x/theta}}{\theta^k \Gamma(k)} dx
496 * where x>0, k,theta>0.
497 * The marsaglia_tsang_gamma algorithm is for x>0, k>=1, theta=1:
498 * f(y;k) = y^{k-1}\frac{e^{-y}}{\Gamma(k)}
499 * Let x = theta*y. dy = dx/theta
500 *
501 * Then:
502 *
503 * f(y;k)dy = (x/theta)^{k-1}\frac{e^{-x/theta}}{\Gamma(k)} dx / \theta
504 * = x^{k-1}\frac{e^{-x/theta}}{\theta^k\Gamma(k)} dx
505 * = f(x; k,theta)
506 *
507 * Also the final note in the Marsaglia paper says that values for
508 * k<1 may be obtained using X_k = X_{k+1}U^{1/k} where X is a random
509 * value produced by the algorithm, and U is a uniform random number
510 * on (0,1).
511 *
512 * Marsaglia, G. and Tsang, W. W. 2000a. A simple method for generating
513 * gamma variables. ACM Transactions on Mathematical Software 26, 363–372.
514 * <http://oldmill.uchicago.edu/~wilder/Code/random/Papers/Marsaglia_00_SMGGV.pdf>
515 *
516 */
marsaglia_tsang_gamma_internal(Rand_Type * rt,double c,double d)517 static double marsaglia_tsang_gamma_internal (Rand_Type *rt, double c, double d)
518 {
519 while (1)
520 {
521 double x, u, v;
522 do
523 {
524 if (rt->one_available)
525 {
526 x = rt->g2;
527 rt->one_available = 0;
528 }
529 else x = gaussian_box_muller (rt);
530 v = 1.0 + c*x;
531 }
532 while (v <= 0.0);
533
534 v = v*v*v;
535 u = open_interval_random (rt);
536 x = x*x;
537
538 if ((u < 1.0 - 0.0331*(x*x))
539 || (log (u) < 0.5*x + d*(1.0 - v + log (v))))
540 return d*v;
541 }
542 }
543
544 /* k > 0. theta > 0 */
rand_gamma(Rand_Type * rt,double k,double theta)545 static double rand_gamma (Rand_Type *rt, double k, double theta)
546 {
547 double c, d;
548
549 #ifdef HAVE_ISNAN
550 if (isnan(k) || isnan(theta))
551 return k*theta;
552 #endif
553
554 if (k < 1.0)
555 {
556 d = k + 2.0/3.0;
557 c = (1.0/3.0)/sqrt(d);
558 return theta * marsaglia_tsang_gamma_internal (rt, c, d)
559 * pow (open_interval_random(rt), 1.0/k);
560 }
561
562 d = k - 1.0/3.0;
563 c = (1.0/3.0)/sqrt(d);
564 return theta * marsaglia_tsang_gamma_internal (rt, c, d);
565 }
566
generate_gamma_randoms(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)567 static void generate_gamma_randoms (Rand_Type *rt, VOID_STAR ap,
568 SLuindex_Type num, VOID_STAR parms)
569 {
570 double *x = (double *)ap;
571 double *xmax = x + num;
572 double k, theta;
573 double c, d;
574
575 k = ((double *)parms)[0];
576 theta = ((double *)parms)[1];
577
578 #ifdef HAVE_ISNAN
579 if (isnan(k) || isnan(theta))
580 {
581 while (x < xmax)
582 *x++ = k*theta;
583 return;
584 }
585 #endif
586
587 if (k < 1.0)
588 {
589 double kinv = 1.0/k;
590 d = k + 2.0/3.0;
591 c = (1.0/3.0)/sqrt(d);
592
593 while (x < xmax)
594 {
595 *x++ = theta * marsaglia_tsang_gamma_internal (rt, c, d)
596 * pow (open_interval_random(rt), kinv);
597 }
598 return;
599 }
600
601 d = k - 1.0/3.0;
602 c = (1.0/3.0)/sqrt(d);
603 while (x < xmax)
604 *x++ = theta * marsaglia_tsang_gamma_internal (rt, c, d);
605 }
606
607 /*}}}*/
608
609 /*{{{ Beta Distribution */
610
knuth_beta(Rand_Type * rt,double alpha,double beta)611 static double knuth_beta (Rand_Type *rt, double alpha, double beta)
612 {
613 if (0.0 == (alpha = rand_gamma (rt, alpha, 1.0)))
614 return 0.0;
615
616 beta = rand_gamma (rt, beta, 1.0);
617
618 return alpha/(alpha+beta);
619 }
620
generate_beta_randoms(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)621 static void generate_beta_randoms (Rand_Type *rt, VOID_STAR ap,
622 SLuindex_Type num, VOID_STAR parms)
623 {
624 double *x = (double *)ap;
625 double *xmax = x + num;
626 double alpha, beta;
627
628 alpha = ((double *)parms)[0];
629 beta = ((double *)parms)[1];
630
631 while (x < xmax)
632 *x++ = knuth_beta (rt, alpha, beta);
633 }
634
635 /*}}}*/
636
637 /*{{{ Binomial Distribution */
638
639 /* This algorithm is from:
640 * Hormann, W. (1993), The generation of binomial random variates,
641 * Journal of Statistical Computation and Simulation 46, 101-110.
642 * I obtained it from the preprint, but not the actual journal article.
643 * The preprint contains a typo in the algorithm that I managed to
644 * find and correct. It pays to read the paper.
645 */
646
647 typedef struct
648 {
649 double a, b, c, vr, alpha, lpq, fm, h, p;
650 unsigned int n;
651 }
652 BTRS_Type;
653
init_btrs(BTRS_Type * btrs,double p,unsigned int n)654 static void init_btrs (BTRS_Type *btrs, double p, unsigned int n)
655 {
656 double spq = sqrt (n*p*(1.0-p));
657
658 btrs->p = p;
659 btrs->n = n;
660
661 btrs->b = 1.15 + 2.53*spq;
662 btrs->a = -0.0873 + 0.0248*btrs->b + 0.01*p;
663 btrs->c = n*p+0.5;
664 btrs->vr = 0.92 - 4.2/btrs->b;
665 btrs->alpha = (2.83+5.1/btrs->b) * spq;
666 btrs->lpq = log (p/(1.0-p));
667 btrs->fm = floor ((n+1)*p);
668 btrs->h = log_factorial (btrs->fm) + log_factorial(n-btrs->fm);
669 }
670
671 /* This algorithm assumes that p <= 0.5 and p*n >= 10.0 */
binomial_btrs(Rand_Type * rt,BTRS_Type * btrs)672 static double binomial_btrs (Rand_Type *rt, BTRS_Type *btrs)
673 {
674 double a = btrs->a;
675 double b = btrs->b;
676 double c = btrs->c;
677 double vr = btrs->vr;
678 double h = btrs->h;
679 double lpq = btrs->lpq;
680 double fm = btrs->fm;
681 double alpha = btrs->alpha;
682 unsigned int n = btrs->n;
683
684 while (1)
685 {
686 double u = open_interval_random (rt) - 0.5;
687 double v = open_interval_random (rt);
688 double us = 0.5 - fabs(u);
689 double fk = floor ((2.0*a/us + b)*u + c);
690
691 if ((fk < 0.0) || ((unsigned int) fk > n))
692 continue;
693
694 if ((us >= 0.07) && (v <= vr))
695 return (unsigned int) fk;
696
697 v = log (v*alpha/(a/(us*us) + b));
698 if (v <= (h - log_factorial(fk) - log_factorial(n-fk)
699 + (fk-fm)*lpq))
700 return (unsigned int) fk;
701 }
702 }
703
704 typedef struct
705 {
706 unsigned int n;
707 double p;
708 }
709 Binomial_Parms_Type;
710
generate_binomial_randoms(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)711 static void generate_binomial_randoms (Rand_Type *rt, VOID_STAR ap,
712 SLuindex_Type num, VOID_STAR parms)
713 {
714 unsigned int *x = (unsigned int *)ap;
715 unsigned *xmax = x + num;
716 Binomial_Parms_Type *s = (Binomial_Parms_Type *)parms;
717 double p, qn, r, q, g;
718 unsigned int n;
719 int swapped = 0;
720
721 n = s->n;
722 p = s->p;
723
724 if (p > 0.5)
725 {
726 p = 1.0 - p;
727 swapped = 1;
728 }
729
730 if (n*p > 10.0)
731 {
732 BTRS_Type btrs;
733 init_btrs (&btrs, p, n);
734 if (swapped)
735 {
736 while (x < xmax)
737 *x++ = n - binomial_btrs (rt, &btrs);
738 }
739 else
740 {
741 while (x < xmax)
742 *x++ = binomial_btrs (rt, &btrs);
743 }
744 return;
745 }
746
747 /* Inverse CDF method (Adapted from RANLIB) */
748 q = 1.0 - p;
749 qn = pow (q, n);
750 r = p/q;
751 g = r*(n+1);
752 while (x < xmax)
753 {
754 #define MAX_INV_BINOMIAL_CDF_LOOPS 110
755 double f = qn;
756 double u = uniform_random (rt);
757 unsigned int k = 0;
758 unsigned kmax = (n > MAX_INV_BINOMIAL_CDF_LOOPS
759 ? MAX_INV_BINOMIAL_CDF_LOOPS : n);
760
761 while (k <= kmax)
762 {
763 if (u < f)
764 {
765 if (swapped)
766 k = n - k;
767 *x++ = k;
768 break;
769 }
770 u -= f;
771 k++;
772 f *= (g/k - r);
773 }
774 }
775 }
776
777 /*}}}*/
778
generate_cauchy_randoms(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)779 static void generate_cauchy_randoms (Rand_Type *rt, VOID_STAR ap,
780 SLuindex_Type num, VOID_STAR parms)
781 {
782 double *x = (double *)ap;
783 double *xmax = x + num;
784 double a = *(double *)parms;
785
786 while (x < xmax)
787 {
788 double u;
789 do
790 {
791 u = uniform_random (rt);
792 }
793 while (u == 0.5);
794 *x++ = a * tan (PI*u);
795 }
796 }
797
generate_geometric_randoms(Rand_Type * rt,VOID_STAR ap,SLuindex_Type num,VOID_STAR parms)798 static void generate_geometric_randoms (Rand_Type *rt, VOID_STAR ap,
799 SLuindex_Type num, VOID_STAR parms)
800 {
801 unsigned int *x = (unsigned int *)ap;
802 unsigned int *xmax = x + num;
803 double a = *(double *)parms;
804
805 if (a == 1.0)
806 {
807 while (x < xmax)
808 *x++ = 1;
809 return;
810 }
811
812 a = 1.0 / log (1.0-a); /* is negative */
813
814 while (x < xmax)
815 *x++ = (unsigned int) (1.0 + a*log(open_interval_random (rt)));
816 }
817
818 /* Interpreter Interface */
819
820 static Rand_Type *Default_Rand;
821 static int Rand_Type_Id = -1;
822
pop_seeds(unsigned long seeds[NUM_SEEDS])823 static int pop_seeds (unsigned long seeds[NUM_SEEDS])
824 {
825 SLang_Array_Type *at;
826 unsigned long *s;
827 SLuindex_Type i;
828
829 if (-1 == SLang_pop_array_of_type (&at, SLANG_ULONG_TYPE))
830 return -1;
831
832 if (at->num_elements == 0)
833 {
834 SLang_verror (SL_InvalidParm_Error, "The seed array has no elements");
835 SLang_free_array (at);
836 return -1;
837 }
838
839 s = (unsigned long *)at->data;
840 i = 0;
841 while (i < NUM_SEEDS)
842 {
843 seeds[i] = *s;
844 i++;
845 if (i < at->num_elements)
846 s++;
847 }
848 SLang_free_array (at);
849 return 0;
850 }
851
generate_seeds(unsigned long seeds[NUM_SEEDS])852 static void generate_seeds (unsigned long seeds[NUM_SEEDS])
853 {
854 unsigned int i;
855 unsigned long s = (unsigned long) time(NULL)*(unsigned long) getpid ();
856 for (i = 0; i < NUM_SEEDS; i++)
857 {
858 s = s*69069UL + 1013904243UL;
859 seeds[i] = s;
860 }
861 }
862
new_rand_intrin(void)863 static void new_rand_intrin (void) /*{{{*/
864 {
865 unsigned long seeds[NUM_SEEDS];
866 Rand_Type *r;
867 SLang_MMT_Type *mmt;
868
869 if (SLang_Num_Function_Args == 1)
870 {
871 if (-1 == pop_seeds (seeds))
872 return;
873 }
874 else generate_seeds (seeds);
875
876 if (NULL == (r = create_random (seeds)))
877 return;
878
879 if (NULL == (mmt = SLang_create_mmt (Rand_Type_Id, (VOID_STAR) r)))
880 {
881 free_random (r);
882 return;
883 }
884
885 if (0 == SLang_push_mmt (mmt))
886 return;
887
888 SLang_free_mmt (mmt);
889 }
890
891 /*}}}*/
892
893 /*{{{ Utility Functions */
894
destroy_rand_type(SLtype type,VOID_STAR vr)895 static void destroy_rand_type (SLtype type, VOID_STAR vr)
896 {
897 (void) type;
898 free_random ((Rand_Type *) vr);
899 }
900
pop_rand_type_and_dims(int argc,SLang_MMT_Type ** mmtp,SLindex_Type * dims,unsigned int * ndims,int * is_scalarp)901 static int pop_rand_type_and_dims (int argc, SLang_MMT_Type **mmtp,
902 SLindex_Type *dims, unsigned int *ndims,
903 int *is_scalarp)
904 {
905 int type;
906 SLang_MMT_Type *mmt;
907
908 *mmtp = NULL;
909
910 switch (argc)
911 {
912 default:
913 SLang_verror (SL_NumArgs_Error, "Expecting 0, 1, or 2 arguments");
914 return -1;
915
916 case 0:
917 *is_scalarp = 1;
918 return 0;
919
920 case 1:
921 type = SLang_peek_at_stack ();
922 if (type == Rand_Type_Id)
923 {
924 if (NULL == (mmt = SLang_pop_mmt (Rand_Type_Id)))
925 return -1;
926 *is_scalarp = 1;
927 *mmtp = mmt;
928 return 0;
929 }
930 break;
931
932 case 2:
933 type = SLang_peek_at_stack ();
934 break;
935 }
936
937 *is_scalarp = 0;
938
939 if (type != SLANG_ARRAY_TYPE)
940 {
941 if (-1 == SLang_pop_array_index (dims))
942 return -1;
943
944 *ndims = 1;
945 }
946 else
947 {
948 SLang_Array_Type *at;
949 unsigned int i, imax;
950
951 if (-1 == SLang_pop_array (&at, 1))
952 return -1;
953
954 *ndims = imax = at->num_dims;
955 for (i = 0; i < imax; i++)
956 dims[i] = at->dims[i];
957 SLang_free_array (at);
958 }
959
960 if (argc == 2)
961 {
962 if (NULL == (mmt = SLang_pop_mmt (Rand_Type_Id)))
963 return -1;
964
965 *mmtp = mmt;
966 }
967 return 0;
968 }
969
do_xxxrand(int argc,SLtype type,void (* func)(Rand_Type *,VOID_STAR,SLuindex_Type,VOID_STAR),VOID_STAR parms,int * is_scalar_p,VOID_STAR scalar_addr)970 static int do_xxxrand (int argc, SLtype type,
971 void (*func)(Rand_Type *, VOID_STAR, SLuindex_Type, VOID_STAR),
972 VOID_STAR parms,
973 int *is_scalar_p, VOID_STAR scalar_addr)
974 {
975 SLang_Array_Type *at;
976 SLindex_Type dims [SLARRAY_MAX_DIMS];
977 unsigned int ndims;
978 SLang_MMT_Type *mmt;
979 Rand_Type *rt;
980 int is_scalar;
981 int status = -1;
982
983 if (-1 == pop_rand_type_and_dims (argc, &mmt, dims, &ndims, &is_scalar))
984 return -1;
985
986 if (mmt != NULL)
987 {
988 if (NULL == (rt = (Rand_Type *) SLang_object_from_mmt (mmt)))
989 goto free_return;
990 }
991 else
992 rt = Default_Rand;
993
994 *is_scalar_p = is_scalar;
995
996 if (is_scalar)
997 {
998 (*func) (rt, scalar_addr, 1, parms);
999 status = 0;
1000 goto free_return;
1001 }
1002
1003 if (NULL == (at = SLang_create_array (type, 0, NULL, dims, ndims)))
1004 goto free_return;
1005
1006 (*func) (rt, at->data, at->num_elements, parms);
1007 status = SLang_push_array (at, 0);
1008 SLang_free_array (at);
1009 /* drop */
1010
1011 free_return:
1012 if (mmt != NULL)
1013 SLang_free_mmt (mmt);
1014
1015 return status;
1016 }
1017
1018 /* The calling syntax for the generators is:
1019 * rand_foo ([Rand_Type,] args... [,num]);
1020 */
check_stack_args(int num_args,int num_parms,SLFUTURE_CONST char * usage,int * nargsp)1021 static int check_stack_args (int num_args, int num_parms, SLFUTURE_CONST char *usage, int *nargsp)
1022 {
1023 if ((num_args < num_parms) || (num_args > num_parms + 2))
1024 goto usage_error;
1025
1026 *nargsp = num_args - num_parms;
1027
1028 if ((num_args == num_parms) || (num_parms == 0))
1029 return 0; /* rand_foo (parms...) */
1030
1031 if (num_args == num_parms + 2)
1032 {
1033 if (Rand_Type_Id != SLang_peek_at_stack_n (num_args-1))
1034 goto usage_error;
1035
1036 /* rand_foo (r, parms..., num) */
1037 }
1038 else if (Rand_Type_Id == SLang_peek_at_stack_n (num_args-1))
1039 return 0; /* rand_foo (r, parms...) */
1040
1041 return SLroll_stack (num_parms + 1);
1042
1043 usage_error:
1044 SLang_verror (SL_Usage_Error, "Usage: %s", usage);
1045 return -1;
1046 }
1047
1048 /*}}}*/
1049
urand_intrin(void)1050 static void urand_intrin (void) /*{{{*/
1051 {
1052 SLFUTURE_CONST char *usage = "r = rand_uniform ([Rand_Type] [num])";
1053 int is_scalar;
1054 int nargs;
1055 double d;
1056
1057 if (-1 == check_stack_args (SLang_Num_Function_Args, 0, usage, &nargs))
1058 return;
1059
1060 if (-1 == do_xxxrand (nargs, SLANG_DOUBLE_TYPE,
1061 generate_random_doubles, NULL, &is_scalar, &d))
1062 return;
1063
1064 if (is_scalar)
1065 (void) SLang_push_double (d);
1066 }
1067
1068 /*}}}*/
1069
urand_pos_intrin(void)1070 static void urand_pos_intrin (void) /*{{{*/
1071 {
1072 SLFUTURE_CONST char *usage = "r = rand_uniform_pos ([Rand_Type] [num])";
1073 int is_scalar;
1074 int nargs;
1075 double d;
1076
1077 if (-1 == check_stack_args (SLang_Num_Function_Args, 0, usage, &nargs))
1078 return;
1079
1080 if (-1 == do_xxxrand (nargs, SLANG_DOUBLE_TYPE,
1081 generate_random_open_doubles, NULL, &is_scalar, &d))
1082 return;
1083
1084 if (is_scalar)
1085 (void) SLang_push_double (d);
1086 }
1087
1088 /*}}}*/
1089
rand_intrin(void)1090 static void rand_intrin (void) /*{{{*/
1091 {
1092 SLFUTURE_CONST char *usage = "r = rand ([Rand_Type] [num])";
1093 int is_scalar;
1094 int nargs;
1095 uint32 u;
1096
1097 if (-1 == check_stack_args (SLang_Num_Function_Args, 0, usage, &nargs))
1098 return;
1099
1100 if (-1 == do_xxxrand (nargs, UINT32_TYPE,
1101 generate_random_uints, NULL, &is_scalar, &u))
1102 return;
1103 if (is_scalar)
1104 (void) PUSH_UINT32 (u);
1105 }
1106
1107 /*}}}*/
1108
rand_gauss_intrin(void)1109 static void rand_gauss_intrin (void) /*{{{*/
1110 {
1111 SLFUTURE_CONST char *usage = "r = rand_gauss ([Rand_Type,] sigma [,num])";
1112 int is_scalar;
1113 double d, sigma;
1114 int nargs;
1115
1116 if (-1 == check_stack_args (SLang_Num_Function_Args, 1, usage, &nargs))
1117 return;
1118
1119 if (-1 == SLang_pop_double (&sigma))
1120 return;
1121 sigma = fabs(sigma);
1122
1123 if (-1 == do_xxxrand (nargs, SLANG_DOUBLE_TYPE,
1124 generate_gaussian_randoms, (VOID_STAR) &sigma,
1125 &is_scalar, &d))
1126 return;
1127 if (is_scalar)
1128 (void) SLang_push_double (d);
1129 }
1130
1131 /*}}}*/
1132
rand_beta_intrin(void)1133 static void rand_beta_intrin (void) /*{{{*/
1134 {
1135 SLFUTURE_CONST char *usage = "r = rand_beta ([Rand_Type,] a, b [,num])";
1136 int is_scalar;
1137 double d;
1138 double parms[2];
1139 int nargs;
1140
1141 if (-1 == check_stack_args (SLang_Num_Function_Args, 2, usage, &nargs))
1142 return;
1143
1144 if ((-1 == SLang_pop_double (parms+1))
1145 || (-1 == SLang_pop_double (parms)))
1146 return;
1147
1148 if ((parms[0] <= 0.0) || (parms[1] <= 0.0))
1149 {
1150 SLang_verror (SL_Domain_Error, "rand_beta parameters must be > 0");
1151 return;
1152 }
1153
1154 if (-1 == do_xxxrand (nargs, SLANG_DOUBLE_TYPE,
1155 generate_beta_randoms, (VOID_STAR) parms,
1156 &is_scalar, &d))
1157 return;
1158 if (is_scalar)
1159 (void) SLang_push_double (d);
1160 }
1161
1162 /*}}}*/
1163
rand_cauchy_intrin(void)1164 static void rand_cauchy_intrin (void) /*{{{*/
1165 {
1166 SLFUTURE_CONST char *usage = "r = rand_cauchy ([Rand_Type,] gamma, [,num])";
1167 int is_scalar;
1168 double d;
1169 double a;
1170 int nargs;
1171
1172 if (-1 == check_stack_args (SLang_Num_Function_Args, 1, usage, &nargs))
1173 return;
1174
1175 if (-1 == SLang_pop_double (&a))
1176 return;
1177
1178 a = fabs(a);
1179
1180 if (-1 == do_xxxrand (nargs, SLANG_DOUBLE_TYPE,
1181 generate_cauchy_randoms, (VOID_STAR) &a,
1182 &is_scalar, &d))
1183 return;
1184 if (is_scalar)
1185 (void) SLang_push_double (d);
1186 }
1187
1188 /*}}}*/
1189
rand_geometric_intrin(void)1190 static void rand_geometric_intrin (void) /*{{{*/
1191 {
1192 SLFUTURE_CONST char *usage = "r = rand_geometric ([Rand_Type,] p, [,num])";
1193 int is_scalar;
1194 double p;
1195 unsigned int d;
1196 int nargs;
1197
1198 if (-1 == check_stack_args (SLang_Num_Function_Args, 1, usage, &nargs))
1199 return;
1200
1201 if (-1 == SLang_pop_double (&p))
1202 return;
1203
1204 if ((p < 0.0) || (p > 1.0))
1205 {
1206 SLang_verror (SL_Domain_Error, "rand_geometric parameter must be beteen 0 and 1");
1207 return;
1208 }
1209
1210 if (-1 == do_xxxrand (nargs, SLANG_UINT_TYPE,
1211 generate_geometric_randoms, (VOID_STAR) &p,
1212 &is_scalar, &d))
1213 return;
1214 if (is_scalar)
1215 (void) SLang_push_uint (d);
1216 }
1217
1218 /*}}}*/
1219
rand_poisson_intrin(void)1220 static void rand_poisson_intrin (void) /*{{{*/
1221 {
1222 SLFUTURE_CONST char *usage = "r = rand_poisson ([Rand_Type,] mu [,num])";
1223 unsigned int p;
1224 int is_scalar;
1225 int nargs;
1226 double mu;
1227
1228 if (-1 == check_stack_args (SLang_Num_Function_Args, 1, usage, &nargs))
1229 return;
1230
1231 if (-1 == SLang_pop_double (&mu))
1232 return;
1233
1234 if (mu < 0.0)
1235 SLang_verror (SL_InvalidParm_Error, "The poisson rate must be non-negative");
1236
1237 if (-1 == do_xxxrand (nargs, SLANG_UINT_TYPE,
1238 generate_poisson_randoms, (VOID_STAR) &mu,
1239 &is_scalar, &p))
1240 return;
1241
1242 if (is_scalar)
1243 (void) SLang_push_uint (p);
1244 }
1245
1246 /*}}}*/
1247
rand_gamma_intrin(void)1248 static void rand_gamma_intrin (void) /*{{{*/
1249 {
1250 SLFUTURE_CONST char *usage = "r = rand_gamma([Rand_Type,] k, theta [,num])";
1251 int is_scalar, nargs;
1252 double parms[2];
1253 double p, k, theta;
1254
1255 if (-1 == check_stack_args (SLang_Num_Function_Args, 2, usage, &nargs))
1256 return;
1257
1258 if ((-1 == SLang_pop_double (&theta))
1259 || (-1 == SLang_pop_double (&k)))
1260 return;
1261
1262 if ((theta <= 0) || (k <= 0))
1263 {
1264 SLang_verror (SL_InvalidParm_Error, "rand_gamma assumes k,theta>0");
1265 return;
1266 }
1267 parms[0] = k;
1268 parms[1] = theta;
1269
1270 if (-1 == do_xxxrand (nargs, SLANG_DOUBLE_TYPE,
1271 generate_gamma_randoms, (VOID_STAR) parms,
1272 &is_scalar, &p))
1273 return;
1274
1275 if (is_scalar)
1276 (void) SLang_push_double (p);
1277 }
1278
1279 /*}}}*/
1280
rand_binomial_intrin(void)1281 static void rand_binomial_intrin (void) /*{{{*/
1282 {
1283 SLFUTURE_CONST char *usage = "r = rand_binomial ([Rand_Type,] p, n [,num])";
1284 int is_scalar, nargs;
1285 Binomial_Parms_Type bp;
1286 unsigned int u;
1287 int n;
1288
1289 if (-1 == check_stack_args (SLang_Num_Function_Args, 2, usage, &nargs))
1290 return;
1291
1292 if ((-1 == SLang_pop_int (&n))
1293 || (-1 == SLang_pop_double (&bp.p)))
1294 return;
1295
1296 if ((n < 0) || (bp.p < 0.0) || (bp.p > 1.0))
1297 {
1298 SLang_verror (SL_InvalidParm_Error, "rand_binomial assumes 0<=p<=1 and n>=0");
1299 return;
1300 }
1301 bp.n = (unsigned int)n;
1302
1303 if (-1 == do_xxxrand (nargs, SLANG_UINT_TYPE,
1304 generate_binomial_randoms, (VOID_STAR) &bp,
1305 &is_scalar, &u))
1306 return;
1307
1308 if (is_scalar)
1309 (void) SLang_push_uint (u);
1310 }
1311
1312 /*}}}*/
1313
srand_intrin(void)1314 static void srand_intrin (void) /*{{{*/
1315 {
1316 SLang_MMT_Type *mmt = NULL;
1317 Rand_Type *r = Default_Rand;
1318 unsigned long seeds[NUM_SEEDS];
1319 int nargs = SLang_Num_Function_Args;
1320
1321 if (-1 == pop_seeds (seeds))
1322 return;
1323
1324 if (nargs == 2)
1325 {
1326 if (NULL == (mmt = SLang_pop_mmt (Rand_Type_Id)))
1327 return;
1328 r = (Rand_Type *)SLang_object_from_mmt (mmt);
1329 }
1330
1331 if (r != NULL)
1332 seed_random (r, seeds);
1333
1334 if (mmt != NULL)
1335 SLang_free_mmt (mmt);
1336 }
1337
1338 /*}}}*/
1339
rand_permutation_intrin(void)1340 static void rand_permutation_intrin (void)
1341 {
1342 int nargs = SLang_Num_Function_Args;
1343 Rand_Type *rt = Default_Rand;
1344 SLang_MMT_Type *mmt = NULL;
1345 SLang_Array_Type *at = NULL;
1346 int *data;
1347 SLindex_Type i, n;
1348
1349 switch (nargs)
1350 {
1351 default:
1352 SLang_verror (SL_Usage_Error, "Usage: p = rand_permutation([Rand_Type,], n)");
1353 return;
1354
1355 case 2:
1356 case 1:
1357 if (-1 == SLang_pop_array_index (&n))
1358 return;
1359 if (nargs == 2)
1360 {
1361 if (NULL == (mmt = SLang_pop_mmt (Rand_Type_Id)))
1362 return;
1363
1364 if (NULL == (rt = (Rand_Type *) SLang_object_from_mmt (mmt)))
1365 goto free_return;
1366 }
1367 }
1368
1369 if (n < 0)
1370 {
1371 SLang_verror (SL_InvalidParm_Error, "rand_permutation: expected n>=0");
1372 goto free_return;
1373 }
1374
1375 if (NULL == (at = SLang_create_array (SLANG_INT_TYPE, 0, NULL, &n, 1)))
1376 goto free_return;
1377
1378 data = (int *) at->data;
1379 for (i = 0; i < n; i++)
1380 data[i] = i;
1381
1382 /* Fisher-Yates */
1383 while (n > 1)
1384 {
1385 int k, p;
1386
1387 k = (int) (n*uniform_random (rt)); /* 0 <= k < n */
1388 n--;
1389 p = data[n];
1390 data[n] = data[k];
1391 data[k] = p;
1392 }
1393
1394 (void) SLang_push_array (at, 0);
1395
1396 free_return:
1397 if (at != NULL)
1398 SLang_free_array (at);
1399 if (mmt != NULL)
1400 SLang_free_mmt (mmt);
1401 }
1402
1403 static SLang_Intrin_Fun_Type Module_Intrinsics [] =
1404 {
1405 MAKE_INTRINSIC_0("rand", rand_intrin, SLANG_VOID_TYPE),
1406 MAKE_INTRINSIC_0("srand", srand_intrin, SLANG_VOID_TYPE),
1407 MAKE_INTRINSIC_0("rand_uniform", urand_intrin, SLANG_VOID_TYPE),
1408 MAKE_INTRINSIC_0("rand_uniform_pos", urand_pos_intrin, SLANG_VOID_TYPE),
1409 MAKE_INTRINSIC_0("rand_gauss", rand_gauss_intrin, SLANG_VOID_TYPE),
1410 MAKE_INTRINSIC_0("rand_poisson", rand_poisson_intrin, SLANG_VOID_TYPE),
1411 MAKE_INTRINSIC_0("rand_gamma", rand_gamma_intrin, SLANG_VOID_TYPE),
1412 MAKE_INTRINSIC_0("rand_binomial", rand_binomial_intrin, SLANG_VOID_TYPE),
1413 MAKE_INTRINSIC_0("rand_beta", rand_beta_intrin, SLANG_VOID_TYPE),
1414 MAKE_INTRINSIC_0("rand_cauchy", rand_cauchy_intrin, SLANG_VOID_TYPE),
1415 MAKE_INTRINSIC_0("rand_geometric", rand_geometric_intrin, SLANG_VOID_TYPE),
1416
1417 MAKE_INTRINSIC_0("rand_permutation", rand_permutation_intrin, SLANG_VOID_TYPE),
1418
1419 MAKE_INTRINSIC_0("rand_new", new_rand_intrin, SLANG_VOID_TYPE),
1420 SLANG_END_INTRIN_FUN_TABLE
1421 };
1422
init_rand_module_ns(char * ns_name)1423 int init_rand_module_ns (char *ns_name)
1424 {
1425 SLang_NameSpace_Type *ns = SLns_create_namespace (ns_name);
1426 if (ns == NULL)
1427 return -1;
1428
1429 if (Default_Rand == NULL)
1430 {
1431 unsigned long seeds[NUM_SEEDS];
1432 generate_seeds (seeds);
1433 Default_Rand = create_random (seeds);
1434 if (Default_Rand == NULL)
1435 return -1;
1436
1437 init_poisson ();
1438 }
1439
1440 if (Rand_Type_Id == -1)
1441 {
1442 SLang_Class_Type *cl;
1443
1444 if (NULL == (cl = SLclass_allocate_class ("Rand_Type")))
1445 return -1;
1446
1447 (void) SLclass_set_destroy_function (cl, destroy_rand_type);
1448
1449 if (-1 == SLclass_register_class (cl, SLANG_VOID_TYPE,
1450 sizeof (Rand_Type),
1451 SLANG_CLASS_TYPE_MMT))
1452 return -1;
1453
1454 Rand_Type_Id = SLclass_get_class_id (cl);
1455 }
1456
1457 if (-1 == SLns_add_intrin_fun_table (ns, Module_Intrinsics, NULL))
1458 return -1;
1459
1460 return 0;
1461 }
1462
1463 /* This function is optional */
deinit_rand_module(void)1464 void deinit_rand_module (void)
1465 {
1466 }
1467