PowerSeries.hs

-- Power series in Haskell from Doug McIlroy's beautiful "Power
-- Series, Power Serious" functional pearl and his subsequent "Music
-- of Streams" paper.

module Main where

import Ratio (Rational)

infixl 7 .*
default (Integer, Rational, Double)

-- constant series
ps0, x:: Num a => [a]
ps0 = []             -- 0 : ps0
x   = [0,1]          -- 0 : 1 : ps0

-- arithmetic
(.*):: Num a => a->[a]->[a]
c .* (f:fs) = c*f : c.*fs
c .* _      = []

instance Num a => Num [a] where
    negate []     = []
    negate (f:fs) = (negate f) : (negate fs)
    fs     + []     = fs
    []     + gs     = gs
    (f:fs) + (g:gs) = f+g : fs+gs
    []     * _      = []
    _      * []     = []
    (f:fs) * (g:gs) = f*g : (f.*gs + fs*(g:gs))
    fromInteger c = [fromInteger c] -- was: fromInteger c : ps0

instance Fractional a => Fractional [a] where
    recip fs = 1/fs
    []     / []     = error "0 / 0"
    []     / (0:_)  = error "0 / 0:_"
    []     / _      = []
    (0:fs) / (0:gs) = fs/gs
    (f:fs) / (g:gs) = let q = f/g in
                      q : (fs - q.*gs)/(g:gs)

-- functional composition
compose:: Num a => [a]->[a]->[a]
compose (f:fs) (0:gs) = f : gs*(compose fs (0:gs))

revert::Fractional a => [a]->[a]
revert (0:fs) = rs where
    rs = 0 : 1/(compose fs rs)

-- calculus
deriv:: Num a => [a]->[a]
deriv (f:fs) = (deriv1 fs 1) where
    deriv1 (g:gs) n = n*g : (deriv1 gs (n+1))

integral:: Fractional a => [a]->[a]
integral fs = 0 : (int1 fs 1) where
    int1 (g:gs) n = g/n : (int1 gs (n+1))

expx, cosx, sinx :: Fractional a => [a]
expx = 1 + (integral expx)
sinx = integral cosx
cosx = 1 - (integral sinx)

instance Fractional a => Floating [a] where
    sqrt (0:0:fs) = 0 : sqrt fs
    sqrt (1:fs)   = qs where
        qs = 1 + integral((deriv (1:fs))/(2.*qs))

arctanx, tanx :: Fractional a => [a]
arctanx = integral $ 1/(1+x^2)
tanx = revert arctanx

-- Pascal's triangle.  The generating function for the binomial
-- coefficients of order n is (1 + y)^n.  To enumerate the binomial
-- coefficients for all n, we may replace z in the geometric series
-- 1 + z + z 2 + z 3 + · · · by (1 + y) * z . That is, we compose the
-- generating function for the geometric series, 1/(1 - z), with the
-- polynomial 0 + (1 + y) * z:
--           (1/[1,-1]) . [[0], [1,1]]
pascal :: (Fractional a) => [[a]]
pascal = compose (recip [1,-1]) [[0], [1,1]]

-- tests
test1 = sinx - sqrt(1-cosx^2)
test2 = sinx/cosx - revert(integral(1/(1+x^2)))
test3 = tanx - sinx / cosx
test4 = [[1], [1,1], [1,2,1], [1,3,3,1], [1,4,6,4,1]] == take 5 pascal

main = putStrLn $ show $ take 5 pascal