我正在Haskell中实现N-Body模拟。 https://github.com/thorlucas/N-Body-Simulation
现在，每个粒子计算其力，然后对彼此粒子进行加速度。换句话说，力的O(n²)计算。如果我要计算每个组合一次，我可以将其减少到 O(n 选择 2)。
let combs = [(a, b) | (a:bs) <- tails ps, b <- bs ]
force = map (comb -> gravitate (fst comb) (snd comb)) combs
但是我不知道如何在不使用状态的情况下将这些应用于粒子。在上面的示例中，ps[Particle]其中
data Particle = Particle Mass Pos Vel Acc deriving (Eq, Show)
从理论上讲，在有状态语言中，我只需循环组合，根据每个a和b的力计算相关加速度，然后在我这样做时更新每个Particleps加速度。
我想过做一些像foldr f ps combs.起始累加器将是当前ps，f是某个函数，它接受每个comb并在ps中更新相关Particle，并返回该累加器。对于如此简单的过程，这似乎确实需要大量内存并且非常复杂。
有什么想法吗？
从 https://github.com/thorlucas/N-Body-Simulation 获取代码
updateParticle :: Model -> Particle -> Particle
updateParticle ps p@(Particle m pos vel acc) =
let accs = map (gravitate p) ps
acc' = foldr ((accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
vel' = (fst vel + fst acc, snd vel + snd acc)
pos' = (fst pos + fst vel, snd pos + snd vel)
in Particle m pos' vel' acc'
step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) ps
并对其进行修改，以便在矩阵(好吧，列表列表...)中分别更新每个粒子的加速度，我们得到......
updateParticle :: Model -> (Particle, [Acc]) -> Particle
updateParticle ps (p@(Particle m pos vel acc), accs) =
let acc' = foldr ((accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
vel' = (fst vel + fst acc, snd vel + snd acc)
pos' = (fst pos + fst vel, snd pos + snd vel)
in Particle m pos' vel' acc'
step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) $ zip ps accsMatrix where
accsMatrix = [map (gravitate p) ps | p <- ps]
。所以问题本质上是如何在计算accsMatrix时减少对gravitate的调用次数，利用gravitate a b=-1 * gravitate b a的事实。
如果我们要打印出accsMatrix，它看起来像...
[[( 0.0, 0.0), ( 1.0, 2.3), (-1.0, 0.0), ...
[[(-1.0, -2.3), ( 0.0, 0.0), (-1.2, 5.3), ...
[[( 1.0, 0.0), ( 1.2, -5.3), ( 0.0, 0.0), ...
...
。所以我们看到了accsMatrix !! i !! j == -1 * accsMatrix !! j !! i.
因此，要使用上述事实，我们需要访问一些索引。首先，我们索引外部列表...
accsMatrix = [map (gravitate p) ps | (i,p) <- zip [0..] ps]
。并将内部列表替换为列表理解...
accsMatrix = [[ gravitate p p' | p' <- ps] | (i,p) <- zip [0..] ps]
。通过zip获取更多可用的索引...
accsMatrix = [[ gravitate p p' | (j, p') <- zip [0..] ps] | (i,p) <- zip [0..] ps]
。然后，关键是，让accsMatrix依赖自己一半的矩阵......
accsMatrix = [[ if i == j then 0 else if i < j then gravitate p p' else -1 * accsMatrix !! j !! i | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]
我们也可以把它分开一点，如下所示...
accsMatrix = [[ accs (j, p') (i, p) | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]
accs (j, p') (i, p)
| i == j = 0
| i < j = gravitate p p'
| otherwise = -1 * accsMatrix !! j !! i
。或使用map避免列表推导式
accsMatrix = map (flip map indexedPs) $ map accs indexedPs
indexedPs = zip [0..] ps
accs (i, p) (j, p')
| i == j = 0
| i < j = gravitate p p'
| otherwise = -1 * accsMatrix !! j !! i
。或者通过使用列表 monad...
accsMatrix = map accs indexedPs >>= (:[]) . flip map indexedPs
。虽然(对我来说)很难看到这些发生了什么。
这种列表列表方法可能存在一些可怕的性能问题，尤其是使用!!，以及由于遍历而仍在运行 O(n²)操作的事实，以及 O (n ·(n – 1))≡ O (n²) 如@leftaroundabout所述，但每次迭代应调用gravitaten * (n-1) / 2次。