An infinite stream of “Pi”

So, after TimToady’s help with my last problem, finishing this is trivial. You just convert the Haskell code without worrying about type safety.

type LFT = (Integer, Integer, Integer, Integer) 
extr :: LFT -> Integer -> Rational 
extr (q,r,s,t) x = ((fromInteger q) * x + (fromInteger r)) / 
                           ((fromInteger s) * x + (fromInteger t)) 
unit :: LFT 
unit = (1,0,0,1) 
comp :: LFT -> LFT -> LFT 
comp (q,r,s,t) (u,v,w,x) = (q*u+r*w,q*v+r*x,s*u+t*w,s*v+t*x) 

becomes

sub extr([$q, $r, $s, $t], $x) { 
    ($q * $x + $r) / ($s * $x + $t); 
}

my $unit = [1, 0, 0, 1];

sub comp([$q,$r,$s,$t], [$u,$v,$w,$x]) {
    [$q * $u + $r * $w, 
     $q * $v + $r * $x, 
     $s * $u + $t * $w, 
     $s * $v + $t * $x];
}

And then the final piece in the puzzle,

pi = stream next safe prod cons init lfts where 
  init = unit 
  lfts = [(k, 4*k+2, 0, 2*k+1) | k<-[1..]] 
  next z = floor (extr z 3) 
  safe z n = (n == floor (extr z 4)) 
  prod z n = comp (10, -10*n, 0, 1) z 
  cons z z’ = comp z z’ 

becomes

sub pi-stream() {
    stream(-> $z { extr($z, 3).floor; },
           -> $z, $n { $n == extr($z, 4).floor; },
           -> $z, $n { comp([10, -10*$n, 0, 1], $z); },
           &comp,
           $unit, 
           (1..*).map({ [$_, 4 * $_ + 2, 0, 2 * $_ + 1] }));
}

It’s a very direct translation.

Does it work?

> my @pi := pi-stream;
> say @pi[^40].join('');
3141592653589793238468163213056056860170

Yay!

Except, according to the Joy of Pi, the first 40 digits of pi are

3.1415926535 8979323846 2643383279 502884197 # pi
3.1415926535 8979323846 8163213056 056860170 # ours

What’s going wrong? I haven’t empirically verified it yet, but I’m pretty sure the issue is Rakudo’s Ints and Rats overflowing. Which means our next post is going to have to dive back into Math::BigInt and Math::FatRat…