[Haskell-cafe] fast Eucl. dist. - Haskell vs C

[Daniel Schüssler]

Hello!

On Monday 18 May 2009 14:37:51 Kenneth Hoste wrote:
> I'm mostly interested in the range 10D to 100D

is the dimension known at compile-time? Then you could consider Template 
Haskell. I wrote up some code for generating the vector types and vector 
subtraction/inner product below, HTH. One problem is that I'm using a 
typeclass and apparently you can't make {-# SPECIALISE #-} pragmas with TH, 
so let's hope it is automatically specialised by GHC.

Greetings,
Daniel

TH.hs
------------------


{-# LANGUAGE TemplateHaskell #-}
{-# OPTIONS -fglasgow-exts #-}

module TH where
    
import Language.Haskell.TH
import Control.Monad

-- Non-TH stuff
class InnerProductSpace v r | v -> r where
    innerProduct :: v -> v -> r

class AbGroup v where
    minus :: v -> v -> v
            
euclidean x y = case minus x y of
                  z -> sqrt $! innerProduct z z

-- TH
noContext :: Q Cxt
noContext = return []
strict :: Q Type -> StrictTypeQ
strict = liftM ((,) IsStrict)
                      
makeVectors :: Int -- ^ Dimension
            -> Q Type -- ^ Component type, assumed to be a 'Num'
            -> String -- ^ Name for the generated type
            -> Q [Dec]
makeVectors n ctyp name0 = do
  -- let's assume ctyp = Double, name = Vector for the comments
  
  -- generate names for the variables we will need
  xs <- replicateM n (newName "x")
  ys <- replicateM n (newName "y")
       
  let
      name = mkName name0
      
      -- shorthands for arithmetic expressions; the first takes expressions,
      -- the others take variable names 
      sumE  e1 e2     = infixE (Just       e1)  [|(+)|] (Just       e2)
      varDiffE e1 e2  = infixE (Just (varE e1)) [|(-)|] (Just (varE e2))
      varProdE e1 e2  = infixE (Just (varE e1)) [|(*)|] (Just (varE e2))

      
      conPat vars = conP name (fmap varP vars)
      
      -- > data Vector = Vector !Double ... !Double
      theDataD = 
          dataD noContext name [] -- no context, no params
                    [normalC name (replicate n (strict ctyp))] 
                    [''Eq,''Ord,''Show] -- 'deriving' clause
                  
      innerProdD = 
          -- > instance InnerProductSpace Vector Double where ...
          instanceD noContext ( conT ''InnerProductSpace 
                                `appT` conT name 
                                `appT` ctyp)  
                    -- > innerProduct = ...
                    [valD 
                     (varP 'innerProduct)
                     (normalB
                      -- \(Vector x1 x2 ... xn) (Vector y1 y2 ... yn) ->
                      (lamE [conPat xs, conPat ys]
                       -- x1*y1 + .... + xn*yn + 0
                       (foldl sumE [|0|] $
                              zipWith varProdE xs ys)
                       ))
                       
                     [] -- no 'where' clause
                    ]
          
      abGroupD = 
          instanceD noContext ( conT ''AbGroup
                                `appT` conT name)  
                    -- > minus = ...
                    [valD 
                     (varP 'minus)
                     (normalB
                      -- \(Vector x1 x2 ... xn) (Vector y1 y2 ... yn) ->
                      (lamE [conPat xs, conPat ys] 
                       -- Vector (x1-y1) ... (xn-yn)
                       (foldl appE (conE name) $
                              zipWith varDiffE xs ys)
                       ))
                       
                     [] -- no 'where' clause
                    ]
  

  sequence [theDataD,innerProdD,abGroupD]



Main.hs
------------------
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE MultiParamTypeClasses #-}

module Main where

import TH

$(makeVectors 3 [t|Double|] "Vec3")

main = print $ euclidean (Vec3 1 1 1) (Vec3 0 0 0)