f# - 在 F# 中展开嵌套循环

Question

我一直在努力使用以下代码。它是 Forward-Euler algorithm 的 F# 实现。用于模拟在引力场中运动的恒星。

let force (b1:Body) (b2:Body) = 
    let r = (b2.Position - b1.Position)
    let rm = (float32)r.MagnitudeSquared + softeningLengthSquared
    if (b1 = b2) then
        VectorFloat.Zero
    else
        r * (b1.Mass * b2.Mass) / (Math.Sqrt((float)rm) * (float)rm)    

member this.Integrate(dT, (bodies:Body[])) = 

    for i = 0 to bodies.Length - 1 do
        for j = (i + 1) to bodies.Length - 1 do
            let f = force bodies.[i] bodies.[j]
            bodies.[i].Acceleration <- bodies.[i].Acceleration + (f / bodies.[i].Mass)
            bodies.[j].Acceleration <- bodies.[j].Acceleration - (f / bodies.[j].Mass)
        bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
        bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT   

虽然这可行，但它并不完全是“功能性的”。它的性能也很糟糕，比等效的 c# 代码慢 2.5 倍。 body 是 Body 类型的结构数组。

我正在努力解决的问题是 force() 是一个昂贵的函数，所以通常你为每对计算一次并依赖 Fij = -Fji 的事实。但这确实会打乱任何展开的循环等。

建议已收到!不，这不是家庭作业...

谢谢，

阿德

更新:澄清 Body 和 VectorFloat 被定义为 C# 结构。这是因为程序在 F#/C# 和 C++/CLI 之间进行互操作。最终，我将在 BitBucket 上获取代码，但这是一项正在进行的工作，我有一些问题需要解决，然后才能将其发布。

[StructLayout(LayoutKind.Sequential)]
public struct Body
{
    public VectorFloat Position;
    public float Size;
    public uint Color;

    public VectorFloat Velocity;
    public VectorFloat Acceleration;
              '''
}   

[StructLayout(LayoutKind.Sequential)]
public partial struct VectorFloat
{
    public System.Single X { get; set; }
    public System.Single Y { get; set; }
    public System.Single Z { get; set; }
}

向量定义了标准向量类所期望的运算符类型。对于这种情况，您可能可以使用 .NET 框架中的 Vector3D 类(我实际上正在研究切换到它)。

更新 2:根据下面的前两个回复改进了代码:

    for i = 0 to bodies.Length - 1 do
    for j = (i + 1) to bodies.Length - 1 do
        let r = ( bodies.[j].Position -  bodies.[i].Position)
        let rm = (float32)r.MagnitudeSquared + softeningLengthSquared
        let f = r / (Math.Sqrt((float)rm) * (float)rm)    
        bodies.[i].Acceleration <- bodies.[i].Acceleration + (f * bodies.[j].Mass)
        bodies.[j].Acceleration <- bodies.[j].Acceleration - (f * bodies.[i].Mass)
    bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
    bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT  

force 函数中覆盖 b1 == b2 情况的分支是最严重的违规者。如果 softeningLength 始终非零，即使它非常小(Epsilon)，您也不需要它。此优化在 C# 代码中，但不在 F# 版本中(doh!)。

Math.Pow(x, -1.5) 似乎比 1/(Math.Sqrt(x) * x) 慢很多。本质上，这个算法有点奇怪，因为它的性能取决于这一步的成本。

将力计算内联并消除一些分歧也带来了一些改进，但性能确实被分支杀死并且由 Sqrt 的成本主导。

WRT 在结构上使用类:在某些情况下(此代码的 CUDA 和 native C++ 实现以及 DX9 渲染器)，我需要将主体数组放入非托管代码或 GPU 上。在这些情况下，能够 memcpy 一个连续的内存块似乎是要走的路。不是我从一组 Body 中得到的东西。

Answer 1

我不确定以功能样式重写此代码是否明智。我见过一些尝试以函数方式编写对交互计算的尝试，其中每一个都比两个嵌套循环更难遵循。

在查看结构与类之前(我相信其他人对此有聪明的说法)，也许您可​​以尝试优化计算本身？

您正在计算两个加速度增量，我们称它们为 dAi 和 dAj:

dAi = r*m1*m2/(rm*sqrt(rm))/m1

dAj = r*m1*m2/(rm*sqrt(rm))/m2

[注:m1 = body.[i].mass，m2=body.[j].mass]]

质量除法取消如下:

dAi = rm2/(rmsqrt(rm))

dAj = rm1/(rmsqrt(rm))

现在您只需计算每对 (i,j) 的 r/(rmsqrt(rm))。
这可以进一步优化，因为 1/(rmsqrt(rm)) = 1/(rm^1.5) = rm^-1.5，所以如果让 r' = r * (rm ** -1.5)，那么编辑:否它不能，那是过早的优化在那里讨论(见评论)。计算 r' = 1.0/(r * sqrt r) 是最快的。

dAi = m2 * r'

dAj = m1 * r'

您的代码将变成类似

member this.Integrate(dT, (bodies:Body[])) = 
    for i = 0 to bodies.Length - 1 do
        for j = (i + 1) to bodies.Length - 1 do
            let r = (b2.Position - b1.Position)
            let rm = (float32)r.MagnitudeSquared + softeningLengthSquared            
            let r' = r * (rm ** -1.5)
            bodies.[i].Acceleration <- bodies.[i].Acceleration + r' * bodies.[j].Mass
            bodies.[j].Acceleration <- bodies.[j].Acceleration - r' * bodies.[i].Mass
        bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
        bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT

看，妈，别再 split 了!

警告:未经测试的代码。尝试自担风险。