log.hs

{-# LANGUAGE FlexibleContexts #-}

-- Version 2.0

import Data.Ratio

(*!) _ 0 = 0
(*!) a b = a*b

(^+) a b = zipWith (+) a b
(^-) a b = zipWith (-) a b

~(a : as) `convolve` (b : bs) = (a *! b) :
    ((map (a *) bs) ^+ (as `convolve` (b : bs)))

compose (f : fs) (0 : gs) = f : (gs `convolve` (compose fs (0 : gs)))

integrate x = 0 : zipWith (/) x (map fromInteger [1..])

sin' = integrate cos'
cos' = let _ : cos'' = map negate (integrate sin') in 1 : cos''

delta (g : gs) h = let g' = delta gs h
                   in (0 : ((1 : h) `convolve` g')) ^+ gs

fsqrt (0 : 1 : fs) =
    let gs = map (/ 2) $ fs ^- (0 : gs `convolve`
                                ((0 : delta gs gs) ^+
                                 ((2 : gs) `convolve` (gs `compose` g))))
        g = 0 : 1 : gs
    in g

type Q = Ratio Integer
type R = Double

-- The function x -> x+x^2
f = 0 : 1 : 1 : repeat 0 :: [Rational]

sqrt' x = 1:rs where rs = map (/2) (xs ^- (rs `convolve` (0:rs)))
                     _ : xs = x

-- The inverse of x -> x+x^2, i.e. x -> (1/2)(-1 + sqrt(1+4x))
f' = fmap (/2) $ sqrt' (1 : 4 : repeat 0) ^- (1 : repeat 0)

church 0 f = id
church n f = church (n-1) f . f

n = 256

g, g' :: [Q]
g = church n fsqrt f ^- (0 : 1 : repeat 0)
g' = church n fsqrt f' ^- (0 : 1 : repeat 0)

main = do
    print "'log' f"
    mapM_ print $ take 20 (map (flip approxRational ((1/10)^20)) $ map ((2^n) *) g :: [Q])
    print "'log' (inverse of f)"
    mapM_ print $ take 20 (map (flip approxRational ((1/10)^20)) $ map ((2^n) *) g' :: [Q])