Matrix Multiplication

use Array: all;
 
double[.,.] matmul (double[.,.] a, double[.,.] b)
{
  // transpose
  bt = { [i,j] -> b[j,i] };
  // multiplication
  return { iv -> sum (a[iv[0]] * bt[iv[1]])
           | iv < [shape(a)[0],shape(b)[1]] };
}
 
int main()
{
  // create matrices
  A = reshape ([6,8], tod (iota (48)));
  B = reshape ([8,2], tod (iota (16)));
 
  C = matmul (A, B);
 
  StdIO::print (C);
 
  return 0;
}

Parallel Scan

int [+], int scan (int [+] a)
{
  // TODO TC multiple return not supported yet
  a, m = { iv -> scan (a[iv ])
           | iv < drop ([ -1] , shape (a)) };
  m, s = scan (m);
  a = a ^+ m;
  return (a, s);
}

\begin{eqnarray} \overset{\tiny k+1}{p_i} &=& \overset{k}{p_i} + \overset{k}{v_i} dt \\ \overset{k+1}{v_i} &=& \overset{k}{v_i} + \overset{k}{a_i} dt \\ \overset{k+1}{a_i} &=& \sum\limits_{j \neq i}^{n} \dfrac{m_j (\overset{k+1}{p_j} - \overset{k+1}{p_i})} {\left|\overset{k+1}{p_j} - \overset{k+1}{p_i} \right|^3} \end{eqnarray}

p = p + v * dt;
v = v + a * dt;
a = {[i] 
      -> vsum ({[j]
                 -> (i == j  ? [0.0, 0.0, 0.0] 
                             : m[j] * (p[j] - p[i])
                    / (l2norm (p[j] - p[i]) ^ 3) 
              })
    };