in reply to Re^7: Boolean math: Fill in the blanks.
in thread Boolean math: Fill in the blanks.

Now that brother jdalbec has shown us the way...

...I observe that the values can be calculated directly from the bit pattern -- right shifting trailing zeros off, and then oring/anding new R's for each '1'/'0', shifting right until have shifted at total of 'n' times, where the denominator is 2**n.

This leads to a minor variation on the theme of constructing a Perl expression to do this. Which takes the bits in reverse order. Taken in this order, the operations want to be executed left to right, so only need to add brackets to cope with the higher precedence of &.

Taken in the forward order brackets are required to ensure that where there is a mix of operations, those at the back end of the expression are done first. This may, or may not, be obvious. (Though within groups of things ored together the order is irrelevant, and similarly groups of things anded together.)

I don't know if this is any simpler, but I wrote it and got it working, so here you are...

FWIW, on my machine (64bit modest Semperon) I benchmarked the code to create values directly from the bit-pattern against an anonymous subroutine built by eval. The subroutine was faster, but only by 75% -- so I guess the running time is dominated by rand. I was using a denominator of 2**33, which mapped a probability of 0.1 to:

  sub { ((((((((R | R) & R & R | R | R) & R & R | R | R) & R & R | R | R)
      & R & R | R | R) & R & R | R | R) & R & R | R | R) & R & R | R | R) & R & R & R }
Easy when you know how :-)

In passing, I discovered that the built-in rand wasn't up to much on my 32-bit machine :-(

use strict ;
use warnings ;

use POSIX qw(ULONG_MAX) ;

my $MB = bits(ULONG_MAX) ;

my $p  = shift(@ARGV) || rand ;    # 1st argument -- probability to tr
+y
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%%  sam
+ple=%d\n", 
                         $count/$sample, 100 * $count/($sample * $B), 
+100 * $p, $sample ;

my $s  = prepare_sub($p, $D, $N) ;    # construct a subroutine equival
+ent
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%%  sam
+ple=%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 denominato
+r
#
# 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 (ef
+fectively) 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 repre
+sentation 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 delim
+it 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 return
+ed 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 ;
} ;

[download]

In Section Seekers of Perl Wisdom