Nested triangular loop in Haskell?

I have the following in Java which basically does a nested triangular loop:

    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);
    }

Class B:

class B {
    double x, v, m;
}

In each iteration, the value at index i and j of the array is updated at the same time thus avoiding to do a complete nested loop. The following gives the same result but it does a complete nested loop (excuse me for the terms i'm using, they may not be correct but i hope it does make sense).

    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;
        }
    }

NOTE: the only change from the previous code is int j = 0; NOT int j = (i+1); and remove开发者_开发知识库d bs[j].v += d * bs[i].m * mag;

I want to do same in Haskell but having difficulty to think about it properly. I have the following code. The array in the Haskell version is represented as a list (xs) which i've initialised to 0.

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

I am going to have very big value for n e.g. 1000, 5000, etc. A complete nested loop when n = 1000 gives length [(i,j)|i<-[0..1000],j<-[0..1000]] = 1002001 but a triangular version gives length [(i,j)|i<-[0..1000],j<-[(i+1)..1000]] = 500500. Doing 2 maps in Haskell is easy to get it to do the complete loops but I want the triangular version. I guess this implies keeping the changes to i and j in a list and then update the original list at the end? Any idea would be much appreciated. Thanks

---

Here's a straightforward translation using unboxed mutable vectors from the vector package. Code is somewhat ugly, but should be very fast:

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 ()

Update: Regarding the "very fast" claim: I benchmarked my solution against the pure one provided by Federico and got the following results (for 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

So the imperative solution is approx. 50 times faster than the functional one (the difference is even more dramatic for smaller n, when everything fits in cache). I tried to make Federico's solution work with unboxed vectors, but apparently it relies on laziness in a crucial way, which makes the unboxed version loop forever. The "pure vector" version uses boxed vectors.

---

I'm not sure this solves your problem because I didn't grasp it completely yet, but the triangular loop itself is very easy to do in Haskell:

triangularLoop :: (a -> a -> b) -> [a] -> [b]
triangularLoop f xs = do
    (x1 : t) <- tails xs
    x2 <- t
    return $ f x1 x2

Or, written without the monadic syntax,

triangularLoop f = concat . map singlePass . tails
    where
        singlePass [] = []
        singlePass (h:t) = map (f h) t

---

A typical, idiomatic way of writing nested loops in Haskell is using list comprehensions.

Here is how I would translate your code:

import Data.Array

import Data.List (tails)

data Body = Body {x::Double,v::Double,m::Double}
            deriving Show

n::Int
n = 9

dt::Double
dt = 0.001

bs_0 :: Array Int Body
bs_0 = array (0,n) [(i,Body {x = i'*0.5,v = i'*2.5,m = i'*5.5}) | 
                    i <- [0..n], let i' = fromIntegral i]

bs :: Array Int Body
bs = accum (\b dv -> b {v = v b + dv}) bs_0 dvs
     where 
       dvs :: [(Int,Double)]
       dvs = concat [[(i,dv_i),(j,dv_j)] | (i:is) <- tails [0..n], 
                                            j <- is,
                                            let d = x(bs!i) - x(bs!j)
                                                sqr =  d * d + 0.01
                                                dist = sqrt sqr
                                                mag = dt / (sqr * dist)
                                                dv_i = -d * m(bs!j) * mag
                                                dv_j =  d * m(bs!i) * mag]

main :: IO()
main = mapM_ print (assocs bs)