...find the probability density function of distribution (pdf) of tossing 20 coins at a time.

Let me add a pedantic point. When asked to find a PDF (e.g., in a statistical theory class), the requester is likely to mean that the function should be found analytically. That said, it is equally valid to compute it exactly by enumeration (as roboticus did).

Someone looking for a Monte Carlo approach will typically ask you to estimate the PDF.

You can use Bit::Vector, but I suspect that would be a good deal slower.

There wouldn't be any need to do that. With 32-bit ints, a full run would take around 7 hours. With 64-bit ints, it would take around 3 million years.

If you switched to using unpack to count the bits:

    $tailCnt[ unpack '%32b*', pack 'N', $collection ]++;
[download]

Then the full 32-bit comes down to a very reasonable 40 minutes, but that still leaves the full 64-bit requiring several hundred thousand years.

BrowserUk:

Heh, yep, that would be a good deal slower, all right! Of course, if speed gets to be an issue, you can always skip counting all the iterations, and directly compute the binomial expansion:

use strict;
use warnings;
use List::Util qw(reduce);

my $numTosses = 20;
#my $runs = (1<<$numTosses)-1;

my @triangle = (0, 1, 0);
for (1 .. $numTosses) {
    my @newTriangle=(0);
    push @newTriangle, $triangle[$_]+$triangle[$_+1] for 0 .. $#triang
+le-1;
    push @newTriangle, 0;
    @triangle = @newTriangle;
}

print <<EOHDR;
Tails      Count     %
-----   ----------  ------
EOHDR
my $runs = reduce { $a + $b } @triangle;
for (my $i = 0; $i < $numTosses+1; $i++) {
    printf "% 4u    % 10u   %5.2f\n",
            $i, $triangle[$i+1], 100*$triangle[$i+1]/$runs;
}
[download]

It runs a good bit faster:

$ time perl 892293.pl
Name "main::a" used only once: possible typo at 892293.pl line 21.
Name "main::b" used only once: possible typo at 892293.pl line 21.
Tails      Count     %
-----   ----------  ------
   0             1    0.00
   1            20    0.00
   2           190    0.02
   3          1140    0.11
   4          4845    0.46
   5         15504    1.48
   6         38760    3.70
   7         77520    7.39
   8        125970   12.01
   9        167960   16.02
  10        184756   17.62
  11        167960   16.02
  12        125970   12.01
  13         77520    7.39
  14         38760    3.70
  15         15504    1.48
  16          4845    0.46
  17          1140    0.11
  18           190    0.02
  19            20    0.00
  20             1    0.00

real    0m0.034s
user    0m0.024s
sys     0m0.012s
[download]

Even if you use 32 bits:

$ time perl 892293.pl
Name "main::a" used only once: possible typo at 892293.pl line 21.
Name "main::b" used only once: possible typo at 892293.pl line 21.
Tails      Count     %
-----   ----------  ------
   0             1    0.00
   1            32    0.00
   2           496    0.00
   3          4960    0.00
   4         35960    0.00
   5        201376    0.00
   6        906192    0.02
   7       3365856    0.08
   8      10518300    0.24
   9      28048800    0.65
  10      64512240    1.50
  11     129024480    3.00
  12     225792840    5.26
  13     347373600    8.09
  14     471435600   10.98
  15     565722720   13.17
  16     601080390   13.99
  17     565722720   13.17
  18     471435600   10.98
  19     347373600    8.09
  20     225792840    5.26
  21     129024480    3.00
  22      64512240    1.50
  23      28048800    0.65
  24      10518300    0.24
  25       3365856    0.08
  26        906192    0.02
  27        201376    0.00
  28         35960    0.00
  29          4960    0.00
  30           496    0.00
  31            32    0.00
  32             1    0.00

real    0m0.034s
user    0m0.028s
sys     0m0.004s
[download]

...roboticus

When your only tool is a hammer, all problems look like your thumb.

Update: When I use bignum;, it prints the PDF for 64 tosses per run in just over a second, but I haven't got it formatted nicely...

use strict;
use warnings;
use List::Util qw(reduce);
use bignum;

my $numTosses = 64;
#my $runs = (1<<$numTosses)-1;

my @triangle = (0, 1, 0);
for (1 .. $numTosses) {
    my @newTriangle=(0);
    push @newTriangle, $triangle[$_]+$triangle[$_+1] for 0 .. $#triang
+le-1;
    push @newTriangle, 0;
    @triangle = @newTriangle;
}


print <<EOHDR;
Tails      Count     %
-----   ----------  ------
EOHDR
my $runs = reduce { $a + $b } @triangle;
for (my $i = 0; $i < $numTosses+1; $i++) {
    print "$i\t$triangle[$i+1]\t",100*$triangle[$i+1]/$runs, "\n";
}
[download]

I guess I'll have to figure out how to make bugnum and printf cooperate better...

Here's a slightly speedier alternative to get Pascal's triangle row for the number of tosses:


sub pascal_tri_row {
    my $r = shift;
    return () if $r < 0;

    my @row = (1) x ($r + 1);

    for my $i (1 .. $r - 1)
    {
        $row[$_] += $row[$_ - 1]
          for reverse 1 .. $i;
    }
    return @row;
}

sub triangle {
    my $numTosses = shift;

    my @triangle = (0, 1, 0);
    for (1 .. $numTosses) {
        my @newTriangle=(0);
        push @newTriangle, $triangle[$_]+$triangle[$_+1] for 0 .. $#tr
+iangle-1;
        push @newTriangle, 0;
        @triangle = @newTriangle;
    }

    return @triangle;
}

use Benchmark qw(cmpthese);

cmpthese -1, {
    robo_tri => sub { triangle(32) },
    repel_tri => sub { 0, pascal_tri_row(32), 0 }, # 0's to match tria
+ngle()'s return
};

__END__

            Rate  robo_tri repel_tri
robo_tri  1128/s        --      -38%
repel_tri 1811/s       60%        --
[download]