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

package Math::Random;

use strict;
use Carp;
use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $AUTOLOAD);

require Exporter;
require DynaLoader;
require AutoLoader;

@ISA = qw(Exporter DynaLoader);
$VERSION = '0.72';

@EXPORT = qw(random_normal
	     random_permutation
	     random_permuted_index
	     random_uniform
	     random_uniform_integer
             random_seed_from_phrase
             random_get_seed
             random_set_seed_from_phrase
             random_set_seed
	     );

@EXPORT_OK = qw(random_beta
		random_chi_square
		random_exponential
		random_f
		random_gamma
		random_multivariate_normal
		random_multinomial
		random_noncentral_chi_square
		random_noncentral_f
		random_normal
		random_permutation
		random_permuted_index
		random_uniform
		random_poisson
		random_uniform_integer
		random_negative_binomial
		random_binomial
                random_seed_from_phrase
                random_get_seed
                random_set_seed_from_phrase
                random_set_seed
		);

%EXPORT_TAGS = ( all => [ @EXPORT_OK ] );

sub AUTOLOAD {
    # This AUTOLOAD is used to 'autoload' constants from the constant()
    # XS function.  If a constant is not found then control is passed
    # to the AUTOLOAD in AutoLoader.

    my $constname;
    ($constname = $AUTOLOAD) =~ s/.*:://;
    croak "& not defined" if $constname eq 'constant';
    my $val = constant($constname, @_ ? $_[0] : 0);
    if ($! != 0) {
	if ($! =~ /Invalid/) {
	    $AutoLoader::AUTOLOAD = $AUTOLOAD;
	    goto &AutoLoader::AUTOLOAD;
	}
	else {
		croak "Your vendor has not defined Math::Random macro $constname";
	}
    }
    *$AUTOLOAD = sub () { $val };
    goto &$AUTOLOAD;
}

bootstrap Math::Random $VERSION;


### set seeds by default
salfph(get_seed() || scalar(localtime));

#####################################################################
#		      RANDOM DEVIATE GENERATORS                     #
#####################################################################

sub random_beta { # Arguments: ($n,$aa,$bb)
    croak "Usage: random_beta(\$n,\$aa,\$bb)" if scalar(@_) < 3;
    my($n, $aa, $bb) = @_;
    croak("($aa = \$aa < 1.0E-37) or ($bb = \$bb < 1.0E-37)\nin ".
	  "random_beta(\$n,\$aa,\$bb)")
	if (($aa < 1.0E-37) or ($bb < 1.0E-37));
    return genbet($aa,$bb) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = genbet($aa,$bb); }
    return @ans;
}

sub random_chi_square { # Arguments: ($n,$df)
    croak "Usage: random_chi_square(\$n,\$df)" if scalar(@_) < 2;
    my($n, $df) = @_;
    croak "$df = \$df <= 0\nin random_chi_square(\$n,\$df)" if ($df <= 0);
    return genchi($df) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = genchi($df); }
    return @ans;
}

sub random_exponential { # Arguments: ($n,$av), defaults (1,1)
    return wantarray() ? (genexp(1)) : genexp(1)
	if scalar(@_) == 0; # default behavior if no arguments
    my($n, $av) = @_;
    $av = 1 unless defined($av); # default $av is 1
    croak "$av = \$av < 0\nin random_exponential(\$n,\$av)" if ($av < 0);
    return genexp($av) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = genexp($av); }
    return @ans;
}

sub random_f { # Arguments: ($n,$dfn,$dfd)
    croak "Usage: random_f(\$n,\$dfn,\$dfd)" if scalar(@_) < 3;
    my($n, $dfn, $dfd) = @_;
    croak("($dfn = \$dfn <= 0) or ($dfd = \$dfd <= 0)\nin ".
	  "random_f(\$n,\$dfn,\$dfd)") if (($dfn <= 0) or ($dfd <= 0));
    return genf($dfn,$dfd) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = genf($dfn,$dfd); }
    return @ans;
}

sub random_gamma { # Arguments: ($n,$a,$r)
    croak "Usage: random_gamma(\$n,\$a,\$r)" if scalar(@_) < 3;
    my($n, $a, $r) = @_;
    croak "($a = \$a <= 0) or ($r = \$r <= 0)\nin random_gamma(\$n,\$a,\$r)"
	if (($a <= 0) or ($r <= 0));
    return gengam($a,$r) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = gengam($a,$r); }
    return @ans;
}

sub random_multivariate_normal { # Arguments: ($n, @mean, @covar(2-dim'l))

    croak "Usage: random_multivariate_normal(\$n,\@mean,\@covar(2-dim'l))"
	if (scalar(@_)) < 3;
    my $n = shift(@_); # first element is number of obs. desired
    my $p = scalar(@_)/2; # best guess at dimension of deviate

    # check outline of arguments
    croak("Sizes of \@mean and \@covar don't match\nin ".
	  "random_multivariate_normal(\$n, \@mean, \@covar(2-dim'l))")
	unless (($p == int($p)) and ("$_[$p - 1]" !~ /^ARRAY/) and
		("$_[$p]" =~ /^ARRAY/));

    # linearize input - it seems faster to push
    my @linear = ();

    push @linear, splice(@_, 0, $p); # fill first $p slots w/ mean

    # expand array references
    my $ref;
    foreach $ref (@_) { # for the rest of the input

	# check length of row of @covariance
	croak("\@covar is not a $p by $p array ($p is size of \@mean)\nin ".
	      "random_multivariate_normal(\$n, \@mean, \@covar(2-dim'l))")
	    unless (scalar(@{$ref}) == $p);

	push @linear, @{$ref};
    }

    # load float working array with linearized input
    putflt(@linear) or
	croak "Unable to allocate memory\nin random_multivariate_normal";

    # initialize parameter array for multivariate normal generator
    psetmn($p) or
	croak "Unable to allocate memory\nin random_multivariate_normal";

    unless (wantarray()) {
	### if called in a scalar context, returns single refernce to obs
	pgenmn();
	return [ getflt($p) ];
    }

    # otherwise return an $n by $p array of obs.
    my @ans = (0) x $n;
    foreach $ref (@ans) {
	pgenmn();
	$ref = [ getflt($p) ];
    }
    return @ans;
}

sub random_multinomial { # Arguments: ($n,@p)
    my($n, @p) = @_;
    my $ncat = scalar(@p); # number of categories
    $n = int($n);
    croak "$n = \$n < 0\nin random_multinomial(\$n,\@p)" if ($n < 0);
    croak "$ncat = (length of \@p) < 2\nin random_multinomial(\$n,\@p)"
	if ($ncat < 2);
    rspriw($ncat) or croak "Unable to allocate memory\nin random_multinomial";
    my($i,$sum,$val) = (0,0,0);
    pop @p;
    rsprfw(scalar(@p)) or
	croak "Unable to allocate memory\nin random_multinomial";
    foreach $val (@p) {
	croak "$val = (some \$p[i]) < 0 or > 1\nin random_multinomial(\$n,\@p)"
	    if (($val < 0) or ($val > 1));
	svprfw($i,$val);
	$i++;
	$sum += $val;
    }
    croak "Sum of \@p > 1\nin random_multinomial(\$n,\@p)" if ($sum > 0.99999);
    pgnmul($n, $ncat);
    ### get the results
    $i = 0;
    foreach $val (@p) {
	$val = gvpriw($i);
	$i++;
    }
    push @p, gvpriw($i);
    return @p;
}

sub random_noncentral_chi_square { # Arguments: ($n,$df,$nonc)
    croak "Usage: random_noncentral_chi_square(\$n,\$df,\$nonc)"
	if scalar(@_) < 3;
    my($n, $df, $nonc) = @_;
    croak("($df = \$df < 1) or ($nonc = \$nonc) < 0\n".
	  "in random_noncentral_chi_square(\$n,\$df,\$nonc)")
	if (($df < 1) or ($nonc < 0));
    return gennch($df,$nonc) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = gennch($df,$nonc); }
    return @ans;
}

sub random_noncentral_f { # Arguments: ($n,$dfn,$dfd,$nonc)
    croak "Usage: random_noncentral_f(\$n,\$dfn,\$dfd,\$nonc)"
	if scalar(@_) < 4;
    my($n, $dfn, $dfd, $nonc) = @_;
    croak("($dfn = \$dfn < 1) or ($dfd = \$dfd <= 0) or ($nonc ".
	  "= \$nonc < 0)\nin random_noncentral_f(\$n,\$dfn,\$dfd,\$nonc)")
	if (($dfn < 1) or ($dfd <= 0) or ($nonc < 0));
    return gennf($dfn,$dfd,$nonc) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = gennf($dfn,$dfd,$nonc); }
    return @ans;
}

sub random_normal { # Arguments: ($n,$av,$sd), defaults (1,0,1)
    return wantarray() ? (gennor(0,1)) : gennor(0,1)
	if scalar(@_) == 0; # default behavior if no arguments
    my($n, $av, $sd) = @_;
    $av = 0 unless defined($av); # $av defaults to 0
    $sd = 1 unless defined($sd); # $sd defaults to 1, even if $av specified
    croak "$sd = \$sd < 0\nin random_normal([\$n[,\$av[,\$sd]]])" if ($sd < 0);
    return gennor($av,$sd) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = gennor($av,$sd); }
    return @ans;
}

sub random_permutation { # Argument: (@array) - array to be permuted.
    my $n = scalar(@_); # number of elements to be permuted
    return () if $n == 0;
    rspriw($n) or
	croak "Unable to allocate memory\nin random_permutation";
    pgnprm($n);
    my($val, $i) = (0,0);
    my @ans = (0) x $n;
    foreach $val (@ans) {
	$val = gvpriw($i);
	$i++;
    }
    return @_[@ans];
}

sub random_permuted_index { # Argument: $n = scalar(@array) (for permutation)
    croak "Usage: random_permuted_index(\$n)" if scalar(@_) < 1;
    my $n = int(shift(@_)); # number of elements to be permuted
    croak "$n = \$n < 0 in random_permuted_index(\$n)" if $n < 0;
    return () if $n == 0;
    rspriw($n) or
	croak "Unable to allocate memory\nin random_permuted_index";
    pgnprm($n);
    my($val, $i) = (0,0);
    my @ans = (0) x $n;
    foreach $val (@ans) {
	$val = gvpriw($i);
	$i++;
    }
    return @ans;
}

sub random_uniform { # Arguments: ($n,$low,$high), defaults (1,0,1)
    return wantarray() ? (genunf(0,1)) : genunf(0,1)
	if scalar(@_) == 0;
    croak "Usage: random_uniform([\$n,[\$low,\$high]])"
	if scalar(@_) == 2; # only default is (0,1) for ($low,$high) both undef
    my($n, $low, $high) = @_;
    $low  = 0 unless defined($low); # default for $low is 0
    $high = 1 unless defined($high); # default for $high is 1
    croak("$low = \$low > \$high = $high\nin ".
	  "random_uniform([\$n,[\$low,\$high]])") if ($low > $high);
    return genunf($low,$high) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = genunf($low,$high); }
    return @ans;
}

sub random_poisson { # Arguments: ($n, $mu)
    croak "Usage: random_poisson(\$n,\$mu)" if scalar(@_) < 2;
    my($n, $mu) = @_;
    croak "$mu = \$mu < 0\nin random_poisson(\$n,\$mu)" if ($mu < 0);
    return ignpoi($mu) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = ignpoi($mu); }
    return @ans;
}

sub random_uniform_integer { # Arguments: ($n,$low,$high)
    croak "Usage: random_uniform_integer(\$n,\$low,\$high)" if scalar(@_) < 3;
    my($n, $low, $high) = @_;
    $low = int($low);
    $high = int($high);
    croak("$low = \$low > \$high = $high\nin ".
	  "random_uniform_integer(\$n,\$low,\$high)") if ($low > $high);
    my $range = $high - $low;
    croak("$range = (\$high - \$low) > 2147483561\nin ".
	  "random_uniform_integer(\$n,\$low,\$high)") if ($range > 2147483561);
    return ($low + ignuin(0,$range)) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = $low + ignuin(0,$range); }
    return @ans;
}

sub random_negative_binomial { # Arguments: ($n,$ne,$p)
    croak "Usage: random_negative_binomial(\$n,\$ne,\$p)" if scalar(@_) < 3;
    my($n, $ne, $p) = @_;
    $ne = int($ne);
    croak("($ne = \$ne <= 0) or ($p = \$p <= 0 or >= 1)\nin ".
	  "random_negative_binomial(\$n,\$ne,\$p)")
	if (($ne <= 0) or (($p <= 0) or ($p >= 1)));
    return ignnbn($ne,$p) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = ignnbn($ne,$p); }
    return @ans;
}

sub random_binomial { # Arguments: ($n,$nt,$p)
    croak "Usage: random_binomial(\$n,\$nt,\$p)" if scalar(@_) < 3;
    my($n, $nt, $p) = @_;
    $nt = int($nt);
    croak("($nt = \$nt < 0) or ($p = \$p < 0 or > 1)\nin ".
	  "random_binomial(\$n,\$nt,\$p)")
	if (($nt < 0) or (($p < 0) or ($p > 1)));
    return ignbin($nt,$p) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = ignbin($nt,$p); }
    return @ans;
}

#####################################################################
#			SEED HANDLER FUNCTIONS                      #
#####################################################################

sub random_seed_from_phrase { # Argument $phrase
    my $phrase = shift(@_);
    $phrase ||= "";
    return phrtsd($phrase);
}

sub random_get_seed { # no argument
    return getsd();
}

sub random_set_seed_from_phrase { # Argument $phrase
    my $phrase = shift(@_);
    $phrase ||= "";
    salfph($phrase);
    return 1;
}

sub random_set_seed { # Argument @seed
    my($seed1,$seed2) = @_;
    croak("Usage: random_set_seed(\@seed)\n\@seed[0,1] must be two integers ".
	  "in the range (1,1) to (2147483562,2147483398)\nand usually comes ".
	  "from a call to random_get_seed() ".
	  "or\nrandom_seed_from_phrase(\$phrase).")
	unless (((($seed1 == int($seed1)) and ($seed2 == int($seed2))) and
		 (($seed1 > 0) and ($seed2 > 0))) and
		(($seed1 < 2147483563) and ($seed2 < 2147483399)));
    setall($seed1,$seed2);
    return 1;
}

#####################################################################
#			   HELPER ROUTINES                          #
#    These use the C work arrays and are not intended for export    #
#	 (Currently only used in random_multivariate_normal)        #
#####################################################################

sub getflt {
    my $n = $_[0];
    my $val;
    my $i = 0;
    my @junk = (0) x $n;
    foreach $val (@junk) {
	$val = gvprfw($i);
	$i++;
    }
    return @junk;
}

sub putflt {
    my $n = scalar(@_);
    rsprfw($n) or return 0;
    my $val;
    my $i = 0;
    foreach $val (@_) { # load up floats
	svprfw($i,$val);
	$i++;
    }
    return 1;
}

# Autoload methods go after =cut, and are processed by the autosplit program.

1;

__END__