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.

Hesitant Steps Forward

October 1, 2010

So, I just dived in and started trying to implement duration as a role.

role ABC::Duration {
has $.ticks; our multi sub duration-from-parse($top) is export {
ABC::Duration.new(:ticks($top.Int || 1)); } our multi sub duration-from-parse($top, $bottom) is export { if +($top // 0) == 0 && +($bottom // 0) == 0 { ABC::Duration.new(:ticks(1/2)); } else { ABC::Duration.new(:ticks(($top.Int || 1) / ($bottom.Int || 1))); } } our method Str() { given$.ticks {
when 1 { "---"; } # for debugging, should be ""
when 1/2 { "/"; }
when Int { .Str; }
when Rat { .perl; }
die "Duration must be Int or Rat, but it's { .WHAT }";
}
}
}


The corresponding action for an explicit note length is

    method note_length($/) { if$<note_length_denominator> {
make duration-from-parse($<top> ??$<top>[0] !! "", $<note_length_denominator>[0]<bottom>[0]); } else { make duration-from-parse($<top> ?? $<top>[0] !! ""); } }  I messed around a bunch trying to make "e" does Duration work as a type, but eventually gave up and just coded an ABC::Note type: class ABC::Note does ABC::Duration { has$.pitch;
has $.is-tie; method new($pitch, ABC::Duration $duration,$is-tie) {
say :$duration.perl; self.bless(*, :$pitch, :ticks($duration.ticks), :$is-tie);
}
}


with corresponding action

    method mnote($/) { make ABC::Note.new(~$<pitch>,
$<note_length> ??$<note_length>[0].ast !! ABC::Duration.new(1),
$<tie> eq '-'); }  So… that seems to work okay so far. But it does raise some issues for me. 1. I’m finding this $<top> ?? $<top>[0] !! "" pattern to be very repetitive. Surely there must be a better way to do it? (errr… wait a minute, could that be just $<top>[0] // ""?)

2. I don’t mind the proliferation of small classes in my code, that corresponds nicely to what we are modeling. But I am starting to mind the corresponding proliferation of small source files. Is there a better way to organize things?

… and that’s all I can remember at the moment. I think I’ll wander off and have breakfast.

Durations vs Roles

September 30, 2010

So, I patched in some dodgy but working for this example code to handle the duration of notes. And then I started to think about how to do it properly.

About half of the musical elements in an ABC tune have a duration. So my first thought was, oooo, define an ABC::Duration role! I was excited, because this seemed like the first really good example I’d seen in my own work of where a role was appropriate and seemed superior to class inheritance.

So I started out defining a duration role with a Real attribute for its length, and two new methods, one which took a Real, and one which took the parsed bits of the note_length regex as strings. But then I started thinking about what that meant. If you define a class which does a role, how do you call that role’s new method? If I’m building up the AST, I’m going to want to create a duration before I have an element to attach it to — how do you build an element around an existing duration? And half of the objects which have durations are perhaps not best represented using a stored duration value — for instance, a triplet element’s duration is the total duration of its three wrapped notes, times 2/3 — making the stored duration value redundant.

So, this is one of those posts which I make before I’ve figured out what the answer is. I’ve verified that you can create a standalone role object in Rakudo:

> role Duration { has $.duration; }; my$a = Duration.new(:duration(10)); say :$a.perl "a" => .new(duration => 10)  So that’s one piece of the second point addressed. (Though surely that .perl is wrong?) But I’m not seeing how to do the rest of it at the moment. Any ideas out there? First Sheet Music! September 19, 2010 There were a number of comments on #perl6 on how to improve the code from my last post, but I’m ignoring them all today and forging ahead to get an actual ABC to Lilypond translator up and running. Here’s the core of the code: sub StemToLilypond(Context$context, $stem) { my$match = ABC::Grammar.parse($stem, :rule<mnote>); # bad in at least two ways.... my$pitch = ~$match<pitch>; my$length = ~$match<note_length>; print " { %note-map{$pitch} }{ %cheat-length-map{$length} } "; } sub BodyToLilypond(Context$context, @elements) {
say "\{";

for @elements -> $element { given$element.key {
when "stem" { StemToLilypond($context,$element.value); }
when "barline" { say " |"; }
}
}

say "\}";
}


This code is wrong for a host of reasons. But the cool thing is, for simple tunes, it very nearly works. Basically, we go through the list of musical elements in the tune. When we see a barline, we output a bar line to Lilypond. When we see a stem, we parse it again (ugh) assuming it is the simple form a a stem, then map the pitch portion to the equivalent Lilypond note (ignoring key signature and accidentals), and the rhythm portion to the equivalent Lilypond duration (assuming that the ABC is notated in 8ths, and ignoring almost all potential complexity). This crude approximation almost suffices to properly notate “The Star of Rakudo”!

Here’s the PDF Lilypond outputs when feed the output of this script. At first glance this looks pretty good, but actually it’s got one major issue (as well as a host of minor ones). I think it might be obvious even to those who don’t read sheet music if it is compared to this proper PDF: the notes are all correct, but they are shifted one eighth note off in the measures! This is because Lilypond apparently treats bar lines in the source code as a mere debugging aid — if they don’t agree with what it thinks they should be, it will issue a warning and then go with what it thinks rather than what you said. I guess that behavior might make sense in big orchestral scores, but it is just annoying at this level.

So, what needs to be done to make this tune work:
1) Detect pick up bars and send the \partial command to Lilypond so that it keeps the bar lines in the correct places.

2) Preprocess the music looking for repeats, so that the repeated section can be properly marked for Lilypond.

3) Handle the key signature correctly.

4) Handle the ~ gracing, which in should generate a turn symbol over the note. (Used here to indicate an Irish-style roll.)

It appears to me that these are straightforward programming problems, so I’m going to forge ahead and see what I can do…