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…