haskell - Haskell中的嵌套三角形循环？

Question

我在 Java 中有以下内容，它基本上执行嵌套的三角循环:

    int n = 10;
    B bs[] = new B[n];
    // some initial values, bla bla
    double dt = 0.001;
    for (int i = 0; i < n; i++) {
        bs[i] = new B();
        bs[i].x = i * 0.5;
        bs[i].v = i * 2.5;
        bs[i].m = i * 5.5;
    }
    for (int i = 0; i < n; i++) {
        for (int j = **(i+1)**; j < n; j++) {
            double d = bs[i].x - bs[j].x;

            double sqr = d * d + 0.01;
            double dist = Math.sqrt(sqr);
            double mag = dt / (sqr * dist);

            bs[i].v -= d * bs[j].m * mag;
            **bs[j].v += d * bs[i].m * mag;**
        }
    }   

    // printing out the value v
    for (int i = 0; i < n; i++) {
        System.out.println(bs[i].v);
    }

B 类:

class B {
    double x, v, m;
}

在每次迭代中，数组索引 i 和 j 处的值同时更新，从而避免进行完整的嵌套循环。以下给出了相同的结果，但它执行了一个完整的嵌套循环(请原谅我使用的术语，它们可能不正确，但我希望它确实有意义)。

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            double d = bs[i].x - bs[j].x;

            double sqr = d * d + 0.01;
            double dist = Math.sqrt(sqr);
            double mag = dt / (sqr * dist);

            bs[i].v -= d * bs[j].m * mag;
        }
    }

注意: 与之前代码的唯一变化是 int j = 0;不是 int j = (i+1);并删除了 bs[j].v += d * bs[i].m * mag;

我想在 Haskell 中做同样的事情，但很难正确地考虑它。我有以下代码。 Haskell 版本中的数组表示为我已初始化为 0 的列表 (xs)。

n = 20
xs = replicate n 0

update = foldl' (update') xs [0..(n-1)]
    where
        update' i = update'' i (i+1) []
        update'' i j acc
            | j == n = acc
            | otherwise = new_acc
                where
                    new_acc = result:acc
                    result = ...do something

我将拥有非常大的值(value)，例如1000、5000 等 n = 1000 时的完整嵌套循环给出 length [(i,j)|i<-[0..1000],j<-[0..1000]] = 1002001但三角形版本给出 length [(i,j)|i<-[0..1000],j<-[(i+1)..1000]] = 500500 .在 Haskell 中做 2 个 map 很容易让它完成完整的循环，但我想要三角形版本。我想这意味着将对 i 和 j 的更改保存在列表中，然后在最后更新原始列表？任何想法将不胜感激。谢谢

Answer 1

这是使用来自 vector 的未装箱可变向量的直接翻译包裹。代码有点难看，但应该很快:

module Main
    where

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as M

numElts :: Int
numElts = 10

dt :: Double
dt = 0.001

loop :: Int -> M.IOVector Double -> M.IOVector Double 
        -> M.IOVector Double -> IO ()
loop n x v m = go 0
  where
    doWork i j = do xI <- M.read x i
                    xJ <- M.read x j
                    vI <- M.read v i
                    vJ <- M.read v j
                    mI <- M.read m i
                    mJ <- M.read m j

                    let d = xI - xJ
                    let sqr = d * d + 0.01
                    let dist = sqrt sqr
                    let mag = dt / (sqr * dist)

                    M.write v i (vI - d * mJ * mag)
                    M.write v j (vJ + d * mI * mag)

    go i | i < n     = do go' (i+1)
                          go  (i+1)
         | otherwise = return ()
      where
        go' j | j < n     = do doWork i j
                               go' (j + 1)
              | otherwise = return ()

main :: IO ()
main = do x <- generateVector 0.5
          v <- generateVector 2.5
          m <- generateVector 5.5
          loop numElts x v m
          v' <- U.unsafeFreeze v
          U.forM_ v' print

    where
      generateVector :: Double -> IO (M.IOVector Double)
      generateVector d = do v <- M.new numElts
                            generateVector' numElts d v
                            return v

      generateVector' :: Int -> Double -> M.IOVector Double -> IO ()
      generateVector' n d v = go 0
        where
          go i | i < n = do M.unsafeWrite v i (fromIntegral i * d)
                            go (i+1)
               | otherwise = return ()

更新:关于“非常快”的说法:我benchmarked my solution与 Federico 提供的 pure 对比并得到以下结果(对于 n = 1000):

benchmarking pureSolution
collecting 100 samples, 1 iterations each, in estimated 334.5483 s
mean: 2.949640 s, lb 2.867693 s, ub 3.005429 s, ci 0.950
std dev: 421.1978 ms, lb 343.8233 ms, ub 539.4906 ms, ci 0.950
found 4 outliers among 100 samples (4.0%)
  3 (3.0%) high severe
variance introduced by outliers: 5.997%
variance is slightly inflated by outliers

benchmarking pureVectorSolution
collecting 100 samples, 1 iterations each, in estimated 280.4593 s
mean: 2.747359 s, lb 2.709507 s, ub 2.803392 s, ci 0.950
std dev: 237.7489 ms, lb 179.3110 ms, ub 311.8813 ms, ci 0.950
found 13 outliers among 100 samples (13.0%)
  7 (7.0%) high mild
  6 (6.0%) high severe
variance introduced by outliers: 2.998%
variance is slightly inflated by outliers

benchmarking imperativeSolution
collecting 100 samples, 1 iterations each, in estimated 5.905104 s
mean: 58.59154 ms, lb 56.79405 ms, ub 60.60033 ms, ci 0.950
std dev: 11.70101 ms, lb 9.120100 ms, ub NaN s, ci 0.950

所以当务之急的解决方案是大约。比功能性快 50 倍(当一切都适合缓存时，对于较小的 n，差异甚至更加显着)。我试图让 Federico 的解决方案与未装箱的向量一起工作，但显然它以一种关键的方式依赖于惰性，这使得未装箱的版本永远循环。 “纯矢量”版本使用盒装矢量。