Power Series in Haskell

I have come to think that Haskell is one of the nicest vehicles today for doing numerics. As a trivial example, I offer the following code fragment:

factorial 0 = 1
factorial n = n * factorial(n-1)

expCoeffs = map (\n -> 1.0 / fromIntegral (factorial n)) [0..]

powerSeries s epsilon x = sum $ reverse $ takeWhile (>epsilon) $ zipWith (*) s xs
    where xs = map (x**) [0..]

myExp = powerSeries expCoeffs 1e-9

powerSeries takes a sequence of coefficients, a \varepsilon (which is the error bound), and x is the point where the function is to be evaluated. This automatically takes as many terms as is necessary for any given x.

Similar tricks with infinite series I will work charmingly to integrate ordinary differential equations. Note the reverse in the powerSeries. This is necessary for certain series, otherwise the small terms at the end make no change in the large values already there.

Leave a comment