[abfsttwf] Optimizing "iterate"

The logistic map becomes interestingly chaotic for multipliers greater than 3.56995.  However, for our purposes, it is just an example of an endomorphic function.

logistic_chaos :: Double -> Double;
logistic_chaos x = 3.57 * x * (1-x);

As an example of a computationally expensive task, we seek the value of the billionth iteration of this function initialized at 0.5.  This post will explore various implementations of iterateN.  Compiler was ghc 7.6.3 on amd64 with compilation flag -O2.

result :: Double;
result = iterateN 1000000000 logistic_chaos 0.5;

iterateN :: Int -> (a -> a) -> a -> a;

The "model" implementation fails due to a space leak:

iterateN_1 n f start = (iterate f start) !! n;

Hopefully, someday the GHC strictness analyzer will be able to prevent the space leak. It seems unfortunate that the obvious implementation of a simple and common task fails so badly.

Making the evaluated function strict still fails with a space leak:

{-# LANGUAGE BangPatterns #-}
strict_chaos :: Double -> Double;
strict_chaos !x = 3.57 * x * (1-x);

We go back to the non-strict version. This one also fails due to a space leak:

iterateN_2 n f = last . (take (1+n)) . (iterate f);

We reimplement the Prelude function "iterate" by inserting "seq".  There were two ways of doing it, so I tried both.  Putting the "seq" in the tail was slower, completing in 17.0 seconds.

iterateN_3 n f start = let {
seq_iterate :: (a -> a) -> a -> [a];
seq_iterate f start = start : (seq start (seq_iterate f (f start)))
} in (seq_iterate f start) !! n;

Putting the "seq" around the whole list was faster, completing in 14.5 seconds:

iterateN_4 n f start = let {
seq_iterate :: (a -> a) -> a -> [a];
seq_iterate f start = seq start (start : (seq_iterate f (f start)))
} in (seq_iterate f start) !! n;

Not passing "f" into the inner function but instead relying on scope cut the runtime to 10.9 and 10.8 seconds. However, this is an undesirable solution as we would like to make seq_iterate a global library function.

{-# LANGUAGE ScopedTypeVariables #-}
iterateN_5 :: forall a . Int -> (a -> a) -> a -> a;
iterateN_5 n f start = let {
seq_iterate :: a -> [a];
seq_iterate start = start : (seq start (seq_iterate (f start)))
} in (seq_iterate start) !! n;

iterateN_6 :: forall a . Int -> (a -> a) -> a -> a;
iterateN_6 n f start = let {
seq_iterate :: a -> [a];
seq_iterate start = seq start (start : (seq_iterate (f start)))
} in (seq_iterate start) !! n;

Manually maintaining a counter runs in 9.7 seconds:

iterateN_7 n f start = case (compare n 0) of {
  LT -> error "cannot iterate less than zero";
  EQ -> start;
  GT -> seq start (iterateN_7 (pred n) f (f start))
};

Using an Integer instead of an Int counter cost an additional 1.14*10^-8 seconds per iteration, so this ran in 21.1 seconds.

iterateN_Integer :: Integer -> (a -> a) -> a -> a;
iterateN_Integer n f start = case (compare n 0) of {
  LT -> error "cannot iterate less than zero";
  EQ -> start;
  GT -> seq start (iterateN_Integer (pred n) f (f start))
};

A pure C++ version (gcc version 4.6.3 compiled with -O2) runs in 4.1 seconds:

#include <iostream>
int main(){
  double x=0.5;
  for(int i=0;i<1000000000;++i){
    x=3.57*x*(1-x);
  }
  std::cout<<x<<std::endl;
}

Complete Haskell source code.