Bulk generate beta distribution random number generator in XS

README.md

Calling a XS function many times makes overhead. I recommend to run a for-loop in XS. It makes 155% faster than original implementation. (I want to generate a lot of β dsitribution random numbers)

benchmark.pl

use strict;
use warnings;
use Math::Random qw/random_beta/;
use Benchmark qw/cmpthese/;
use Carp;

cmpthese(
    -1, {
        orig => sub {
            for my $n (1..1000) {
                random_beta_072(1, 0.1, 0.2);
            }
        },
        xs => sub {
            random_beta(1000, 0.1, 0.2);
        },
    }
);


# Math::Random-0.72's original implementation
sub random_beta_072 { # 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 Math::Random::genbet($aa,$bb) unless wantarray();
    my $val;
    my @ans = (0) x $n;
    foreach $val (@ans) { $val = Math::Random::genbet($aa,$bb); }
    return @ans;
}

__END__

Benchmark result on my machine:

           Rate orig   xs
    orig 1001/s   -- -61%
    xs   2549/s 155%   --

bulk-run-random_beta-in-xs.patch

diff --git Random.pm Random.pm
index c268f56..8cfdc17 100755
--- Random.pm
+++ Random.pm
@@ -79,19 +79,6 @@ 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) = @_;
diff --git Random.xs Random.xs
index eb4734c..afe7b50 100644
--- Random.xs
+++ Random.xs
@@ -67,6 +67,24 @@ get_seed()
        OUTPUT:
        RETVAL
 
+void
+random_beta(n, aa, bb)
+	PREINIT:
+	IV i;
+	INPUT:
+	IV n
+	NV aa
+	NV bb
+	PPCODE:
+		if ((aa < 1.0E-37) || (bb < 1.0E-37)) {
+			croak("($aa = %f < 1.0E-37) or ($bb = %f < 1.0E-37)\nin random_beta($n,$aa,$bb)",
+                aa, bb);
+		}
+		for (i=0; i<n; ++i) {
+            mXPUSHn(genbet(aa, bb));
+		}
+        XSRETURN(n);
+
 double
 genbet (aa,bb)
 	INPUT: