Inspired by a Stack Overflow question, I decided to practice my Haskell by taking a crack at the Google Code Jam's Minimum Scalar Product problem:
Given two vectors \$\mathbf{v_1}=(x_1,x_2,\ldots,x_n)\$ and \$\mathbf{v_2}=(y_1,y_2,\ldots,y_n)\$. If you could permute the coordinates of each vector, what is the minimum scalar product \$\mathbf{v_1} \cdot \mathbf{v_2}\$?
Constraints:
\$100 \le n \le 800\$
\$-100000 \le x_i, y_i \le 100000\$
I'm not claiming any algorithmic awesomeness (this is just a reference implementation for checking correctness later).
This is my first time using Vectors and the ST monad, so what I really want is a sanity check that I'm using both correctly, and that I'm using the correct tools for the job.
module MinimumScalarProduct where
import Control.Monad (replicateM, forM_, (>=>))
import Control.Monad.ST (runST, ST)
import Data.Vector (thaw, fromList, MVector, Vector, zipWith, foldl', length, (!))
import Data.Vector.Generic.Mutable (read, swap)
import Prelude hiding (read, zipWith, length)
-- sequnce of transpoitions to yield all permutations of n elts
-- http://en.wikipedia.org/wiki/Steinhaus-Johnson-Trotter_algorithm
transpositions :: Int -> [(Int, Int)]
transpositions n = runST $ do
-- p maps index to element
p <- thaw $ fromList [0 .. n-1]
-- q maps element to index
q <- thaw $ fromList [0 .. n-1]
-- let the prefixes define themselves recursively
foldr ((>=>) . extend p q) return [2..n] []
extend :: MVector s Int -> MVector s Int -> Int -> [(Int,Int)] -> ST s [(Int,Int)]
extend p q n ts = fmap ((ts++) . concat) . replicateM (n-1) $ do
-- replace the element in the (n-1)'th position with its predecessor
a <- read p (n-1)
i <- read q (a-1)
swap p (n-1) i
swap q a (a-1)
-- replay the earlier transpositions
forM_ ts $ \(m,j) -> do
b <- read p m
c <- read p j
swap p m j
swap q b c
return $ (n-1,i) : ts
-- reference implementation, takes O(n!)
msp :: Vector Int -> Vector Int -> Int
msp u v | length u /= length v = 0
msp u v = runST $ do
let x = foldl' (+) 0 $ zipWith (*) u v
let n = length u
u' <- thaw u
-- check each permutation of u'
let steps = map (adjust u' v) $ transpositions n
fmap minimum . sequence $ scanl (>>=) (return x) steps
-- adjust the current scalar product for the transposition of u
adjust :: MVector s Int -> Vector Int -> (Int, Int) -> Int -> ST s Int
adjust u v (i,j) x = do
a <- read u i
b <- read u j
let c = v ! i
let d = v ! j
swap u i j
return $ x - (a*c + b*d) + (a*d + b*c)
