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

Here's my idea how to do that: 
#!/usr/bin/perl

use warnings;
use strict;

use Math::Polynomial::Solve qw(quadratic_roots);
use Math::Numbers;
use Data::Dumper;

my $level = 64;
my @powers_of_two = (1);
push (@powers_of_two, $powers_of_two[-1]*2) while $powers_of_two[-1] <
+ $level;

my @numbers = @ARGV;
@numbers = (1..$level-1) unless @numbers;

for my $number (@numbers) {
        print "E($number/$level)";
        for (1) {
                my @divisors = Math::Numbers->new($number)->get_diviso
+rs();
                if (grep({$_ == 2} @divisors)) {
                        # trivial
                        print " = E(" . $number/2 . "/" . $level/2 . "
+)";
                } 
                if ($number < $level/2) {
                        # trivial
                        print " = E($number/" . $level/2 . ") & R";
                }
                if (@divisors > 2) {
                        # not a prime, decompose
                        my ($numerator1) = grep({$_ != 1 && $_ != $num
+ber} @divisors);
                        my $numerator2 = $number/$numerator1;
                        my $denominator1 = 2;
                        $denominator1 *= 2 while ($denominator1 < $num
+erator1);
                        my $denominator2 = $level/$denominator1;
                        if ($denominator2 > $numerator2) {
                                print " = E($numerator1/$denominator1)
+ & E($numerator2/$denominator2)";
                        }
                }
                # prime or failed to decompose nicely, try possible ov
+erlaps
                for (my $overlap = 1; $overlap + $number < 2 * $level;
+ $overlap += 2) {
                        my @o_divisors = sort {$a <=> $b} Math::Number
+s->new($overlap)->get_divisors();
                        for my $numerator1 (@o_divisors) {
                                last if $numerator1 > $overlap/2;
                                my $numerator2 = $overlap/$numerator1;
                                # try to find integer d1, d2 such that
+ n1*n2 = level and n1*d2+n2*d1 = number + overlap
                                my $quad_a = $numerator1;
                                my $quad_b = - ($number + $overlap);
                                my $quad_c = $level * $numerator2;
                                my @roots = quadratic_roots($quad_a, $
+quad_b, $quad_c);
                                #print "overlap $overlap, n1 $numerato
+r1, roots " . join(',', @roots) . "\n";
                                for my $root (@roots) {
                                        next if ref $root; # complex r
+oot
                                        $root += 0;
                                        next unless grep({$root == $_}
+ @powers_of_two);
                                        my ($denominator1, $denominato
+r2) = sort {$a <=> $b} ($root, $level/$root);
                                        next if $denominator1 < $numer
+ator1;
                                        next if $denominator2 < $numer
+ator2;
                                        print " = E($numerator1/$denom
+inator1) | E($numerator2/$denominator2)";
                                }
                        }
                }
        }
        print "\n";
}

[download]
For a given power of two ($level = 2^n), it iterates through all the targets in the form p = k/2^n and tries to reduce the expression E(k/2^n) to an expression consisting of lower-level expressions: 
E(1/32) = E(1/16) & R
E(2/32) = E(1/16) = E(2/16) & R
E(3/32) = E(3/16) & R
E(4/32) = E(2/16) = E(4/16) & R = E(2/2) & E(2/16)
E(5/32) = E(5/16) & R
E(6/32) = E(3/16) = E(6/16) & R = E(2/2) & E(3/16)
E(7/32) = E(7/16) & R
E(8/32) = E(4/16) = E(8/16) & R = E(2/2) & E(4/16)
E(9/32) = E(9/16) & R = E(3/4) & E(3/8)
E(10/32) = E(5/16) = E(10/16) & R = E(2/2) & E(5/16)
E(11/32) = E(11/16) & R
E(12/32) = E(6/16) = E(12/16) & R = E(2/2) & E(6/16)
E(13/32) = E(13/16) & R
E(14/32) = E(7/16) = E(14/16) & R = E(2/2) & E(7/16)
E(15/32) = E(15/16) & R = E(3/4) & E(5/8)
E(16/32) = E(8/16) = E(2/2) & E(8/16)
E(17/32) = E(1/4) | E(3/8)
E(18/32) = E(9/16) = E(2/2) & E(9/16)
E(19/32) = E(1/2) | E(3/16)
E(20/32) = E(10/16) = E(2/2) & E(10/16)
E(21/32) = E(3/4) & E(7/8) = E(1/2) | E(5/16)
E(22/32) = E(11/16) = E(2/2) & E(11/16)
E(23/32) = E(1/4) | E(5/8) = E(1/2) | E(7/16)
E(24/32) = E(12/16) = E(2/2) & E(12/16)
E(25/32) = E(1/4) | E(3/8) = E(1/2) | E(9/16)
E(26/32) = E(13/16) = E(2/2) & E(13/16)
E(27/32) = E(3/4) | E(3/8) = E(3/4) | E(3/8) = E(1/2) | E(11/16)
E(28/32) = E(14/16) = E(2/2) & E(14/16)
E(29/32) = E(1/4) | E(7/8) = E(1/2) | E(13/16) = E(3/4) | E(5/8)
E(30/32) = E(15/16) = E(2/2) & E(15/16)
E(31/32) = E(1/2) | E(15/16) = E(3/4) | E(7/8)

[download]
I hope it works, don't have more time to check it now. And I don't have a proof that it will always give results, but it seems to :-)

Replies are listed 'Best First'.
Re^6: Boolean math: Fill in the blanks.
by jdalbec (Deacon) on Oct 11, 2008 at 13:53 UTC
    I wonder if you're not making this too complicated. If I'm not mistaken, 
    E(X & R) = E(X)/2
E(X | R) = (E(X)+1)/2

    [download]
    so it should be possible to build the expression from the binary expansion of the desired probability: 
    #! /usr/bin/perl

use strict;
use warnings;

my $bits = shift;
my $den = 2**$bits;
foreach my $num (1 .. $den-1) {
        my $str = sprintf "%0${bits}b", $num;
        $str =~ s/0*$//;
        my $len = length $str;
        $str =~ s/0/R & (/g;
        $str =~ s/1/R | (/g;
        $str .= ")" x $len;
        $str =~ s/ \| \(\)//;
        $str =~ s/\(R\)/R/;
        print "E($num/$den) = $str\n";
}

    [download]
    The output has more parentheses than strictly necessary, but odds are there's a CPAN module that can fix that. 
    E(1/32) = R & (R & (R & (R & R)))
E(2/32) = R & (R & (R & R))
E(3/32) = R & (R & (R & (R | R)))
E(4/32) = R & (R & R)
E(5/32) = R & (R & (R | (R & R)))
E(6/32) = R & (R & (R | R))
E(7/32) = R & (R & (R | (R | R)))
E(8/32) = R & R
E(9/32) = R & (R | (R & (R & R)))
E(10/32) = R & (R | (R & R))
E(11/32) = R & (R | (R & (R | R)))
E(12/32) = R & (R | R)
E(13/32) = R & (R | (R | (R & R)))
E(14/32) = R & (R | (R | R))
E(15/32) = R & (R | (R | (R | R)))
E(16/32) = R
E(17/32) = R | (R & (R & (R & R)))
E(18/32) = R | (R & (R & R))
E(19/32) = R | (R & (R & (R | R)))
E(20/32) = R | (R & R)
E(21/32) = R | (R & (R | (R & R)))
E(22/32) = R | (R & (R | R))
E(23/32) = R | (R & (R | (R | R)))
E(24/32) = R | R
E(25/32) = R | (R | (R & (R & R)))
E(26/32) = R | (R | (R & R))
E(27/32) = R | (R | (R & (R | R)))
E(28/32) = R | (R | R)
E(29/32) = R | (R | (R | (R & R)))
E(30/32) = R | (R | (R | R))
E(31/32) = R | (R | (R | (R | R)))

    [download]
[reply]
[d/l]
[select]

      The output has more parentheses than strictly necessary

      Don't add parens within runs of | and runs of &, just around them. 

      #!/usr/bin/perl

use strict;
use warnings;

my $bits = shift;
my $den = 2**$bits;
for my $num (1 .. $den-1) {
   for ( sprintf "%0${bits}b", $num ) {
      s/10*$/R/;
      $_ .= ')' x s/(?<=1)(?=0)|(?<=0)(?=1)/(/g;
      s/0/R & /g;
      s/1/R | /g;
      print "E($num/$den) = $_\n";
   }
}

      [download]
      E(1/32) = R & R & R & R & R
E(2/32) = R & R & R & R
E(3/32) = R & R & R & (R | R)
E(4/32) = R & R & R
E(5/32) = R & R & (R | (R & R))
E(6/32) = R & R & (R | R)
E(7/32) = R & R & (R | R | R)
E(8/32) = R & R
E(9/32) = R & (R | (R & R & R))
E(10/32) = R & (R | (R & R))
E(11/32) = R & (R | (R & (R | R)))
E(12/32) = R & (R | R)
E(13/32) = R & (R | R | (R & R))
E(14/32) = R & (R | R | R)
E(15/32) = R & (R | R | R | R)
E(16/32) = R
E(17/32) = R | (R & R & R & R)
E(18/32) = R | (R & R & R)
E(19/32) = R | (R & R & (R | R))
E(20/32) = R | (R & R)
E(21/32) = R | (R & (R | (R & R)))
E(22/32) = R | (R & (R | R))
E(23/32) = R | (R & (R | R | R))
E(24/32) = R | R
E(25/32) = R | R | (R & R & R)
E(26/32) = R | R | (R & R)
E(27/32) = R | R | (R & (R | R))
E(28/32) = R | R | R
E(29/32) = R | R | R | (R & R)
E(30/32) = R | R | R | R
E(31/32) = R | R | R | R | R

      [download]

      More parens can be removed since | and & don't have the same precedence, but I think that would make it hard to read.

      Update: I originally had this longer code:

      Read more... (680 Bytes)
[reply]
[d/l]
[select]
      Looks like you're right. The formulas are a nice application of those given by psini, and the pattern is much better visible here. Elegant solution!
[reply]

      This just gets better and better! I was heading in this direction, but hadn't got there yet++

      The really nice thing about this analysis is that I no longer have to produce the entire sequence for any given depth. Once I have chosen the optimum depth (denominator) and fraction (enumerator) to meet the required ratio, and reduced it to its lowest form, I can just construct a sub to produce my initialisers.

      That works out to be very efficient.

      Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
      "Science is about questioning the status quo. Questioning authority".
      In the absence of evidence, opinion is indistinguishable from prejudice.
[reply]

        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 :-(

        Read more... Lump of code... (6 kB)
[reply]
[d/l]
[select]
Re^6: Boolean math: Fill in the blanks.
by BrowserUk (Patriarch) on Oct 10, 2008 at 16:20 UTC

    This is very cool code. I love that it produces the reductions.

    Update: Even if they are sometimes a little ecclectic: 

    E(64/128) = E(32/64) = E(2/2) & E(32/64)

    [download]

    I'm having a little trouble wrapping it into my test harness, but it seem to produce exactly the results I am after. Thank you.

    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.
[reply]
[d/l]
      Even if they are sometimes a little ecclectic

      Yes, originally, the code stopped searching when it found the first reduction. Given the order in which these are searched, it didn't produce this garbage. I only changed this behavior before pasting the code, because it may be interesting to see multiple possibilities (like with E(21/32) = E(3/4) & E(7/8) = E(1/2) | E(5/16)), and I didn't notice the noise.

      To fix it, change
      if ($denominator2 > $numerator2) {
      to
      if ($denominator1 > $numerator1 && $denominator2 > $numerator2) {
      but I'm sure you know that.

[reply]
[d/l]
[select]

In Section Seekers of Perl Wisdom