gistfile1.hs

module Main where

import Ratio

stream :: (b->c) -> (b->c->Bool) -> (b->c->b) -> (b->a->b) -> b -> [a] -> [c]
stream next safe prod cons z (x:xs)
  = if safe z y
    then y : stream next safe prod cons (prod z y) (x:xs)
    else stream next safe prod cons (cons z x) xs
      where y = next z

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)

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'

piL = stream next safe prod cons init lfts where
  init = ((0,4,1,0), 1)
  lfts = [(2*i-1, i*i, 1, 0) | i<-[1..]]
  next ((q,r,s,t),i) = floor ((q*x+r) % (s*x+t)) where x=2*i-1
  safe ((q,r,s,t),i) n = (n == floor ((q*x+2*r) % (s*x+2*t)))
                         where x=5*i-2
  prod (z,i) n = (comp (10, -10*n, 0, 1) z, i)
  cons (z,i) z' = (comp z z', i+1)