Haskell – 优化微分方程求解器

.

> C代码使用GCC 4.5.4和-O3编译.它运行在1.166秒.

> Haskell代码使用GHC 7.4.1和-O3编译.它运行21.3秒.

>如果我使用-O3 -fllvm编译Haskell,它将在4.022秒内运行.

那么,我是否遗漏了一些优化我的Haskell代码的东西？

PS.：我使用了以下参数：1e-8 5.

C代码：

#include <stdio.h>
double p, v, a, t;

double func(double t) {
  return t * t;
}

void euler(double dt) {
  double nt = t + dt;
  double na = func(nt);
  double nv = v + na * dt;
  double np = p + nv * dt;

  p = np;
  v = nv;
  a = na;
  t = nt;
}

int main(int argc, char ** argv) {
  double dt, limit;
  sscanf(argv[1], "%lf", &dt);
  sscanf(argv[2], "%lf", &limit);
  p = 0.0;
  v = 0.0;
  a = 0.0;
  t = 0.0;

  while(t < limit) euler(dt);
  printf("%f %f %f %f\n", p, v, a, t);
  return 0;
}

Haskell代码：

import System.Environment (getArgs)

data EulerState = EulerState !Double !Double !Double !Double deriving(Show)
type EulerFunction = Double -> Double

main = do
  [dt, l] <- fmap (map read) getArgs
  print $runEuler (EulerState 0 0 0 0) (**2) dt l

runEuler :: EulerState -> EulerFunction -> Double -> Double -> EulerState
runEuler s@(EulerState _ _ _ t) f dt limit = let s' = euler s f dt
                                             in case t `compare` limit of
                                                  LT -> s' `seq` runEuler s' f dt limit
                                                  _ -> s'

euler :: EulerState -> EulerFunction -> Double -> EulerState
euler (EulerState p v a t) f dt = (EulerState p' v' a' t')
    where t' = t + dt
          a' = f t'
          v' = v + a'*dt
          p' = p + v'*dt

应用于runEuler,我得到了很好的提升.

runEuler :: EulerState -> EulerFunction -> Double -> Double -> EulerState
runEuler s f dt limit = go s
  where go s@(EulerState _ _ _ t) = if t < limit then go (euler s f dt) else s

这有助于f内联到循环中(这可能也发生在C版本中),摆脱了大量的开销.