File Coverage

blib/lib/Math/Random.pm
Criterion Covered Total %
statement 204 244 83.6
branch 59 146 40.4
condition 11 58 18.9
subroutine 21 28 75.0
pod 21 23 91.3
total 316 499 63.3


line stmt bran cond sub pod time code
1             package Math::Random;
2              
3 1     1   792 use strict;
  1         2  
  1         23  
4 1     1   4 use Carp;
  1         2  
  1         78  
5 1     1   4 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $AUTOLOAD);
  1         2  
  1         3090  
6              
7             require Exporter;
8             require DynaLoader;
9             require AutoLoader;
10              
11             @ISA = qw(Exporter DynaLoader);
12             $VERSION = '0.72';
13              
14             @EXPORT = qw(random_normal
15             random_permutation
16             random_permuted_index
17             random_uniform
18             random_uniform_integer
19             random_seed_from_phrase
20             random_get_seed
21             random_set_seed_from_phrase
22             random_set_seed
23             );
24              
25             @EXPORT_OK = qw(random_beta
26             random_chi_square
27             random_exponential
28             random_f
29             random_gamma
30             random_multivariate_normal
31             random_multinomial
32             random_noncentral_chi_square
33             random_noncentral_f
34             random_normal
35             random_permutation
36             random_permuted_index
37             random_uniform
38             random_poisson
39             random_uniform_integer
40             random_negative_binomial
41             random_binomial
42             random_seed_from_phrase
43             random_get_seed
44             random_set_seed_from_phrase
45             random_set_seed
46             );
47              
48             %EXPORT_TAGS = ( all => [ @EXPORT_OK ] );
49              
50             sub AUTOLOAD {
51             # This AUTOLOAD is used to 'autoload' constants from the constant()
52             # XS function. If a constant is not found then control is passed
53             # to the AUTOLOAD in AutoLoader.
54              
55 0     0   0 my $constname;
56 0         0 ($constname = $AUTOLOAD) =~ s/.*:://;
57 0 0       0 croak "& not defined" if $constname eq 'constant';
58 0 0       0 my $val = constant($constname, @_ ? $_[0] : 0);
59 0 0       0 if ($! != 0) {
60 0 0       0 if ($! =~ /Invalid/) {
61 0         0 $AutoLoader::AUTOLOAD = $AUTOLOAD;
62 0         0 goto &AutoLoader::AUTOLOAD;
63             }
64             else {
65 0         0 croak "Your vendor has not defined Math::Random macro $constname";
66             }
67             }
68 0     0   0 *$AUTOLOAD = sub () { $val };
  0         0  
69 0         0 goto &$AUTOLOAD;
70             }
71              
72             bootstrap Math::Random $VERSION;
73              
74              
75             ### set seeds by default
76             salfph(get_seed() || scalar(localtime));
77              
78             #####################################################################
79             # RANDOM DEVIATE GENERATORS #
80             #####################################################################
81              
82             sub random_beta { # Arguments: ($n,$aa,$bb)
83 1 50   1 1 43 croak "Usage: random_beta(\$n,\$aa,\$bb)" if scalar(@_) < 3;
84 1         2 my($n, $aa, $bb) = @_;
85 1 50 33     8 croak("($aa = \$aa < 1.0E-37) or ($bb = \$bb < 1.0E-37)\nin ".
86             "random_beta(\$n,\$aa,\$bb)")
87             if (($aa < 1.0E-37) or ($bb < 1.0E-37));
88 1 50       2 return genbet($aa,$bb) unless wantarray();
89 1         1 my $val;
90 1         3 my @ans = (0) x $n;
91 1         7 foreach $val (@ans) { $val = genbet($aa,$bb); }
  3         19  
92 1         4 return @ans;
93             }
94              
95             sub random_chi_square { # Arguments: ($n,$df)
96 1 50   1 1 50 croak "Usage: random_chi_square(\$n,\$df)" if scalar(@_) < 2;
97 1         6 my($n, $df) = @_;
98 1 50       3 croak "$df = \$df <= 0\nin random_chi_square(\$n,\$df)" if ($df <= 0);
99 1 50       3 return genchi($df) unless wantarray();
100 1         1 my $val;
101 1         3 my @ans = (0) x $n;
102 1         2 foreach $val (@ans) { $val = genchi($df); }
  3         14  
103 1         4 return @ans;
104             }
105              
106             sub random_exponential { # Arguments: ($n,$av), defaults (1,1)
107 1 0   1 1 45 return wantarray() ? (genexp(1)) : genexp(1)
    50          
108             if scalar(@_) == 0; # default behavior if no arguments
109 1         2 my($n, $av) = @_;
110 1 50       2 $av = 1 unless defined($av); # default $av is 1
111 1 50       3 croak "$av = \$av < 0\nin random_exponential(\$n,\$av)" if ($av < 0);
112 1 50       2 return genexp($av) unless wantarray();
113 1         2 my $val;
114 1         2 my @ans = (0) x $n;
115 1         2 foreach $val (@ans) { $val = genexp($av); }
  3         9  
116 1         4 return @ans;
117             }
118              
119             sub random_f { # Arguments: ($n,$dfn,$dfd)
120 1 50   1 1 44 croak "Usage: random_f(\$n,\$dfn,\$dfd)" if scalar(@_) < 3;
121 1         2 my($n, $dfn, $dfd) = @_;
122 1 50 33     7 croak("($dfn = \$dfn <= 0) or ($dfd = \$dfd <= 0)\nin ".
123             "random_f(\$n,\$dfn,\$dfd)") if (($dfn <= 0) or ($dfd <= 0));
124 1 50       4 return genf($dfn,$dfd) unless wantarray();
125 1         1 my $val;
126 1         3 my @ans = (0) x $n;
127 1         2 foreach $val (@ans) { $val = genf($dfn,$dfd); }
  3         24  
128 1         4 return @ans;
129             }
130              
131             sub random_gamma { # Arguments: ($n,$a,$r)
132 1 50   1 1 44 croak "Usage: random_gamma(\$n,\$a,\$r)" if scalar(@_) < 3;
133 1         2 my($n, $a, $r) = @_;
134 1 50 33     6 croak "($a = \$a <= 0) or ($r = \$r <= 0)\nin random_gamma(\$n,\$a,\$r)"
135             if (($a <= 0) or ($r <= 0));
136 1 50       3 return gengam($a,$r) unless wantarray();
137 1         7 my $val;
138 1         3 my @ans = (0) x $n;
139 1         2 foreach $val (@ans) { $val = gengam($a,$r); }
  3         10  
140 1         3 return @ans;
141             }
142              
143             sub random_multivariate_normal { # Arguments: ($n, @mean, @covar(2-dim'l))
144              
145 1 50   1 1 46 croak "Usage: random_multivariate_normal(\$n,\@mean,\@covar(2-dim'l))"
146             if (scalar(@_)) < 3;
147 1         2 my $n = shift(@_); # first element is number of obs. desired
148 1         3 my $p = scalar(@_)/2; # best guess at dimension of deviate
149            
150             # check outline of arguments
151 1 50 33     22 croak("Sizes of \@mean and \@covar don't match\nin ".
      33        
152             "random_multivariate_normal(\$n, \@mean, \@covar(2-dim'l))")
153             unless (($p == int($p)) and ("$_[$p - 1]" !~ /^ARRAY/) and
154             ("$_[$p]" =~ /^ARRAY/));
155            
156             # linearize input - it seems faster to push
157 1         2 my @linear = ();
158            
159 1         3 push @linear, splice(@_, 0, $p); # fill first $p slots w/ mean
160            
161             # expand array references
162 1         1 my $ref;
163 1         26 foreach $ref (@_) { # for the rest of the input
164            
165             # check length of row of @covariance
166             croak("\@covar is not a $p by $p array ($p is size of \@mean)\nin ".
167             "random_multivariate_normal(\$n, \@mean, \@covar(2-dim'l))")
168 2 50       3 unless (scalar(@{$ref}) == $p);
  2         10  
169            
170 2         2 push @linear, @{$ref};
  2         6  
171             }
172            
173             # load float working array with linearized input
174 1 50       6 putflt(@linear) or
175             croak "Unable to allocate memory\nin random_multivariate_normal";
176            
177             # initialize parameter array for multivariate normal generator
178 1 50       8 psetmn($p) or
179             croak "Unable to allocate memory\nin random_multivariate_normal";
180            
181 1 50       3 unless (wantarray()) {
182             ### if called in a scalar context, returns single refernce to obs
183 0         0 pgenmn();
184 0         0 return [ getflt($p) ];
185             }
186            
187             # otherwise return an $n by $p array of obs.
188 1         3 my @ans = (0) x $n;
189 1         2 foreach $ref (@ans) {
190 2         7 pgenmn();
191 2         5 $ref = [ getflt($p) ];
192             }
193 1         5 return @ans;
194             }
195              
196             sub random_multinomial { # Arguments: ($n,@p)
197 1     1 1 41 my($n, @p) = @_;
198 1         2 my $ncat = scalar(@p); # number of categories
199 1         1 $n = int($n);
200 1 50       11 croak "$n = \$n < 0\nin random_multinomial(\$n,\@p)" if ($n < 0);
201 1 50       3 croak "$ncat = (length of \@p) < 2\nin random_multinomial(\$n,\@p)"
202             if ($ncat < 2);
203 1 50       12 rspriw($ncat) or croak "Unable to allocate memory\nin random_multinomial";
204 1         14 my($i,$sum,$val) = (0,0,0);
205 1         1 pop @p;
206 1 50       5 rsprfw(scalar(@p)) or
207             croak "Unable to allocate memory\nin random_multinomial";
208 1         3 foreach $val (@p) {
209 2 50 33     12 croak "$val = (some \$p[i]) < 0 or > 1\nin random_multinomial(\$n,\@p)"
210             if (($val < 0) or ($val > 1));
211 2         5 svprfw($i,$val);
212 2         2 $i++;
213 2         4 $sum += $val;
214             }
215 1 50       3 croak "Sum of \@p > 1\nin random_multinomial(\$n,\@p)" if ($sum > 0.99999);
216 1         6 pgnmul($n, $ncat);
217             ### get the results
218 1         2 $i = 0;
219 1         2 foreach $val (@p) {
220 2         5 $val = gvpriw($i);
221 2         3 $i++;
222             }
223 1         3 push @p, gvpriw($i);
224 1         4 return @p;
225             }
226              
227             sub random_noncentral_chi_square { # Arguments: ($n,$df,$nonc)
228 0 0   0 1 0 croak "Usage: random_noncentral_chi_square(\$n,\$df,\$nonc)"
229             if scalar(@_) < 3;
230 0         0 my($n, $df, $nonc) = @_;
231 0 0 0     0 croak("($df = \$df < 1) or ($nonc = \$nonc) < 0\n".
232             "in random_noncentral_chi_square(\$n,\$df,\$nonc)")
233             if (($df < 1) or ($nonc < 0));
234 0 0       0 return gennch($df,$nonc) unless wantarray();
235 0         0 my $val;
236 0         0 my @ans = (0) x $n;
237 0         0 foreach $val (@ans) { $val = gennch($df,$nonc); }
  0         0  
238 0         0 return @ans;
239             }
240              
241             sub random_noncentral_f { # Arguments: ($n,$dfn,$dfd,$nonc)
242 0 0   0 1 0 croak "Usage: random_noncentral_f(\$n,\$dfn,\$dfd,\$nonc)"
243             if scalar(@_) < 4;
244 0         0 my($n, $dfn, $dfd, $nonc) = @_;
245 0 0 0     0 croak("($dfn = \$dfn < 1) or ($dfd = \$dfd <= 0) or ($nonc ".
      0        
246             "= \$nonc < 0)\nin random_noncentral_f(\$n,\$dfn,\$dfd,\$nonc)")
247             if (($dfn < 1) or ($dfd <= 0) or ($nonc < 0));
248 0 0       0 return gennf($dfn,$dfd,$nonc) unless wantarray();
249 0         0 my $val;
250 0         0 my @ans = (0) x $n;
251 0         0 foreach $val (@ans) { $val = gennf($dfn,$dfd,$nonc); }
  0         0  
252 0         0 return @ans;
253             }
254              
255             sub random_normal { # Arguments: ($n,$av,$sd), defaults (1,0,1)
256 1 0   1 1 35 return wantarray() ? (gennor(0,1)) : gennor(0,1)
    50          
257             if scalar(@_) == 0; # default behavior if no arguments
258 1         2 my($n, $av, $sd) = @_;
259 1 50       4 $av = 0 unless defined($av); # $av defaults to 0
260 1 50       3 $sd = 1 unless defined($sd); # $sd defaults to 1, even if $av specified
261 1 50       3 croak "$sd = \$sd < 0\nin random_normal([\$n[,\$av[,\$sd]]])" if ($sd < 0);
262 1 50       8 return gennor($av,$sd) unless wantarray();
263 1         6 my $val;
264 1         3 my @ans = (0) x $n;
265 1         1 foreach $val (@ans) { $val = gennor($av,$sd); }
  3         11  
266 1         4 return @ans;
267             }
268              
269             sub random_permutation { # Argument: (@array) - array to be permuted.
270 1     1 1 49 my $n = scalar(@_); # number of elements to be permuted
271 1 50       3 return () if $n == 0;
272 1 50       9 rspriw($n) or
273             croak "Unable to allocate memory\nin random_permutation";
274 1         7 pgnprm($n);
275 1         2 my($val, $i) = (0,0);
276 1         3 my @ans = (0) x $n;
277 1         2 foreach $val (@ans) {
278 8         14 $val = gvpriw($i);
279 8         10 $i++;
280             }
281 1         34 return @_[@ans];
282             }
283              
284             sub random_permuted_index { # Argument: $n = scalar(@array) (for permutation)
285 1 50   1 1 32 croak "Usage: random_permuted_index(\$n)" if scalar(@_) < 1;
286 1         2 my $n = int(shift(@_)); # number of elements to be permuted
287 1 50       2 croak "$n = \$n < 0 in random_permuted_index(\$n)" if $n < 0;
288 1 50       3 return () if $n == 0;
289 1 50       4 rspriw($n) or
290             croak "Unable to allocate memory\nin random_permuted_index";
291 1         4 pgnprm($n);
292 1         2 my($val, $i) = (0,0);
293 1         3 my @ans = (0) x $n;
294 1         2 foreach $val (@ans) {
295 9         14 $val = gvpriw($i);
296 9         10 $i++;
297             }
298 1         5 return @ans;
299             }
300              
301             sub random_uniform { # Arguments: ($n,$low,$high), defaults (1,0,1)
302 1 0   1 1 28 return wantarray() ? (genunf(0,1)) : genunf(0,1)
    50          
303             if scalar(@_) == 0;
304 1 50       3 croak "Usage: random_uniform([\$n,[\$low,\$high]])"
305             if scalar(@_) == 2; # only default is (0,1) for ($low,$high) both undef
306 1         3 my($n, $low, $high) = @_;
307 1 50       3 $low = 0 unless defined($low); # default for $low is 0
308 1 50       3 $high = 1 unless defined($high); # default for $high is 1
309 1 50       3 croak("$low = \$low > \$high = $high\nin ".
310             "random_uniform([\$n,[\$low,\$high]])") if ($low > $high);
311 1 50       4 return genunf($low,$high) unless wantarray();
312 1         1 my $val;
313 1         5 my @ans = (0) x $n;
314 1         3 foreach $val (@ans) { $val = genunf($low,$high); }
  3         18  
315 1         7 return @ans;
316             }
317              
318             sub random_poisson { # Arguments: ($n, $mu)
319 1 50   1 1 46 croak "Usage: random_poisson(\$n,\$mu)" if scalar(@_) < 2;
320 1         2 my($n, $mu) = @_;
321 1 50       3 croak "$mu = \$mu < 0\nin random_poisson(\$n,\$mu)" if ($mu < 0);
322 1 50       3 return ignpoi($mu) unless wantarray();
323 1         1 my $val;
324 1         3 my @ans = (0) x $n;
325 1         2 foreach $val (@ans) { $val = ignpoi($mu); }
  3         14  
326 1         4 return @ans;
327             }
328              
329             sub random_uniform_integer { # Arguments: ($n,$low,$high)
330 1 50   1 1 68 croak "Usage: random_uniform_integer(\$n,\$low,\$high)" if scalar(@_) < 3;
331 1         3 my($n, $low, $high) = @_;
332 1         2 $low = int($low);
333 1         2 $high = int($high);
334 1 50       4 croak("$low = \$low > \$high = $high\nin ".
335             "random_uniform_integer(\$n,\$low,\$high)") if ($low > $high);
336 1         3 my $range = $high - $low;
337 1 50       3 croak("$range = (\$high - \$low) > 2147483561\nin ".
338             "random_uniform_integer(\$n,\$low,\$high)") if ($range > 2147483561);
339 1 50       2 return ($low + ignuin(0,$range)) unless wantarray();
340 1         2 my $val;
341 1         3 my @ans = (0) x $n;
342 1         2 foreach $val (@ans) { $val = $low + ignuin(0,$range); }
  3         8  
343 1         4 return @ans;
344             }
345              
346             sub random_negative_binomial { # Arguments: ($n,$ne,$p)
347 1 50   1 1 67 croak "Usage: random_negative_binomial(\$n,\$ne,\$p)" if scalar(@_) < 3;
348 1         3 my($n, $ne, $p) = @_;
349 1         3 $ne = int($ne);
350 1 50 33     10 croak("($ne = \$ne <= 0) or ($p = \$p <= 0 or >= 1)\nin ".
      33        
351             "random_negative_binomial(\$n,\$ne,\$p)")
352             if (($ne <= 0) or (($p <= 0) or ($p >= 1)));
353 1 50       9 return ignnbn($ne,$p) unless wantarray();
354 1         3 my $val;
355 1         4 my @ans = (0) x $n;
356 1         3 foreach $val (@ans) { $val = ignnbn($ne,$p); }
  3         13  
357 1         5 return @ans;
358             }
359              
360             sub random_binomial { # Arguments: ($n,$nt,$p)
361 1 50   1 1 42 croak "Usage: random_binomial(\$n,\$nt,\$p)" if scalar(@_) < 3;
362 1         2 my($n, $nt, $p) = @_;
363 1         2 $nt = int($nt);
364 1 50 33     14 croak("($nt = \$nt < 0) or ($p = \$p < 0 or > 1)\nin ".
      33        
365             "random_binomial(\$n,\$nt,\$p)")
366             if (($nt < 0) or (($p < 0) or ($p > 1)));
367 1 50       3 return ignbin($nt,$p) unless wantarray();
368 1         2 my $val;
369 1         2 my @ans = (0) x $n;
370 1         2 foreach $val (@ans) { $val = ignbin($nt,$p); }
  3         13  
371 1         4 return @ans;
372             }
373              
374             #####################################################################
375             # SEED HANDLER FUNCTIONS #
376             #####################################################################
377              
378             sub random_seed_from_phrase { # Argument $phrase
379 0     0 1 0 my $phrase = shift(@_);
380 0   0     0 $phrase ||= "";
381 0         0 return phrtsd($phrase);
382             }
383              
384             sub random_get_seed { # no argument
385 0     0 1 0 return getsd();
386             }
387              
388             sub random_set_seed_from_phrase { # Argument $phrase
389 1     1 1 35 my $phrase = shift(@_);
390 1   50     4 $phrase ||= "";
391 1         25 salfph($phrase);
392 1         2 return 1;
393             }
394              
395             sub random_set_seed { # Argument @seed
396 0     0 1 0 my($seed1,$seed2) = @_;
397 0 0 0     0 croak("Usage: random_set_seed(\@seed)\n\@seed[0,1] must be two integers ".
      0        
      0        
      0        
      0        
398             "in the range (1,1) to (2147483562,2147483398)\nand usually comes ".
399             "from a call to random_get_seed() ".
400             "or\nrandom_seed_from_phrase(\$phrase).")
401             unless (((($seed1 == int($seed1)) and ($seed2 == int($seed2))) and
402             (($seed1 > 0) and ($seed2 > 0))) and
403             (($seed1 < 2147483563) and ($seed2 < 2147483399)));
404 0         0 setall($seed1,$seed2);
405 0         0 return 1;
406             }
407              
408             #####################################################################
409             # HELPER ROUTINES #
410             # These use the C work arrays and are not intended for export #
411             # (Currently only used in random_multivariate_normal) #
412             #####################################################################
413              
414             sub getflt {
415 2     2 0 4 my $n = $_[0];
416 2         8 my $val;
417 2         3 my $i = 0;
418 2         4 my @junk = (0) x $n;
419 2         4 foreach $val (@junk) {
420 4         9 $val = gvprfw($i);
421 4         8 $i++;
422             }
423 2         7 return @junk;
424             }
425              
426             sub putflt {
427 1     1 0 3 my $n = scalar(@_);
428 1 50       4 rsprfw($n) or return 0;
429 1         1 my $val;
430 1         3 my $i = 0;
431 1         3 foreach $val (@_) { # load up floats
432 6         8 svprfw($i,$val);
433 6         9 $i++;
434             }
435 1         9 return 1;
436             }
437              
438             # Autoload methods go after =cut, and are processed by the autosplit program.
439              
440             1;
441              
442             __END__