The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.

pp_bless('PDL::GSL::RNG'); # make the functions generated go into our namespace, and 
				# not PDL's namespace

pp_addpm({At=>Top},<<'EOD');
=head1 NAME

PDL::GSL::RNG - PDL interface to RNG and randist routines in GSL

=head1 DESCRIPTION

This is an interface to the rng and randist packages present
in the GNU Scientific Library.

=head1 SYNOPSIS

   use PDL;
   use PDL::GSL::RNG;

   $rng = PDL::GSL::RNG->new('taus');

   $rng->set_seed(time());

   $a=zeroes(5,5,5)

   $rng->get_uniform($a); # inplace

   $b=$rng->get_uniform(3,4,5); # creates new pdl

=head1 FUNCTIONS

=head2 new

=for ref

The new method initializes a new instance of the RNG.

The available RNGs are:

 coveyou cmrg fishman18 fishman20 fishman2x gfsr4 knuthran
 knuthran2 knuthran2002 lecuyer21 minstd mrg mt19937 mt19937_1999
 mt19937_1998 r250 ran0 ran1 ran2 ran3 rand rand48 random128_bsd
 random128_glibc2 random128_libc5 random256_bsd random256_glibc2
 random256_libc5 random32_bsd random32_glibc2 random32_libc5
 random64_bsd random64_glibc2 random64_libc5 random8_bsd
 random8_glibc2 random8_libc5 random_bsd random_glibc2
 random_libc5 randu ranf ranlux ranlux389 ranlxd1 ranlxd2 ranlxs0
 ranlxs1 ranlxs2 ranmar slatec taus taus2 taus113 transputer tt800
 uni uni32 vax waterman14 zuf default

The last one (default) uses the environment variable GSL_RNG_TYPE.

Note that only a few of these rngs are recommended for general
use. Please check the GSL documentation for more information.

=for usage

Usage:

   $blessed_ref = PDL::GSL::RNG->new($RNG_name);

Example:

=for example

   $rng = PDL::GSL::RNG->new('taus');

=head2 set_seed

=for ref

Sets the RNG seed.

Usage:

=for usage

   $rng->set_seed($integer);
   # or
   $rng = PDL::GSL::RNG->new('taus')->set_seed($integer);

Example:

=for example

   $rng->set_seed(666);

=head2 min

=for ref

Return the minimum value generable by this RNG.

Usage: 

=for usage

   $integer = $rng->min();

Example:

=for example

   $min = $rng->min(); $max = $rng->max();

=head2 max

=for ref

Return the maximum value generable by the RNG.

Usage: 

=for usage

   $integer = $rng->max();

Example:

=for example

   $min = $rng->min(); $max = $rng->max();

=head2 name

=for ref

Returns the name of the RNG.

Usage: 

=for usage

   $string = $rng->name();

Example:

=for example

   $name = $rng->name();

=head2 get

=for ref

This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get() returns integer values
between a minimum and a maximum specific to every RNG.

Usage: 

=for usage

   $piddle = $rng->get($list_of_integers)
   $rng->get($piddle);

Example:

=for example

   $a = zeroes 5,6;
   $o = $rng->get(10,10); $rng->get($a);

=head2 get_int

=for ref

This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get_int() returns integer values
between 0 and $max.

Usage: 

=for usage

   $piddle = $rng->get($max, $list_of_integers)
   $rng->get($max, $piddle);

Example:

=for example

   $a = zeroes 5,6; $max=100;
   $o = $rng->get(10,10); $rng->get($a);

=head2 get_uniform

=for ref

This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get_uniform() returns values 0<=x<1,

Usage: 

=for usage

   $piddle = $rng->get_uniform($list_of_integers)
   $rng->get_uniform($piddle);

Example:

=for example

   $a = zeroes 5,6; $max=100;
   $o = $rng->get_uniform(10,10); $rng->get_uniform($a);

=head2 get_uniform_pos

=for ref

This function creates a piddle with given dimensions or accept an
existing piddle and fills it. get_uniform_pos() returns values 0<x<1,

Usage: 

=for usage

   $piddle = $rng->get_uniform_pos($list_of_integers)
   $rng->get_uniform_pos($piddle);

Example:

=for example

   $a = zeroes 5,6;
   $o = $rng->get_uniform_pos(10,10); $rng->get_uniform_pos($a);

=head2 ran_shuffle

=for ref

Shuffles values in piddle

Usage:

=for usage

   $rng->ran_shuffle($piddle);

=head2 ran_shuffle_vec

=for ref

Shuffles values in piddle

Usage:

=for usage

   $rng->ran_shuffle_vec(@vec);

=head2 ran_choose

=for ref

Chooses values from C<$inpiddle> to C<$outpiddle>.

Usage:

=for usage

   $rng->ran_choose($inpiddle,$outpiddle);

=head2 ran_choose_vec

=for ref

Chooses C<$n> values from C<@vec>.

Usage:

=for usage

   @chosen = $rng->ran_choose_vec($n,@vec);

=head2 ran_gaussian

=for ref

Fills output piddle with random values from Gaussian distribution with mean zero and standard deviation C<$sigma>.

Usage:

=for usage

 $piddle = $rng->ran_gaussian($sigma,[list of integers = output piddle dims]);
 $rng->ran_gaussian($sigma, $output_piddle);

Example:

=for example

  $o = $rng->ran_gaussian($sigma,10,10);
  $rng->ran_gaussian($sigma,$a);

=head2 ran_gaussian_var

=for ref

This method is similar to L<ran_gaussian|/ran_gaussian> except that it takes
the parameters of the distribution as a piddle and returns a piddle of equal
dimensions.

Usage:

=for usage

   $piddle = $rng->ran_gaussian_var($sigma_piddle);
   $rng->ran_gaussian_var($sigma_piddle, $output_piddle);

Example:

=for example

   $sigma_pdl = rvals zeroes 11,11;
   $o = $rng->ran_gaussian_var($sigma_pdl);

=head2 ran_additive_gaussian

=for ref

Add Gaussian noise of given sigma to a piddle.

Usage:

=for usage

   $rng->ran_additive_gaussian($sigma,$piddle);

Example:

=for example

   $rng->ran_additive_gaussian(1,$image);

=head2 ran_bivariate_gaussian

=for ref

Generates C<$n> bivariate gaussian random deviates.

Usage:

=for usage

   $piddle = $rng->ran_bivariate_gaussian($sigma_x,$sigma_y,$rho,$n);

Example:

=for example

   $o = $rng->ran_bivariate_gaussian(1,2,0.5,1000);

=head2 ran_poisson

=for ref

Fills output piddle by with random integer values from the Poisson distribution with mean C<$mu>.

Usage:

=for usage

   $piddle = $rng->ran_poisson($mu,[list of integers = output piddle dims]);
   $rng->ran_poisson($mu,$output_piddle);

=head2 ran_poisson_var

=for ref

Similar to L<ran_poisson|/ran_poisson> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_poisson_var($mu_piddle);

=head2 ran_additive_poisson

=for ref

Add Poisson noise of given C<$mu> to a C<$piddle>.

Usage:

=for usage

   $rng->ran_additive_poisson($mu,$piddle);

Example:

=for example

   $rng->ran_additive_poisson(1,$image);

=head2 ran_feed_poisson

=for ref

This method simulates shot noise, taking the values of piddle as
values for C<$mu> to be fed in the poissonian RNG.

Usage:

=for usage

   $rng->ran_feed_poisson($piddle);

Example:

=for example

   $rng->ran_feed_poisson($image);

=head2 ran_bernoulli

=for ref

Fills output piddle with random values 0 or 1, the result of a Bernoulli trial with probability C<$p>.

Usage:

=for usage

   $piddle = $rng->ran_bernoulli($p,[list of integers = output piddle dims]);
   $rng->ran_bernoulli($p,$output_piddle);

=head2 ran_bernoulli_var

=for ref

Similar to L<ran_bernoulli|/ran_bernoulli> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_bernoulli_var($p_piddle);

=head2 ran_beta

=for ref

Fills output piddle with random variates from the beta distribution with parameters C<$a> and C<$b>.

Usage:

=for usage

   $piddle = $rng->ran_beta($a,$b,[list of integers = output piddle dims]);
   $rng->ran_beta($a,$b,$output_piddle);

=head2 ran_beta_var

=for ref

Similar to L<ran_beta|/ran_beta> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_beta_var($a_piddle, $b_piddle);

=head2 ran_binomial

=for ref

Fills output piddle with random integer values from the binomial distribution, the number of
successes in C<$n> independent trials with probability C<$p>.

Usage:

=for usage

   $piddle = $rng->ran_binomial($p,$n,[list of integers = output piddle dims]);
   $rng->ran_binomial($p,$n,$output_piddle);

=head2 ran_binomial_var

=for ref

Similar to L<ran_binomial|/ran_binomial> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_binomial_var($p_piddle, $n_piddle);

=head2 ran_cauchy

=for ref

Fills output piddle with random variates from the Cauchy distribution with scale parameter C<$a>.

Usage:

=for usage

   $piddle = $rng->ran_cauchy($a,[list of integers = output piddle dims]);
   $rng->ran_cauchy($a,$output_piddle);

=head2 ran_cauchy_var

=for ref

Similar to L<ran_cauchy|/ran_cauchy> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_cauchy_var($a_piddle);

=head2 ran_chisq

=for ref

Fills output piddle with random variates from the chi-squared distribution with
C<$nu> degrees of freedom.

Usage:

=for usage

   $piddle = $rng->ran_chisq($nu,[list of integers = output piddle dims]);
   $rng->ran_chisq($nu,$output_piddle);

=head2 ran_chisq_var

=for ref

Similar to L<ran_chisq|/ran_chisq> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_chisq_var($nu_piddle);

=head2 ran_exponential

=for ref

Fills output piddle with random variates from the exponential distribution with mean C<$mu>.

Usage:

=for usage

   $piddle = $rng->ran_exponential($mu,[list of integers = output piddle dims]);
   $rng->ran_exponential($mu,$output_piddle);

=head2 ran_exponential_var

=for ref

Similar to L<ran_exponential|/ran_exponential> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_exponential_var($mu_piddle);

=head2 ran_exppow

=for ref

Fills output piddle with random variates from the exponential power distribution with scale
parameter C<$a> and exponent C<$b>.

Usage:

=for usage

   $piddle = $rng->ran_exppow($mu,$a,[list of integers = output piddle dims]);
   $rng->ran_exppow($mu,$a,$output_piddle);

=head2 ran_exppow_var

=for ref

Similar to L<ran_exppow|/ran_exppow> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_exppow_var($mu_piddle, $a_piddle);

=head2 ran_fdist

=for ref

Fills output piddle with random variates from the F-distribution with degrees
of freedom C<$nu1> and C<$nu2>.

Usage:

=for usage

   $piddle = $rng->ran_fdist($nu1, $nu2,[list of integers = output piddle dims]);
   $rng->ran_fdist($nu1, $nu2,$output_piddle);

=head2 ran_fdist_var

=for ref

Similar to L<ran_fdist|/ran_fdist> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_fdist_var($nu1_piddle, $nu2_piddle);

=head2 ran_flat

=for ref

Fills output piddle with random variates from the flat (uniform) distribution from C<$a> to C<$b>.

Usage:

=for usage

   $piddle = $rng->ran_flat($a,$b,[list of integers = output piddle dims]);
   $rng->ran_flat($a,$b,$output_piddle);

=head2 ran_flat_var

=for ref

Similar to L<ran_flat|/ran_flat> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_flat_var($a_piddle, $b_piddle);

=head2 ran_gamma

=for ref

Fills output piddle with random variates from the gamma distribution.

Usage:

=for usage

   $piddle = $rng->ran_gamma($a,$b,[list of integers = output piddle dims]);
   $rng->ran_gamma($a,$b,$output_piddle);

=head2 ran_gamma_var

=for ref

Similar to L<ran_gamma|/ran_gamma> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_gamma_var($a_piddle, $b_piddle);

=head2 ran_geometric

=for ref

Fills output piddle with random integer values from the geometric distribution,
the number of independent trials with probability C<$p> until the first success.

Usage:

=for usage

   $piddle = $rng->ran_geometric($p,[list of integers = output piddle dims]);
   $rng->ran_geometric($p,$output_piddle);

=head2 ran_geometric_var

=for ref

Similar to L<ran_geometric|/ran_geometric> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_geometric_var($p_piddle);

=head2 ran_gumbel1

=for ref

Fills output piddle with random variates from the Type-1 Gumbel distribution.

Usage:

=for usage

   $piddle = $rng->ran_gumbel1($a,$b,[list of integers = output piddle dims]);
   $rng->ran_gumbel1($a,$b,$output_piddle);

=head2 ran_gumbel1_var

=for ref

Similar to L<ran_gumbel1|/ran_gumbel1> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_gumbel1_var($a_piddle, $b_piddle);

=head2 ran_gumbel2

=for ref

Fills output piddle with random variates from the Type-2 Gumbel distribution.

Usage:

=for usage

   $piddle = $rng->ran_gumbel2($a,$b,[list of integers = output piddle dims]);
   $rng->ran_gumbel2($a,$b,$output_piddle);

=head2 ran_gumbel2_var

=for ref

Similar to L<ran_gumbel2|/ran_gumbel2> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_gumbel2_var($a_piddle, $b_piddle);

=head2 ran_hypergeometric

=for ref

Fills output piddle with random integer values from the hypergeometric distribution.
If a population contains C<$n1> elements of type 1 and C<$n2> elements of
type 2 then the hypergeometric distribution gives the probability of obtaining
C<$x> elements of type 1 in C<$t> samples from the population without replacement.

Usage:

=for usage

   $piddle = $rng->ran_hypergeometric($n1, $n2, $t,[list of integers = output piddle dims]);
   $rng->ran_hypergeometric($n1, $n2, $t,$output_piddle);

=head2 ran_hypergeometric_var

=for ref

Similar to L<ran_hypergeometric|/ran_hypergeometric> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_hypergeometric_var($n1_piddle, $n2_piddle, $t_piddle);

=head2 ran_laplace

=for ref

Fills output piddle with random variates from the Laplace distribution with width C<$a>.

Usage:

=for usage

   $piddle = $rng->ran_laplace($a,[list of integers = output piddle dims]);
   $rng->ran_laplace($a,$output_piddle);

=head2 ran_laplace_var

=for ref

Similar to L<ran_laplace|/ran_laplace> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_laplace_var($a_piddle);

=head2 ran_levy

=for ref

Fills output piddle with random variates from the Levy symmetric stable
distribution with scale C<$c> and exponent C<$alpha>.

Usage:

=for usage

   $piddle = $rng->ran_levy($mu,$a,[list of integers = output piddle dims]);
   $rng->ran_levy($mu,$a,$output_piddle);

=head2 ran_levy_var

=for ref

Similar to L<ran_levy|/ran_levy> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_levy_var($mu_piddle, $a_piddle);

=head2 ran_logarithmic

=for ref

Fills output piddle with random integer values from the logarithmic distribution.

Usage:

=for usage

   $piddle = $rng->ran_logarithmic($p,[list of integers = output piddle dims]);
   $rng->ran_logarithmic($p,$output_piddle);

=head2 ran_logarithmic_var

=for ref

Similar to L<ran_logarithmic|/ran_logarithmic> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_logarithmic_var($p_piddle);

=head2 ran_logistic

=for ref

Fills output piddle with random random variates from the logistic distribution.

Usage:

=for usage

   $piddle = $rng->ran_logistic($m,[list of integers = output piddle dims]u)
   $rng->ran_logistic($m,$output_piddle)

=head2 ran_logistic_var

=for ref

Similar to L<ran_logistic|/ran_logistic> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_logistic_var($m_piddle);

=head2 ran_lognormal

=for ref

Fills output piddle with random variates from the lognormal distribution with
parameters C<$mu> (location) and C<$sigma> (scale).

Usage:

=for usage

   $piddle = $rng->ran_lognormal($mu,$sigma,[list of integers = output piddle dims]);
   $rng->ran_lognormal($mu,$sigma,$output_piddle);

=head2 ran_lognormal_var

=for ref

Similar to L<ran_lognormal|/ran_lognormal> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_lognormal_var($mu_piddle, $sigma_piddle);

=head2 ran_negative_binomial

=for ref

Fills output piddle with random integer values from the negative binomial
distribution, the number of failures occurring before C<$n> successes in
independent trials with probability C<$p> of success. Note that C<$n> is
not required to be an integer.

Usage:

=for usage

   $piddle = $rng->ran_negative_binomial($p,$n,[list of integers = output piddle dims]);
   $rng->ran_negative_binomial($p,$n,$output_piddle);

=head2 ran_negative_binomial_var

=for ref

Similar to L<ran_negative_binomial|/ran_negative_binomial> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_negative_binomial_var($p_piddle, $n_piddle);

=head2 ran_pareto

=for ref

Fills output piddle with random variates from the Pareto distribution of
order C<$a> and scale C<$b>.

Usage:

=for usage

   $piddle = $rng->ran_pareto($a,$b,[list of integers = output piddle dims]);
   $rng->ran_pareto($a,$b,$output_piddle);

=head2 ran_pareto_var

=for ref

Similar to L<ran_pareto|/ran_pareto> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_pareto_var($a_piddle, $b_piddle);

=head2 ran_pascal

=for ref

Fills output piddle with random integer values from the Pascal distribution.
The Pascal distribution is simply a negative binomial distribution
(see L<ran_negative_binomial|/ran_negative_binomial>) with an integer value of C<$n>.

Usage:

=for usage

   $piddle = $rng->ran_pascal($p,$n,[list of integers = output piddle dims]);
   $rng->ran_pascal($p,$n,$output_piddle);

=head2 ran_pascal_var

=for ref

Similar to L<ran_pascal|/ran_pascal> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_pascal_var($p_piddle, $n_piddle);

=head2 ran_rayleigh

=for ref

Fills output piddle with random variates from the Rayleigh distribution with scale parameter C<$sigma>.

Usage:

=for usage

   $piddle = $rng->ran_rayleigh($sigma,[list of integers = output piddle dims]);
   $rng->ran_rayleigh($sigma,$output_piddle);

=head2 ran_rayleigh_var

=for ref

Similar to L<ran_rayleigh|/ran_rayleigh> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_rayleigh_var($sigma_piddle);

=head2 ran_rayleigh_tail

=for ref

Fills output piddle with random variates from the tail of the Rayleigh distribution
with scale parameter C<$sigma> and a lower limit of C<$a>.

Usage:

=for usage

   $piddle = $rng->ran_rayleigh_tail($a,$sigma,[list of integers = output piddle dims]);
   $rng->ran_rayleigh_tail($a,$sigma,$output_piddle);

=head2 ran_rayleigh_tail_var

=for ref

Similar to L<ran_rayleigh_tail|/ran_rayleigh_tail> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_rayleigh_tail_var($a_piddle, $sigma_piddle);

=head2 ran_tdist

=for ref

Fills output piddle with random variates from the t-distribution (AKA Student's
t-distribution) with C<$nu> degrees of freedom.

Usage:

=for usage

   $piddle = $rng->ran_tdist($nu,[list of integers = output piddle dims]);
   $rng->ran_tdist($nu,$output_piddle);

=head2 ran_tdist_var

=for ref

Similar to L<ran_tdist|/ran_tdist> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_tdist_var($nu_piddle);

=head2 ran_ugaussian_tail

=for ref

Fills output piddle with random variates from the upper tail of a Gaussian
distribution with C<standard deviation = 1> (AKA unit Gaussian distribution).

Usage:

=for usage

   $piddle = $rng->ran_ugaussian_tail($tail,[list of integers = output piddle dims]);
   $rng->ran_ugaussian_tail($tail,$output_piddle);

=head2 ran_ugaussian_tail_var

=for ref

Similar to L<ran_ugaussian_tail|/ran_ugaussian_tail> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_ugaussian_tail_var($tail_piddle);

=head2 ran_weibull

=for ref

Fills output piddle with random variates from the Weibull distribution.

Usage:

=for usage

   $piddle = $rng->ran_weibull($mu,$a,[list of integers = output piddle dims]);
   $rng->ran_weibull($mu,$a,$output_piddle);

=head2 ran_weibull_var

=for ref

Similar to L<ran_weibull|/ran_weibull> except that it takes the distribution
parameters as a piddle and returns a piddle of equal dimensions.

Usage:

=for usage

   $piddle = $rng->ran_weibull_var($mu_piddle, $a_piddle);

=head2 ran_dir

=for ref

Returns C<$n> random vectors in C<$ndim> dimensions.

Usage:

=for usage

   $piddle = $rng->ran_dir($ndim,$n);

Example:

=for example

   $o = $rng->ran_dir($ndim,$n);

=head2 ran_discrete_preproc

=for ref

This method returns a handle that must be used when calling
L<ran_discrete|/ran_discrete>. You specify the probability of the integer number
that are returned by L<ran_discrete|/ran_discrete>.

Usage:

=for usage

   $discrete_dist_handle = $rng->ran_discrete_preproc($double_piddle_prob);

Example:

=for example

   $prob = pdl [0.1,0.3,0.6];
   $ddh = $rng->ran_discrete_preproc($prob);
   $o = $rng->ran_discrete($discrete_dist_handle,100);

=head2 ran_discrete

=for ref

Is used to get the desired samples once a proper handle has been
enstablished (see ran_discrete_preproc()).

Usage:

=for usage

   $piddle = $rng->ran_discrete($discrete_dist_handle,$num);

Example:

=for example

   $prob = pdl [0.1,0.3,0.6];
   $ddh = $rng->ran_discrete_preproc($prob);
   $o = $rng->ran_discrete($discrete_dist_handle,100);

=head2 ran_ver

=for ref

Returns a piddle with C<$n> values generated by the Verhulst map from C<$x0> and
parameter C<$r>.

Usage:

=for usage

   $rng->ran_ver($x0, $r, $n);

=head2 ran_caos

=for ref

Returns values from Verhuls map with C<$r=4.0> and randomly chosen
C<$x0>. The values are scaled by C<$m>.

Usage:

=for usage

   $rng->ran_caos($m,$n);

=head1 BUGS

Feedback is welcome. Log bugs in the PDL bug database (the
database is always linked from L<http://pdl.perl.org/>).

=head1 SEE ALSO

L<PDL>

The GSL documentation is online at
L<http://www.gnu.org/software/gsl/manual/html_node/>

=head1 AUTHOR

This file copyright (C) 1999 Christian Pellegrin <chri@infis.univ.trieste.it>
Docs mangled by C. Soeller. All rights reserved. There
is no warranty. You are allowed to redistribute this software /
documentation under certain conditions. For details, see the file
COPYING in the PDL distribution. If this file is separated from the
PDL distribution, the copyright notice should be included in the file.

The GSL RNG and randist modules were written by James Theiler.

=cut


EOD

# PP interface to RNG


##############################
# 
# make_get_sub generates a wrapper PDL subroutine that handles the 
# fill-a-PDL and create-a-PDL cases for each of the GSL functions.
#  --CED 
#
sub make_get_sub {
  my ($fname,$par) =@_;
  my $s;

  $s = '
sub ' . $fname . ' {
my ($obj,' . $par . '@var) = @_;';
   
    if ($par ne '') {
       my $ss=$par;

       $ss =~ s/,//;
       $s .= 'if (!(' . $ss . '>0)) {barf("first parameter must be an int >0")};';
    }

$s .= 'if (ref($var[0]) eq \'PDL\') {
    gsl_' . $fname . '_meat($var[0],' . $par . '$$obj);
    return $var[0]; 
}
else {
    my $p;

    $p = zeroes @var;
    gsl_' . $fname . '_meat($p,' . $par . '$$obj);
    return $p;
}
}
'
}

pp_addpm(<<'EOPM');

use strict;

# PDL::GSL::RNG::nullcreate just creates a null PDL. Used
#  for the GSL functions that create PDLs
sub nullcreate{
	
	my ($type,$arg) = @_;

	PDL->nullcreate($arg);
}
EOPM

pp_addpm(make_get_sub('get_uniform',''));
pp_addpm(make_get_sub('get_uniform_pos',''));
pp_addpm(make_get_sub('get',''));
pp_addpm(make_get_sub('get_int','$n,'));

pp_addhdr('
#include <string.h>

#include "gsl/gsl_rng.h"
#include "gsl/gsl_randist.h"


');

sub pp_defnd { # hide the docs
  my ($name, %hash) = @_;
  pp_def($name,%hash,Doc=>undef);
}

pp_defnd(
       'gsl_get_uniform_meat',
       Pars => '[o]a()',
       GenericTypes => [F,D],
       OtherPars => 'IV rng',
       Code => '
$a() = gsl_rng_uniform(INT2PTR(gsl_rng *, $COMP(rng)));');

pp_defnd(
       'gsl_get_uniform_pos_meat',
       Pars => '[o]a()',
       GenericTypes => [F,D],
       OtherPars => 'IV rng',
       Code => '
$a() = gsl_rng_uniform_pos(INT2PTR(gsl_rng *, $COMP(rng)));');

pp_defnd(
       'gsl_get_meat',
       Pars => '[o]a()',
       OtherPars => 'IV rng',
       Code => '
$a() = gsl_rng_get(INT2PTR(gsl_rng *, $COMP(rng)));');

pp_defnd(
       'gsl_get_int_meat',
       Pars => '[o]a()',
       OtherPars => 'int n; IV rng',
       Code => '
$a() = gsl_rng_uniform_int(INT2PTR(gsl_rng *, $COMP(rng)),$COMP(n));');

# randist stuff

sub add_randist {
  my ($name,$npar) = @_;
  my ($pars1,$fcall1,$arglist);
  
  if ($npar==1) {
    $pars1='double a; IV rng';
    $fcall1='$COMP(a)';
    $arglist='$a,';
    $pars2='a()';
    $fcall2='$a()';
  }
  if ($npar==2) {
    $pars1='double a; double b; IV rng';
    $fcall1='$COMP(a),$COMP(b)';
    $arglist='$a,$b,';
    $pars2='a();b()';
    $fcall2='$a(),$b()';
  }
  if ($npar==3) {
    $pars1='double a; double b; double c; IV rng';
    $fcall1='$COMP(a),$COMP(b),$COMP(c)';
    $arglist='$a,$b,$c,';
    $pars2='a();b();c()';
    $fcall2='$a(),$b(),$c()';
  }

  pp_defnd(
	 'ran_' . $name . '_meat',
	 Pars => '[o]x()',
	 OtherPars => $pars1,
	 Code =>'
$x() = gsl_ran_' . $name . '(INT2PTR(gsl_rng *, $COMP(rng)),' . $fcall1 . ');');

  pp_addpm('
sub ran_' . $name . ' {
my ($obj,' . $arglist . '@var) = @_;
if (ref($var[0]) eq \'PDL\') {
    ran_' . $name . '_meat($var[0],' . $arglist . '$$obj);
    return $var[0]; 
}
else {
    my $p;

    $p = zeroes @var;
    ran_' . $name . '_meat($p,' . $arglist . '$$obj);
    return $p;
}
}
');

  pp_defnd(
	 'ran_' . $name . '_var_meat',
	 Pars => $pars2 . ';[o]x()',
	 OtherPars => 'IV rng',
	 Code =>'
$x() = gsl_ran_' . $name . '(INT2PTR(gsl_rng *, $COMP(rng)),' . $fcall2 . ');');

  pp_addpm('
sub ran_' . $name . '_var {
my ($obj,@var) = @_;
    if (scalar(@var) != ' . $npar . ') {barf("Bad number of parameters!");}
    return ran_' . $name . '_var_meat(@var,$$obj);
}
');

#  pp_defnd(
#	 'ran_' . $name . '_add_meat',
#	 Pars => '[o]x()',
#	 OtherPars => $pars1,
#	 Code =>'
#$x() += gsl_ran_' . $name . '(INT2PTR(gsl_rng *, $COMP(rng)),' . $fcall1 . ');');

#  pp_addpm('
#sub ran_' . $name . '_add {
#my ($obj,' . $arglist . '@var) = @_;
#if (ref($var[0]) eq \'PDL\') {
#    PDL::ran_' . $name . '_add_meat($var[0],' . $arglist . '$$obj);
#    return $var[0]; 
#}
#else {
#    barf("In add mode you must specify a piddle!");
#}
#}
#');

#  if ($npar==1) {
#    pp_defnd(
#	   'ran_' . $name . '_feed_meat',
#	   Pars => '[o]x()',
#	   OtherPars => 'IV rng',
#	   Code =>'
#$x() = gsl_ran_' . $name . '(INT2PTR(gsl_rng *, $COMP(rng)), $x());');

#    pp_addpm('
#sub ran_' . $name . '_feed {
#my ($obj, @var) = @_;
#if (ref($var[0]) eq \'PDL\') {
#    PDL::ran_' . $name . '_feed_meat($var[0], $$obj);
#    return $var[0]; 
#}
#else {
#    barf("In feed mode you must specify a piddle!");
#}
#}
#');
#  }
  
}

add_randist('gaussian',1);
add_randist('ugaussian_tail',1);
add_randist('exponential',1);
add_randist('laplace',1);
add_randist('exppow',2);
add_randist('cauchy',1);
add_randist('rayleigh',1);
add_randist('rayleigh_tail',2);
add_randist('levy',2);
add_randist('gamma',2);
add_randist('flat',2);
add_randist('lognormal',2);
add_randist('chisq',1);
add_randist('fdist',2);
add_randist('tdist',1);
add_randist('beta',2);
add_randist('logistic',1);
add_randist('pareto',2);
add_randist('weibull',2);
add_randist('gumbel1',2);
add_randist('gumbel2',2);
add_randist('poisson',1);
add_randist('bernoulli',1);
add_randist('binomial',2);
add_randist('negative_binomial',2);
add_randist('pascal',2);
add_randist('geometric',1);
add_randist('hypergeometric',3);
add_randist('logarithmic',1);

# specific rnadist

pp_defnd(
       'ran_additive_gaussian_meat',
       Pars => ';[o]x()',
       OtherPars => 'double sigma; IV rng',
       Code =>'$x() += gsl_ran_gaussian(INT2PTR(gsl_rng *, $COMP(rng)), $COMP(sigma));');

pp_addpm('
       sub ran_additive_gaussian {
	 my ($obj,$sigma,@var) = @_;
	 if (ref($var[0]) eq \'PDL\') {
	   ran_additive_gaussian_meat($var[0],$sigma,$$obj);
	   return $var[0]; 
	 }
	 else {
	   barf("In additive gaussian mode you must specify a piddle!");
	 }
       }
       ');

pp_defnd(
       'ran_additive_poisson_meat',
       Pars => ';[o]x()',
       OtherPars => 'double sigma; IV rng',
       Code =>'$x() += gsl_ran_poisson(INT2PTR(gsl_rng *, $COMP(rng)), $COMP(sigma));');

pp_addpm('
       sub ran_additive_poisson {
	 my ($obj,$sigma,@var) = @_;
	 if (ref($var[0]) eq \'PDL\') {
	   ran_additive_poisson_meat($var[0],$sigma,$$obj);
	   return $var[0]; 
	 }
	 else {
	   barf("In additive poisson mode you must specify a piddle!");
	 }
       }
       ');

pp_defnd(
       'ran_feed_poisson_meat',
       Pars => ';[o]x()',
       OtherPars => 'IV rng',
       Code =>'$x() = gsl_ran_poisson(INT2PTR(gsl_rng *, $COMP(rng)), $x());');

pp_addpm('
       sub ran_feed_poisson {
	 my ($obj,@var) = @_;
	 if (ref($var[0]) eq \'PDL\') {
	   ran_feed_poisson_meat($var[0],$$obj);
	   return $var[0]; 
	 }
	 else {
	   barf("In poisson mode you must specify a piddle!");
	 }
       }
       ');

pp_defnd(
	'ran_bivariate_gaussian_meat',
	Pars => ';[o]x(n)',
	OtherPars => 'double sigma_x; double sigma_y; double rho; IV rng',
        Code =>'
double xx,yy;
gsl_ran_bivariate_gaussian(INT2PTR(gsl_rng *, $COMP(rng)), $COMP(sigma_x), $COMP(sigma_y),$COMP(rho), &xx, &yy);
$x(n=>0)=xx;
$x(n=>1)=yy;
');	

pp_addpm('
       sub ran_bivariate_gaussian {
	 my ($obj,$sigma_x,$sigma_y,$rho,$n) = @_;
	 if ($n>0) {
	   my $p = zeroes(2,$n);
	   ran_bivariate_gaussian_meat($p,$sigma_x,$sigma_y,$rho,$$obj);
	   return $p; 
	 }
	 else {
	   barf("Not enough parameters for gaussian bivariate!");
	 }
       }
       ');

pp_defnd(
	'ran_dir_2d_meat',
	Pars => ';[o]x(n)',
	OtherPars => 'IV rng',
        Code =>'
double xx,yy;
gsl_ran_dir_2d(INT2PTR(gsl_rng *, $COMP(rng)), &xx, &yy);
$x(n=>0)=xx;
$x(n=>1)=yy;
');	

pp_defnd(
	'ran_dir_3d_meat',
	Pars => ';[o]x(n)',
	OtherPars => 'IV rng',
        Code =>'
double xx,yy,zz;
gsl_ran_dir_3d(INT2PTR(gsl_rng *, $COMP(rng)), &xx, &yy, &zz);
$x(n=>0)=xx;
$x(n=>1)=yy;
$x(n=>2)=zz;
');	

$MAX_DIMENSIONS = 100;

pp_defnd(
	'ran_dir_nd_meat',
	Pars => ';[o]x(n)',
	OtherPars => 'int ns => n; IV rng',
        Code =>'
double xxx[' . $MAX_DIMENSIONS .']; 
gsl_ran_dir_nd(INT2PTR(gsl_rng *, $COMP(rng)), $COMP(ns), xxx);
loop (n) %{ $x() = xxx[n]; %}');	

pp_addpm('
       sub ran_dir {
	 my ($obj,$ndim,$n) = @_;
	 if ($n>0) {
	   my $p = zeroes($ndim,$n);
	   if ($ndim==2) { ran_dir_2d_meat($p,$$obj); }
	   elsif ($ndim==3) { ran_dir_3d_meat($p,$$obj); }
	   elsif ($ndim>=4 && $ndim<=' . $MAX_DIMENSIONS . ') { ran_dir_nd_meat($p,$ndim,$$obj); }
	   else { barf("Bad number of dimensions!"); }
	   return $p; 
	 }
	 else {
	   barf("Not enough parameters for random vectors!");
	 }
       }
       ');

pp_defnd(
	'ran_discrete_meat',
	Pars => ';[o]x()',
	OtherPars => 'IV rng_discrete; IV rng',
        Code =>'
$x()=gsl_ran_discrete(INT2PTR(gsl_rng *, $COMP(rng)), INT2PTR(gsl_ran_discrete_t *, $COMP(rng_discrete))); ');	

  pp_addpm('
sub ran_discrete {
my ($obj, $rdt, @var) = @_;
if (ref($var[0]) eq \'PDL\') {
    ran_discrete_meat($var[0], $$rdt, $$obj);
    return $var[0]; 
}
else {
    my $p;

    $p = zeroes @var;
    ran_discrete_meat($p, $$rdt, $$obj);
    return $p;
}
}
');

  pp_addpm('
sub ran_shuffle_vec {
my ($obj,@in) = @_;
my (@out,$i,$p);

$p = long [0..$#in];
$obj->ran_shuffle($p);
for($i=0;$i<scalar(@in);$i++) {
$out[$p->at($i)]=$in[$i];
}
return @out;
}
');

  pp_addpm('
sub ran_choose_vec {
my ($obj,$nout,@in) = @_;
my (@out,$i,$pin,$pout);

$pin = long [0..$#in];
$pout = long [0..($nout-1)];
$obj->ran_choose($pin,$pout);
for($i=0;$i<$nout;$i++) {
$out[$i]=$in[$pout->at($i)];
}
return @out;
}
');


pp_defnd(
	'ran_ver_meat',
	Pars => ';[o]x(n)',
	OtherPars => 'double x0; double r;int ns => n; IV rng',
        Code =>'
double xx=$COMP(x0);

loop (n) %{ $x() = xx; xx = $COMP(r)*(1-xx)*xx; %}');	

pp_defnd(
	'ran_caos_meat',
	Pars => ';[o]x(n)',
	OtherPars => 'double m; int ns => n; IV rng',
        Code =>'
double xx=gsl_ran_gaussian(INT2PTR(gsl_rng *, $COMP(rng)),0.1)+0.5;

loop (n) %{ $x() = (xx-0.5)*$COMP(m); xx = 4.0*(1-xx)*xx; %}');	

pp_addpm('
       sub ran_ver {
	 my ($obj,$x0,$r,$n) = @_;
	 if ($n>0) {
	   my $p = zeroes($n);
	   ran_ver_meat($p,$x0,$r,$n,$$obj);
	   return $p; 
	 }
	 else {
	   barf("Not enough parameters for ran_ver!");
	 }
       }
       ');

pp_addpm('
       sub ran_caos {
	 my ($obj,$m,$n) = @_;
	 if ($n>0) {
	   my $p = zeroes($n);
	   ran_caos_meat($p,$m,$n,$$obj);
	   return $p; 
	 }
	 else {
	   barf("Not enough parameters for ran_caos!");
	 }
       }
       ');

# XS function for the RNG object

pp_addxs('','
MODULE = PDL::GSL::RNG PACKAGE = PDL::GSL::RNG

#define DEF_RNG(X) if (!strcmp(TYPE,#X)) rng=gsl_rng_alloc( gsl_rng_ ## X ); strcat(rngs,#X ", ");

gsl_rng *
new (CLASS,TYPE)
  char *CLASS
  char *TYPE
 CODE:
  gsl_rng * rng = NULL;
  char rngs[5000];
strcpy(rngs,"");
DEF_RNG(borosh13)
DEF_RNG(coveyou)
DEF_RNG(cmrg)
DEF_RNG(fishman18)
DEF_RNG(fishman20)
DEF_RNG(fishman2x)
DEF_RNG(gfsr4)
DEF_RNG(knuthran)
DEF_RNG(knuthran2)
DEF_RNG(knuthran2002)
DEF_RNG(lecuyer21)
DEF_RNG(minstd)
DEF_RNG(mrg)
DEF_RNG(mt19937)
DEF_RNG(mt19937_1999)
DEF_RNG(mt19937_1998)
DEF_RNG(r250)
DEF_RNG(ran0)
DEF_RNG(ran1)
DEF_RNG(ran2)
DEF_RNG(ran3)
DEF_RNG(rand)
DEF_RNG(rand48)
DEF_RNG(random128_bsd)
DEF_RNG(random128_glibc2)
DEF_RNG(random128_libc5)
DEF_RNG(random256_bsd)
DEF_RNG(random256_glibc2)
DEF_RNG(random256_libc5)
DEF_RNG(random32_bsd)
DEF_RNG(random32_glibc2)
DEF_RNG(random32_libc5)
DEF_RNG(random64_bsd)
DEF_RNG(random64_glibc2)
DEF_RNG(random64_libc5)
DEF_RNG(random8_bsd)
DEF_RNG(random8_glibc2)
DEF_RNG(random8_libc5)
DEF_RNG(random_bsd)
DEF_RNG(random_glibc2)
DEF_RNG(random_libc5)
DEF_RNG(randu)
DEF_RNG(ranf)
DEF_RNG(ranlux)
DEF_RNG(ranlux389)
DEF_RNG(ranlxd1)
DEF_RNG(ranlxd2)
DEF_RNG(ranlxs0)
DEF_RNG(ranlxs1)
DEF_RNG(ranlxs2)
DEF_RNG(ranmar)
DEF_RNG(slatec)
DEF_RNG(taus)
DEF_RNG(taus2)
DEF_RNG(taus113)
DEF_RNG(transputer)
DEF_RNG(tt800)
DEF_RNG(uni)
DEF_RNG(uni32)
DEF_RNG(vax)
DEF_RNG(waterman14)
DEF_RNG(zuf)
DEF_RNG(default)
  if (rng==NULL) {
    barf("Unknown RNG, plese use one of the following: %s", rngs);
  }
  else
  RETVAL = rng;
 OUTPUT:
  RETVAL

void
set_seed(rng, seed)
  gsl_rng * rng
  int seed
 PPCODE:
  gsl_rng_set(rng,seed);
  XPUSHs(ST(0)); /* return self */

unsigned int
min(rng)
  gsl_rng * rng
 CODE:
  RETVAL = gsl_rng_min(rng);
 OUTPUT:
  RETVAL

unsigned int
max(rng)
  gsl_rng * rng
 CODE:
  RETVAL = gsl_rng_max(rng);
 OUTPUT:
  RETVAL

char*
name(rng)
  gsl_rng * rng
 CODE:
  RETVAL = (char *) gsl_rng_name(rng);
 OUTPUT:
  RETVAL

void
DESTROY(sv)
  SV * sv
 CODE:
  gsl_rng *rng = INT2PTR(gsl_rng *, SvIV((SV*)SvRV(sv)));
  /* fprintf(stderr,"Freeing %d\n",rng); */
  gsl_rng_free((gsl_rng *) rng);

gsl_ran_discrete_t *
ran_discrete_preproc(rng, p)
  gsl_rng * rng
  pdl * p
     CODE:
       int n;

       if (p->ndims!=1 || p->datatype!=PDL_D) {
	 barf("Bad input to ran_discrete_preproc!");
       }
       n = p->dims[0];
       PDL->make_physical(p);
       RETVAL = gsl_ran_discrete_preproc(n,(double *) p->data);
     OUTPUT:
       RETVAL

void
ran_shuffle(rng, in)
  gsl_rng * rng
  pdl * in
 CODE:
  int size, n;

  n = in->nvals;
  PDL->make_physical(in);
  switch(in->datatype) {
   case PDL_INVALID: barf("ran_shuffle was passed a piddle of type PDL_INVALID"); break;
   case PDL_B: size=sizeof(PDL_Byte); break;
   case PDL_S: size=sizeof(PDL_Short); break;
   case PDL_US: size=sizeof(PDL_Ushort); break;
   case PDL_L: size=sizeof(PDL_Long); break;
   case PDL_IND: size=sizeof(PDL_Indx); break;
   case PDL_LL: size=sizeof(PDL_LongLong); break;
   case PDL_F: size=sizeof(PDL_Float); break;
   case PDL_D: size=sizeof(PDL_Double); break;
  }
  gsl_ran_shuffle(rng,(double *) in->data,n,size);

void
ran_choose(rng, in, out)
  gsl_rng * rng
  pdl * in
  pdl * out
 CODE:
  int size, n,m;

  n = in->nvals;
  m = out->nvals;
  if (in->datatype != out->datatype) barf("Data Types must match for ran_chooser");
  PDL->make_physical(in);
  PDL->make_physical(out);
  switch(in->datatype) {
   case PDL_INVALID: barf("ran_choose was passed a piddle of type PDL_INVALID"); break;
   case PDL_B: size=sizeof(PDL_Byte); break;
   case PDL_S: size=sizeof(PDL_Short); break;
   case PDL_US: size=sizeof(PDL_Ushort); break;
   case PDL_L: size=sizeof(PDL_Long); break;
   case PDL_IND: size=sizeof(PDL_Indx); break;
   case PDL_LL: size=sizeof(PDL_LongLong); break;
   case PDL_F: size=sizeof(PDL_Float); break;
   case PDL_D: size=sizeof(PDL_Double); break;
  }
  gsl_ran_choose(rng,(double *) out->data, m, (double *) in->data,n,size);

');


pp_core_importList(' qw/ zeroes long barf  /');  # import just a named list to our namespace, so we don't get warning
				     # messages like 'warning 'min' redefined at line ...'

pp_export_nothing;  # set to not export anything. (This is a OO package, it doesn't need to export any methods.)

pp_add_boot('gsl_set_error_handler_off();
');


pp_done();