I am using Haskell to make a Verlet integrator to model gravity. The integrator uses the first two positions of the object as seeds and generates the rest after this.
I thought a nice way of making this in Haskell would be to use an infinite list. However, when implemented I find that it runs very slowly for large times (Haskell 1700 time steps: 12 seconds, Python 1700 time steps: < 1 second)
Here is the relevant code for a 1d integrator that has similar performance:
verletStep dt acc xn xn1 = 2*xn1 - xn + (acc xn1)*dt*dt
verlet dt acc x0 x1 = x0 : x1 : next (verlet dt acc x0 x1)
where
next (xn : xs@(xn1:_)) = (verletStep dt acc xn xn1) : next xs
I also tried using zipWith to generate the infinite list but it has similar performance.
Why does this take so long? The garbage collection itself is around 5 seconds. Is there a nice way to make this run faster?
This definition...
verlet dt acc x0 x1 = x0 : x1 : next (verlet dt acc x0 x1)
where
next (xn : xs@(xn1:_)) = (verletStep dt acc xn xn1) : next xs
... leads to verlet dt acc x0 x1 being calculated many times unnecessarily, thus building a lot of unneeded lists. That can be seen by working out a time step by hand:
verlet dt acc x0 x1
x0 : x1 : next (verlet dt acc x0 x1)
x0 : x1 : next (x0 : x1 : next (verlet dt acc x0 x1))
x0 : x1 : (verletStep dt acc x0 x1) : next (x1 : next (verlet dt acc x0 x1))
The solution is to eliminate the unnecessary list-building:
verlet dt acc x0 x1 = x0 : x1 : x2 : drop 2 (verlet dt acc x1 x2)
where
x2 = verletStep dt acc x0 x1
drop 2 removes the first two elements of a list (in this case, x1 and x2, which we have already prepended). verlet is called recursively with the second position, x1, and the newly calculated third one, x2. (Compare with the original definition, in which verlet is called recursively with the same arguments. That should raise suspicion.)