use strict ;
use warnings ;

use POSIX qw(ULONG_MAX) ;

my $MB = bits(ULONG_MAX) ;

my $p  = shift(@ARGV) || rand ;	# 1st argument -- probability to try
my $D  = shift(@ARGV) || '1K' ; # 2nd argument -- denominator to use
my $B  = shift(@ARGV) || 32 ;   # 3rd argument -- bits in random value

if ($D =~ /^(\d+)k$/i) { $D = $1 * 1024 ; } ;
if ($D =~ /^(\d+)m$/i) { $D = $1 * 1024 * 1024 ; } ;
if ($D =~ /^(\d+)g$/i) { $D = $1 * 1024 * 1024 * 1024 ; } ;
if ($B > $MB) { die "Asked for $B bits, but only $MB available" ; } ;

my $N  = 2**$B ;
my $DB = bits($D, 1) ;

my $SR = rand(time) ;

sub R() { int(rand($N)) ; } ;

for my $i (0..32) {
  printf "%4d: %s\n", $i, prepare_sub($i/32, 32, $N) ;
} ;

printf "Probability %.6g, denominator 2**$DB, $B bit numbers\n", $p ;

my $w  = prepare($p, $D) ;		# construct working value

my $sample = 10000 ;
my $count ;

$count = 0 ;
srand $SR ;
for (1..$sample) {
  my $r = make($w, $N) ;
  while ($r) { $count += $r & 1 ; $r >>= 1 ; } ;
} ;

printf "'make': average '1's %.2f/$B or %.5g%% cf required %.5g%%  sample=%d\n", 
                         $count/$sample, 100 * $count/($sample * $B), 100 * $p, $sample ;

my $s  = prepare_sub($p, $D, $N) ;	# construct a subroutine equivalent
my $rs = eval $s ;

$count = 0 ;
srand $SR ;
for (1..$sample) {
  my $r = &$rs ;
  while ($r) { $count += $r & 1 ; $r >>= 1 ; } ;
} ;

printf "'sub' : average '1's %.2f/$B or %.5g%% cf required %.5g%%  sample=%d\n", 
                         $count/$sample, 100 * $count/($sample * $B), 100 * $p, $sample ;
print $s, "\n" ;


use Benchmark qw(timethese) ;

my $bm_count = 20 ;
timethese($bm_count, {
	'make' => sub { srand $SR ; make($w, $N) for (1..$sample) ; },
	'sub ' => sub { srand $SR ; &$rs         for (1..$sample) ; },
});

#=========================================================================================
# prepare_sub: make a subroutine from given probability and denominator
#
# Requires: $p    -- probability: 0..1
#           $D    -- denominator -- 2**n, where $D <= ULONG_MAX + 1
#           $N    -- $N-1 is maximum random value -- in case $p is (effectively) 1
#
# Returns:  ref:sub

sub prepare_sub {
  my ($p, $D, $N) = @_ ;

  my $w = prepare($p, $D) ;
  my $s ;

  if ($w) {
    $s = reverse(sprintf("%b", $w)) ;
    $s =~ s/1$// ;		# Discard the "guard"

    my $b = $s =~ s/(0*1+)0/$1)0/g ;
    $s =~ s/0/ & R/g ;
    $s =~ s/1/ | R/g ;

    $s = '(' x $b .'R'. $s ;
  }
  elsif (defined($w)) {
    $s = '0' ;
  }
  else {
    $s = $N - 1 ;
  } ;

  return "sub { $s }" ;
} ;

#=========================================================================================
# prepare: given a probability and a denominator, return working representation for 'make'
#
# Working representation is:
#
#          0 => probability is effectively 0.
#      undef => probability is effectively 1.
#
#    otherwise: rounded(probability * denominator) + denominator
#               shifted right so that the ls '1' bit falls off 
#
#                 -- the denominator introduces a 'guard' bit to delimit leading '0's.
#                 -- the ls '0's are irrelevant.
#                 -- the ls '1' is implied.
#
# Requires: $p    -- probability: 0..1
#           $D    -- denominator -- 2**n, where $D <= ULONG_MAX + 1
#
# Returns:  working value for 'make'

sub prepare {
  my ($p, $D) = @_ ;

  if (($p > 1) or ($p < 0))    { die "probability $p is >1 or <0" ;        } ;
  if (  $D > (ULONG_MAX + 1) ) { die "denominator $D is > ULONG_MAX + 1" ; } ;
  if ( ($D < (ULONG_MAX + 1))
       && ($D & ($D - 1)))     { die "denominator $D is not some 2**n" ;   } ;

  $p = int($p * $D + 0.5) ;

  if ($p == 0)  { return 0 ;     } ;
  if ($p == $D) { return undef ; } ;

  return int(($D + $p) / (($p ^ ($p - 1)) + 1)) ;
} ;

#=========================================================================================
# make: make a random value according to working representation returned by 'prepare'
#
# Requires: $w    -- working representation
#           $N    -- want random value in 0..$N-1
#
# Returns:  working value for 'make'

sub make {
  my ($w, $N) = @_ ;

  if ($w) {
    my $r = R ;
    while ($w > 1) {
      my $n = R ;
      if ($w & 1) { $r |= $n ; }
             else { $r &= $n ; } ;
      $w >>= 1 ;
    } ;
    return $r ;
  }
  elsif (defined($p)) {
    return 0 ;
  }
  else {
    return int($N - 1) ;
  } ;
} ;

#=========================================================================================
# bits: return number of bits to represent value, or value - 1.

sub bits {
  my ($n, $o) = @_ ;

  if ($o) { $n -= 1 ; } ;
  my $b = 0 ; while ($n) { $b++ ; $n >>= 1 ; } ;

  return $b ;
} ;