import qualified Data.Stream as S
import Data.Stream (Stream(..))
-- Approximates c + integral f(t) dt for t=0 to infinity using simple Euler steps
-- Both the function and the steps are provided as streams
eulerIntegral :: Double -> Stream Double -> Stream Double -> Stream Double
eulerIntegral c ft dt = S.scan' (+) c $ S.zipWith (*) ft dt
-- Since scan' is strict it doesn't behave well regarding fixed points,
-- e.g. the computation of the approximation of e blocks
e = es S.!! (10^6)
where
es = eulerIntegral 1 es dts
dts = S.repeat 1e-6
-- Replacing S.scan' with myscan' fixes it
myscan' f b ~(a `Cons` as) = b `seq` (b `Cons` myscan' f (f b a) as)