Perl 6 Fibonacci versus Haskell

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.

Advertisements

Series and Sequences

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…