Archive for the ‘Uncategorized’ Category

A Word for the Slow

January 20, 2011

So, my solution for Masak’s p1 has the distinction of being by far the least efficient working solution. Which is a shame, because I think this one was my favorite solution of the contest. It may be slow, and I’d never advise anyone to use it for anything practical (given the other, much more efficient solutions), but in my opinion it’s a charming piece of code.

The key organizational piece for my solution is the Ordering class. (BTW, I didn’t like that name at all, it’s probably my least favorite part of my solution. My theory is Masak was too awestruck by my inefficiency to quibble about the name.) I was striking about for how to represent a series of matrix multiplications, and hit on the idea of using a very simple stack-based language to do it. The language has two operations: an Int represents putting the matrix with that index on the stack. The string "*" represents popping the top two matrices on the stack, multiplying them, and pushing the result back on the stack. Here’s the code for making that happen while tracking the total number of multiplications:

    method calculate-multiplications(@matrices) {
        my @stack;
        my $total-multiplications = 0;
        for @.ops {
            when "*" { 
                my $a = @stack.pop; 
                my $b = @stack.pop;
                my ($multiplications, $matrix) = multiply($a, $b);
                $total-multiplications += $multiplications;
                @stack.push($matrix);
            }
            when Int {
                @stack.push(@matrices[$_]);
            }
        }

        $total-multiplications;
    }

I’m an old Forth programmer from way back, and I can’t begin to say how much I love how easy p6 makes it to implement a simple stack machine!

Getting the string version of this is equally easy:

    method Str() {
        my @stack;
        for @.ops {
            when "*" { 
                my $a = @stack.pop; 
                my $b = @stack.pop;
                @stack.push("($a$b)");
            }
            when Int {
                @stack.push("A{$_ + 1}");
            }
        }

        @stack[*-1];
    }

This time instead of a stack of Pairs (for the matrix size), we have a stack of Str representing each sub-matrix’s name. At the end we pop the last thing on the stack, and it’s the string representing the entire multiplication. And by making this Ordering.Str, any time you print an Ordering you get this nice string form — handy both for the desired output of the program and for debugging.

I won’t comment on the guts of the generate-orderings function, which is heavily borrowed from HOP via List::Utils. Just note that given the number of matrices, it lazily generates all the possible permutations — both the source of my code’s elegance and its extreme inefficiency.

Oonce you’ve got the array @matrices set up, calculating and reporting the best ordering (very slowly!) is as simple as

say generate-orderings(+@matrices).min(*.calculate-multiplications(@matrices));

(Note that I split the main body of this off into a function in the actual code, to make it easier to test internally.)

So clearly, I badly need to study more dynamic programming. But at the same time, I think there may be useful bits in my code that can be put to better use somewhere else.

A Piece of Pi

January 19, 2011

Between masak’s Perl 6 contest and the Perl 6 Advent calendar, I haven’t done a lot of blogging here lately. Apologies!

So, a discussion on #perl6 the other day got me thinking about another possible interesting use of the Real role. Here’s the basic class:

class TimesPi does Real {
    has Real $.x;
    
    method new($x) { self.bless(*, :$x); }
    method Bridge() { $.x * pi; }
    method Str() { $.x ~ "π"; }
    method perl() { "({ $.x })π"; }
}

sub postfix:<π>($x) {
    TimesPi.new($x);
}

So, that’s simple enough, right? TimesPi stores a Real number $.x internally, and the value it represents is that number times pi. There’s a postfix operator π to make it really easy to construct these numbers. Because we’ve defined a Bridge method, this class has access to all the normal methods and operators of Real. Still, as presented above it is pretty useless, but defining some operators hints at a useful purposes for this class.

multi sub infix:<+>(TimesPi $lhs, TimesPi $rhs) {
    TimesPi.new($lhs.x + $rhs.x);
}

multi sub infix:<->(TimesPi $lhs, TimesPi $rhs) {
    TimesPi.new($lhs.x - $rhs.x);
}

multi sub prefix:<->(TimesPi $n) {
    TimesPi.new(- $n.x);
}

multi sub infix:<*>(TimesPi $lhs, Real $rhs) {
    TimesPi.new($lhs.x * $rhs);
}

multi sub infix:<*>(Real $lhs, TimesPi $rhs) {
    TimesPi.new($lhs * $rhs.x);
}

multi sub infix:<*>(TimesPi $lhs, TimesPi $rhs) {
    $lhs.Bridge * $rhs.Bridge;
}

multi sub infix:</>(TimesPi $lhs, Real $rhs) {
    TimesPi.new($lhs.x / $rhs);
}

multi sub infix:</>(TimesPi $lhs, TimesPi $rhs) {
    $lhs.x / $rhs.x;
}

With these operators in place, basic arithmetic involving TimesPi numbers will stay in the TimesPi class when appropriate. For instance, if you add two TimesPi numbers, the result will be a TimesPi. The cool thing about this is that it is as exact $.x values allow, rather than forcing everything to be a floating point calculation of limited accuracy.

We can even take things a step further, using this to perform exact trig calculations:

    multi method sin(TimesPi $x: $base = Radians) {
        return $x.Bridge.sin($base) unless $base == Radians;
        given $x.x {
            when Int { 0; }
            when Rat {
                given $x.x.denominator {
                    when 1 { 0; }
                    when 2 { ($x.x.numerator - 1) %% 4 ?? 1 !! -1 }
                    # could have a case for 6 as well...
                    $x.Bridge.sin; 
                }
            }
            default {
                $x.Bridge.sin;
            }
        }
    }

This checks for cases where we know the exact value of the result, and returns that if it can, otherwise falling back to the standard Real.sin method.

Of course, just when I was feeling like I might be on to something here, I realized that $.x was just the number of degrees in the angle divided by 180. Sigh.

Perl 6 Fibonacci versus Haskell

December 29, 2010

There’s been some discussion on reddit today about whether

my @fib := 1, 1, *+* ...^ * >= 100;

is unreadable gibberish or not, with the following Haskell suggested as an easier-to-understand version.

fib = 1 : 1 : zipWith (+) fib (tail fib)

(I’ve “corrected” both versions so they start the sequence with 1, 1.)

The first thing to observe here is that this are not the same at all! The Perl 6 version is the Fibonacci numbers less than 100, while the Haskell version lazily generates the entire infinite sequence. If we simplify the Perl 6 to also be the (lazy) infinite Fibonacci sequence, we get the noticeably simpler

my @fib := 1, 1, *+* ... *;

To my (admittedly used to Perl 6) eye, this sequence is about as clean and straightforward as it is possible to get. We have the first two elements of the sequence:

1, 1

We have the operation to apply repeatedly to get the further elements of the sequence:

*+*

And we are told the sequence will go on forever:

... *

The *+* construct may be unfamiliar to people who aren’t Perl 6 programmers, but I hardly think it is more conceptually difficult than referring to two recursive copies of the sequence you are building, as the Haskell version does. Instead, it directly represents the simple understanding of how to get the next element in the Fibonacci sequence in source code.

Of course, this being Perl, there is more than one way to do it. Here’s a direct translation of the Haskell version into idiomatic Perl 6:

my @fib := 1, 1, (@fib Z+ @fib[1..*]);

Well, allowing the use of operators and metaoperators, that is, as zipWith (+) becomes Z+ and tail fib becomes @fib[1..*]. To the best of my knowledge no current Perl 6 implementation actually supports this. I’d be surprised if any Perl 6 programmer would prefer this version, but it is out there.

If you’re insistent on writing function calls instead of operators, you could also say it

my @fib := 1, 1, zipwith(&[+], @fib, tail(@fib));

presuming you write a tail function, but that’s an easy one-liner.

Series and Sequences

December 29, 2010

rokoteko asked on #perl6 about using sequences to calculate arbitrarily accurate values. I’m not sure why he thought I was the expert on this, but I do have some ideas, and thought they should be blogged for future reference.

Say we want to calculate pi using a sequence. The Taylor series for atan is

atan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9...

We can represent those terms easily as a lazy list in Perl 6:

sub atan-taylor($x) { 
    (1..*).map(1 / (* * 2 - 1)) Z* ($x, * * $x * -$x ... *) 
};

Note that we don’t try to do this exclusively as a sequence; getting 1, 1/3, 1/5, 1/7... is much easier using map, and then we mix that with the sequence of powers of $x using Z*.

So, one nice thing about the sequence operator is that you can easily use it to get all the Rats in the terms of a Taylor series, because once the terms get small enough, they will switch to Nums. So we can say

> (atan-taylor(1/5) ...^ Num).perl
(1/5, -1/375, 1/15625, -1/546875, 1/17578125, -1/537109375)

This doesn’t help us get Rat approximation to the series, however, because summing those values results in a Num:

> ([+] (atan-taylor(1/5) ...^ Num)).WHAT
Num()

However, we can use the same idea with the triangle-reduce operator to easily get a version that does work:

> (([\+] atan-taylor(1/5)) ...^ Num).perl
(1/5, 74/375, 9253/46875, 323852/1640625, 24288907/123046875)

We’re mostly interested in the last element there, which is easily split off from the rest:

> (([\+] atan-taylor(1/5)) ...^ Num)[*-1].perl
24288907/123046875

So, having laid that groundwork, how do we calculate pi?

My first answer was the classic pi = 4*atan(1). Unfortunately, as sorear++ pointed out, it is terrible for these purposes. Why?

atan(1) = 1 - 1/3 + 1/5 - 1/7 + 1/9...

A little thought there shows that if you want to get to the denominator of 537109375 that took six terms for atan(1/5), it will take 268,554,687 terms for atan(1). Yeah, that’s not very practical.

Luckily, the above-linked web page has a much better formula to use:

atan(1) = 4 * atan(1/5) - atan(1/239)

It takes a little playing around, but it’s reasonably clean to implement this in p6:

use List::Utils;

sub atan-taylor($x) { 
    (1..*).map(1 / (* * 2 - 1)) Z* ($x, * * $x * -$x ... *) 
};

my @fifth-times-four := atan-taylor(1/5) Z* (4, 4 ... *);
my @neg-two-three-ninth := atan-taylor(1/239) Z* (-1, -1 ... *);
my @terms := sorted-merge(@fifth-times-four, @neg-two-three-ninth, *.abs R<=> *.abs);
my @pi-over4 = ([\+] @terms) ...^ Num;
say (@pi-over4[*-1] * 4).perl;
say (@pi-over4[*-1] * 4);

The results are

1231847548/392109375
3.14159167451684

I use sorted-merge from List::Utils to merge the two sequences of Taylor series terms into one strictly decreasing (in magnitude) sequence. That, in turn, makes it easy to use the triangle-reduce metaoperator to stop summing the terms when they’ve gotten so small they are no longer representable by a Rat.

What is this good for? Well right now, not much. Sure, we’re gotten a fairly accurate Rat version of pi, but we could have gotten that more quickly and accurately by just saying pi.Rat(1e-15).

But once we have working FatRats, this approach will let us get arbitrarily accurate rational approximations to pi. Indeed, it suggests we could have slow but very accurate ways of calculating all sorts of transcendental functions…

Moving Right Along

November 2, 2010

Of the six goals in my last post, four of them now work well. Here’s the latest PDF my code has produced. I’ve added three new tunes to it: another of my tunes, “Tali Foster’s”, a quick snippet from our ceili band’s notes from a few years back, and a double from the repertoire of Rufus Guinchard, “Sydney Pittman’s”. Together they demonstrate time signatures (6/8 and a single bar of 9/8), key changes, broken rhythms, and of course more than one tune on a page. (Note that 1st and 2nd endings are still unimplemented, and look very wrong in “Tali Foster’s”.)

I have to say that this PDF really impressed me with Lilypond’s output. It’s hard to put my finger on it, but something about having all four tunes together like that on the page looks really good, IMO. I’m getting excited about the prospect of being able to produce entire books of tunes with this tool.

In terms of Perl 6, this last nine days of work has been very straightforward. I did whine enough about the undeclared variable error message bug to convince pmichaud++ to fix it. I’d worried a lot about nested levels of Contexts — the key, time signature, etc — but I realized as I looked at it there are no nested levels, all changes affect the rest of the piece (until there is another change). I refactored to change several of the subs into methods, and made the context an attribute of their class. I’ve set things up so that Context objects are read-only, and you just create a new one when you need to change the context of the tune. So far this seems to work well.

I guess at this point I really need to implement 1st and 2nd endings, and then push on by throwing more tunes at it to see where they don’t work.

Update: Just in case you wanted to see what “Tali Foster’s” was supposed to look like with proper repeats, here’s the latest PDF. That is to say, first and second endings now work, at least in this one case. I’ve also realized that while key signatures work at the moment, accidentals don’t properly modify the key signature until the end of the bar — so that’s another thing that needs doing.

ABC Update

October 23, 2010

I’ve been slowly poking away at the ABC module without reporting on it here, as the changes have been pretty straightforward. All of the goals at the end of my First Sheet Music post have now been achieved. Here’s the latest version of “The Star of Rakudo” PDF. If you compare it to the output of my old sheet music generator, you’ll see every musical element except the time signature is now correctly rendered by the combination of the Perl 6 ABC code and Lilypond. (In addition, I’ve also added support for rests and triplets, which are now found in the module’s version of the Star of Rakudo ABC file as an example.)

Where to go from here?

1) I guess I ought to fix the time signature thing. That should be trivial.

2) Support for ABC files whose base note duration is something other than an eighth note. (Right now we’ve just hardcoded things to assume eighth notes are the base unit of time.)

3) Broken rhythms.

4) In-line time signature and key signature changes.

5) Handling more than one ABC tune at a time in the input.

I don’t see any major challenges here, other than finding the time to work on this project!

Update: Later in the same day, I’ve already got #1 and #5 working, but I just realized I left out one important (and possibly tricky) one:

6) Handling first and second endings.

Fibonacci and Primes

October 20, 2010

The middle challenge was to find the first prime Fibonacci number greater than 227000, add one to it, and then sum the prime numbers which were its factors. Here’s my first implementation:

sub is-prime($a) {
    return Bool::True if $a == 2;
    ! [||] $a <<%%<< (2, 3, *+2 ... * > $a.sqrt);
}

my @fib := (1, 1, *+* ... *);

my $cutoff = 227000;
my $least-prime = 0;
for @fib -> $f {
    next if $f <= $cutoff;
    next unless is-prime($f);
    $least-prime = $f;
    last;
}

my $x = $least-prime + 1;
say [+] (2, 3, *+2 ... * > $x.sqrt).grep({ $x %% $^a && is-prime($a) });

Despite what seems like an obvious inefficiency (approximating the prime numbers with the odd numbers), this is pretty snappy, executing in 12.5 seconds.

I was planning to go on and talk about my new Math::Prime module here, but looking at this code, I think it can be expressed rather more nicely with a tweak or two here. Let’s see.


sub is-prime($a) {
    return Bool::True if $a == 2;
    ! [||] $a <<%%<< (2, 3, *+2 ... * > $a.sqrt);
}

my @fib := (1, 1, *+* ... *);
my $cutoff = 227000;
my $least-prime = @fib.first({ $_ > $cutoff && is-prime($_) });
my $x = $least-prime + 1;
say [+] (2, 3, *+2 ... * > $x.sqrt).grep({ $x %% $^a && is-prime($a) });

So that’s what the first method is good for!

I did indeed write Math::Prime just so I could use it here. It’s not a huge change from the previous version, really:

use Math::Prime;

my @fib := (1, 1, *+* ... *);
my $cutoff = 227000;
my $least-prime = @fib.first({ $_ > $cutoff && is-prime($_) });
my $x = $least-prime + 1;
say [+] (primes() ... * > $x.sqrt).grep({ $x %% $^a });

Unfortunately, Math::Prime isn’t optimized yet, and so this version, while a bit nicer, is actually slower than the previous version.

Update: With Moritz’s help I did some optimizations to Math::Prime, and the script using it is now very significantly faster than the others.

Summing Subsets

October 9, 2010

So, the third challenge was to count the number of subsets of a set of numbers such that the largest number in the subset is the sum of the rest of the numbers in the subset. My first attempt was very straightforward: create all the subsets and check to see if they have the desired property:

my @a = 3, 4, 9, 14, 15, 19, 28, 37, 47, 50, 54, 56, 59, 61, 70, 73, 78, 81, 92,
 95, 97, 99;

my $hits = 0;

for 1..(2 ** +@a) -> $key {
    my @b = gather for 0..+@a -> $i { take @a[$i] if $key +& (2 ** $i); }
    my $test = @b.pop;
    next if $test != [+] @b;
    $hits++;
    say (@b, $test).join(' ');
}

say "$hits hits";

I think this works correctly, but it will be a long time before we know — as I type this, it’s been running for ten hours on my fastest computer, and I don’t anticipate it finishing any time soon.

My second attempt relies on recursion, the fact the list is sorted, and skipping fruitless branches to get the job done much faster — 47 seconds, to be precise.

sub CountSumSets($target, @a) {
    my $sets = 0;
    for ^ +@a -> $i {
        if @a[$i] < $target {
            $sets += CountSumSets($target - @a[$i], @a[($i + 1) .. (@a - 1)]);
        } elsif @a[$i] == $target {
            $sets += 1;
        }
    }
    $sets;
}


my @a = 3, 4, 9, 14, 15, 19, 28, 37, 47, 50, 54, 56, 59, 61, 70, 73, 78, 81, 92, 95, 97, 99;
@a .= reverse;

my $hits = 0;

for ^ +@a -> $i {
    $hits += CountSumSets(@a[$i], @a[($i + 1) .. (@a - 1)]);
}

say $hits;

Longest palindrome

October 9, 2010

Via Hacker News, I found the Greplin Programming Challenge, and couldn’t resist trying it in Perl 6. If you’re the sort of person who enjoys that sort of thing, I highly encourage you to stop reading here and go try it yourself!

I’m going to blog about all three challenges one-by-one. I suppose the third will be the most interesting, if for no other reason than my current implementation looks like it will take ~32 hours to run, so I’m probably going to need to find a more clever solution.

Basically, the first challenge is to find the longest palindrome in a string without spaces. As stated, I thought the challenge implied it was case-sensitive; obviously it’s easiest enough to add a call to lc to get a case-insensitive version.

I was pretty sure a naive version would bring Rakudo to its knees, so I tried to be slightly clever. My solution is still O(N^2), but N is the number of occurrences of a given letter rather than the full length of the string, so that’s fairly reasonable.

sub longest-palindrome($string) {
    my @c = $string.comb(/./);
    my %hc;
    for ^ +@c -> $i {
        if %hc{@c[$i]} {
            %hc{@c[$i]}.push($i);
        } else {
            %hc{@c[$i]} = [$i];
        }
    }

    my @candidates := gather for %hc.keys.sort -> $c {
        say "Checking $c";

        my $list = %hc{$c};
        say :$list.perl;
        for $list.list -> $i1 {
            for $list.reverse -> $i2 {
                last if $i2 <= $i1;
                my $j1 = $i1;
                my $j2 = $i2;
                my $candidate = Bool::True;
                while ++$j1 < --$j2 {
                    if @c[$j1] ne @c[$j2] {
                        $candidate = Bool::False;
                        last;
                    }
                }
                if $candidate {
                    say @c[$i1..$i2];
                    take @c[$i1..$i2].join('');
                }
            }
        }
    };

    @candidates.sort({$^a.chars <=> $^b.chars}).perl.say;
}

Basically, my notion was to store all the occurrences of each letter in a hash of arrays, then pair up every two occurrences of the same letter and see if they are a palindrome. This probably isn’t the most elegant solution, nor the fastest, but it was fast enough to get solve the challenge problem in a minute or two, and easy enough to code I could do it in under an hour at three in the morning.

Interestingly, I think there might be a way to solve this using regular expressions…

Update: Moritz has a solution which blows this one out of the water, both much faster and more elegant. (The key idea was using the regex engine to find the center of potential palindromes.) I’ll let him tell you about it

Taking a Rest

October 4, 2010

Once I was reasonably happy with ABC::Duration and ABC::Note, I went ahead and added ABC::Rest, which was truly easy thanks to the magic of code reuse. It actually turned out to be quite a struggle to implement, but that was only because I accidentally typed is $.type rather than has $.type and spent a couple hours trying to sort out what the LTA error message meant.

class ABC::Rest does ABC::Duration {
    has $.type;

    method new($type, ABC::Duration $duration) {
        self.bless(*, :$type, :ticks($duration.ticks));
    }

    method Str() {
        $.type ~ self.duration-to-str;
    }
}

There are several sorts of rests in ABC, which is why ABC::Rest has a $type attribute. In practice, I’ve only implemented “z” rests so far, as they are the basic type.

Rests were already in the grammar, and adding an action to support them was dead easy:

    method rest($/) {
        make ABC::Rest.new(~$<rest_type>, $<note_length>.ast);
    }

After a trivial refactor to make the Lilypond duration handling code easily available, it was just took one line to add rests to the Lilypond output:

when “rest” { print ” r{ DurationToLilypond($context, $element.value) } ” }

That’s for the output section of the code. No changes were required to make rests work in the “figure out the overall duration” section of the code, because it looks for elements which do ABC::Duration, and so it automatically found the rests and correctly handled them.

To me, this looks like a solid win for treating duration as a role.


Follow

Get every new post delivered to your Inbox.