从
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]
所以问题本质上是如何通过使用gravitate a b = -1 * gravitate b a
的事实,在计算accsMatrix
时减少对gravitate
的调用次数。
如果我们打印出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
...或者通过使用列表单子(list monad)...
accsMatrix = map accs indexedPs >>= (:[]) . flip map indexedPs
虽然对我来说,这些东西更难理解。
这种列表嵌套的方法可能会存在一些性能问题,特别是使用!!
,以及由于遍历而仍然运行O(n²)操作的事实,正如@leftaroundabout所提到的那样,O(n·(n–1)) ≡ O(n²),但每次迭代应该调用gravitate
n*(n-1)/2
次。
O
(发音为“大O”)符号实际上并不是某个输入所需计算次数。它表示随着n趋近于无穷大时计算次数的增长率。另一种思考方式是“当n变得很大时,哪个因素最终占主导地位”。虽然在评论中深入讨论可能有些困难! - Michal Charemza