Haskell中力的共享计算
我要在Haskell实现一个N体模拟。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的力的相关加速度，然后在ps加速中更新每个Particle。
我想过要做一些像foldr f ps combs这样的事情。启动累加器将是当前的ps，f是一些函数，它接受每个comb并在ps中更新相关的Particle，并返回该累加器。对于这样一个简单的过程来说，这看起来确实是内存密集型的，并且相当复杂。
有什么想法吗？
回答 1
Stack Overflow用户
发布于 2017-06-19 21:50:57
从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 = -1 * gravitate b a这一事实计算accsMatrix时减少对gravitate a b的调用数。
如果我们打印出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..。或者用单子..。
accsMatrix = map accs indexedPs >>= (:[]) . flip map indexedPs..。虽然(对我来说)很难看到这里面发生了什么。
这个列表方法可能存在一些可怕的性能问题，特别是使用!!，以及由于遍历而仍在运行O (n 2)操作，以及O (n·(n-1)≡O(N 2))为@leftaroundabout所述的事实，但每次迭代都应该调用gravitate n * (n-1) / 2 times。
https://stackoverflow.com/questions/44622757
复制相似问题
腾讯云开发者
