diff --git a/app/src/BoundingBox.hs b/app/src/BoundingBox.hs
--- a/app/src/BoundingBox.hs
+++ b/app/src/BoundingBox.hs
@@ -4,6 +4,7 @@
 
 module BoundingBox where
 
+import PolymorphicNum
 import Vector
 import GHC.Types
 import GHC.Prim
diff --git a/app/src/Camera.hs b/app/src/Camera.hs
--- a/app/src/Camera.hs
+++ b/app/src/Camera.hs
@@ -2,21 +2,22 @@
 
 module Camera where
 
+import PolymorphicNum
 import Vector
 import Matrix
 
-data Camera = Camera { worldToCamera :: Matrix, fieldOfView :: !Double, position :: Vector } deriving (Show)
+data Camera = Camera { worldToCamera :: Matrix, fieldOfView :: !Double, position :: Vector, farClip :: !Double } deriving (Show)
 
-lookAt :: Position -> Position -> Direction -> Double -> Camera
-lookAt pos target up fov = 
-    Camera matrix fov pos
+lookAt :: Position -> Position -> Direction -> Double -> Double -> Camera
+lookAt pos target up fov dist = 
+    Camera matrix fov pos dist
     where
-      forward = normalise $ target - pos
+      forward = normalise $ target <-> pos
       right = normalise $ up `cross` forward
       up' = right `cross` forward
       matrix = buildMatrix right up' forward pos
 
-withVectors :: Position -> Direction -> Direction -> Direction -> Double -> Camera
-withVectors pos basisX basisY basisZ fov = Camera matrix fov pos
+withVectors :: Position -> Direction -> Direction -> Direction -> Double -> Double -> Camera
+withVectors pos basisX basisY basisZ fov dist = Camera matrix fov pos dist
     where
       matrix = buildMatrix basisX basisY basisZ (Vector.negate pos)
diff --git a/app/src/Colour.hs b/app/src/Colour.hs
--- a/app/src/Colour.hs
+++ b/app/src/Colour.hs
@@ -1,11 +1,12 @@
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 module Colour where
 import Vector hiding (min, max)
 import Misc
 import Data.Word
 import Control.DeepSeq
+import PolymorphicNum
 
 -- Normalised RGBA colour
 data Colour = Colour { red :: {-# UNPACK #-} !Double, 
@@ -13,41 +14,26 @@
                        blue :: {-# UNPACK #-} !Double, 
                        alpha :: {-# UNPACK #-} !Double } deriving (Show, Read, Ord, Eq)
 
-instance Num Colour where
-    {-# SPECIALIZE INLINE (+) :: Colour -> Colour -> Colour #-}
-    (Colour !r1 !g1 !b1 !a1) + (Colour !r2 !g2 !b2 !a2) = Colour (r1 + r2) (g1 + g2) (b1 + b2) (a1 + a2)
-    {-# SPECIALIZE INLINE (-) :: Colour -> Colour -> Colour #-}
-    (Colour !r1 !g1 !b1 !a1) - (Colour !r2 !g2 !b2 !a2) = Colour (r1 - r2) (g1 - g2) (b1 - b2) (a1 - a2)
-    {-# SPECIALIZE INLINE (*) :: Colour -> Colour -> Colour #-}
-    (Colour !r1 !g1 !b1 !a1) * (Colour !r2 !g2 !b2 !a2) = Colour (r1 * r2) (g1 * g2) (b1 * b2) (a1 * a2)
-    abs (Colour r g b a) = Colour (abs r) (abs g) (abs b) (abs a)
-    signum (Colour r g b a) = Colour (signum r) (signum g) (signum b) (signum a)
-    fromInteger x = Colour (fromInteger x) (fromInteger x) (fromInteger x) (fromInteger x)
-
-instance Fractional Colour where
-    {-# SPECIALIZE INLINE (/) :: Colour -> Colour -> Colour #-}
-    (Colour !r1 !g1 !b1 !a1) / (Colour !r2 !g2 !b2 !a2) = Colour (r1 / r2) (g1 / g2) (b1 / b2) (a1 / a2)
-    fromRational x = Colour (fromRational x) (fromRational x) (fromRational x) (fromRational x)
-
 instance NFData Colour where
     rnf (Colour r g b a) = rnf r `seq` rnf g `seq` rnf b `seq` rnf a
 
-infixl 7 <*>
-infixl 7 </>
-infixl 6 <+>
-infixl 6 <->
-
-(<*>) :: Colour -> Double -> Colour
-(Colour !r !g !b !a) <*> k = Colour (r * k) (g * k) (b * k) (a * k)
-
-(</>) :: Colour -> Double -> Colour
-(Colour !r !g !b !a) </> k = Colour (r / k) (g / k) (b / k) (a / k)
+instance PolymorphicNum Colour Colour Colour where
+    (Colour !r !g !b !a) <*> (Colour !r' !g' !b' !a') = Colour (r * r') (g * g') (b * b') (a * a')
+    (Colour !r !g !b !a) </> (Colour !r' !g' !b' !a') = Colour (r / r') (g / g') (b / b') (a / a')
+    (Colour !r !g !b !a) <-> (Colour !r' !g' !b' !a') = Colour (r - r') (g - g') (b - b') (a - a')
+    (Colour !r !g !b !a) <+> (Colour !r' !g' !b' !a') = Colour (r + r') (g + g') (b + b') (a + a')
 
-(<+>) :: Colour -> Double -> Colour
-(Colour !r !g !b !a) <+> k = Colour (r + k) (g + k) (b + k) (a + k)
+instance PolymorphicNum Colour Double Colour where
+    (Colour !r !g !b !a) <*> k = Colour (r * k) (g * k) (b * k) (a * k)
+    (Colour !r !g !b !a) </> k = Colour (r / k) (g / k) (b / k) (a / k)
+    (Colour !r !g !b !a) <-> k = Colour (r - k) (g - k) (b - k) (a - k)
+    (Colour !r !g !b !a) <+> k = Colour (r + k) (g + k) (b + k) (a + k)
 
-(<->) :: Colour -> Double -> Colour
-(Colour !r !g !b !a) <-> k = Colour (r - k) (g - k) (b - k) (a - k)
+instance PolymorphicNum Double Colour Colour where
+    k <*> (Colour !r !g !b !a) = Colour (k * r) (k * g) (k * b) (k * a)
+    k </> (Colour !r !g !b !a) = Colour (k / r) (k / g) (k / b) (k / a)
+    k <-> (Colour !r !g !b !a) = Colour (k - r) (k - g) (k - b) (k - a)
+    k <+> (Colour !r !g !b !a) = Colour (k + r) (k + g) (k + b) (k + a)
 
 clamp :: Colour -> Colour
 clamp (Colour !r !g !b !a) = Colour (max 0 (min r 1)) (max 0 (min g 1)) (max 0 (min b 1)) (max 0 (min a 1))
@@ -70,6 +56,9 @@
 
 colBlack :: Colour
 colBlack = Colour 0 0 0 1
+
+colZero :: Colour
+colZero = Colour 0 0 0 0
 
 colGrey :: Colour
 colGrey = Colour 0.5 0.5 0.5 1
diff --git a/app/src/CornellBox.hs b/app/src/CornellBox.hs
--- a/app/src/CornellBox.hs
+++ b/app/src/CornellBox.hs
@@ -1,5 +1,4 @@
 -- Cornell box reference data
-{-# LANGUAGE MagicHash #-}
 
 module CornellBox(cornellBox, cornellBoxCamera, cornellBoxLights) where
 
@@ -27,6 +26,7 @@
 backWallObject :: Object
 tallBlockObject :: Object
 shortBlockObject :: Object
+--lightObject :: Object
 
 leftWallVertices :: [Vector]
 rightWallVertices :: [Vector]
@@ -41,18 +41,27 @@
 whiteMaterial :: Material
 redMaterial :: Material
 greenMaterial :: Material
+--lightMaterial :: Material
 
 cornellBoxLights = [ 
     QuadLight (CommonLightData (Colour 500 500 500 0) True) (Vector 213.0 548.0 227.0 1.0) 600 (Vector 130.0 0.0 0.0 0.0) (Vector 0.0 0.0 105.0 0.0)
     ]
 
 cameraPosition = Vector 278.0 273.0 (-800.0) 1.0
-cornellBoxCamera = withVectors cameraPosition xaxis yaxis zaxis 45.0
+cornellBoxCamera = withVectors cameraPosition xaxis yaxis zaxis 45.0 10000
 
-whiteMaterial = Material (Colour 0.5 0.5 0.5 1) (Colour 0.5 0.5 0.5 1) colBlack 0 0 0 iorAir NullShader
-redMaterial   = Material (Colour 0.5 0.0 0.0 1) (Colour 0.5 0.0 0.0 1) colBlack 0 0 0 iorAir NullShader
-greenMaterial = Material (Colour 0.0 0.5 0.0 1) (Colour 0.0 0.5 0.0 1) colBlack 0 0 0 iorAir NullShader
+whiteMaterial = Material (Colour 0.5 0.5 0.5 1) (Colour 0.5 0.5 0.5 1) colBlack colBlack 0 0 0 iorAir NullShader
+redMaterial   = Material (Colour 0.5 0.0 0.0 1) (Colour 0.5 0.0 0.0 1) colBlack colBlack 0 0 0 iorAir NullShader
+greenMaterial = Material (Colour 0.0 0.5 0.0 1) (Colour 0.0 0.5 0.0 1) colBlack colBlack 0 0 0 iorAir NullShader
+--lightMaterial = Material colBlack colBlack colBlack (Colour 1000 1000 1000 1) 0 0 0 iorAir NullShader
 
+--lightVertices = [
+--    Vector 343.0 548.0 227.0 1.0,
+--    Vector 343.0 548.0 342.2 1.0,
+--    Vector 213.0 548.0 342.0 1.0,
+--    Vector 213.0 548.0 227.2 1.0
+--    ]
+
 floorVertices = [
     Vector 556.0 0.0   0.0 1.0,
     Vector   0.0 0.0   0.0 1.0,
@@ -157,6 +166,7 @@
 backWallObject = Object (TriangleMesh (quadsToTriangles backWallVertices)) whiteMaterial identity
 shortBlockObject = Object (TriangleMesh (quadsToTriangles shortBlockVertices)) whiteMaterial identity
 tallBlockObject = Object (TriangleMesh (quadsToTriangles tallBlockVertices)) whiteMaterial identity
+--lightObject = Object (TriangleMesh (quadsToTriangles lightVertices)) lightMaterial identity
 
 cornellBox :: [Object]
-cornellBox = [ceilingObject, floorObject, leftWallObject, rightWallObject, backWallObject, frontWallObject, tallBlockObject, shortBlockObject]
+cornellBox = [ceilingObject, floorObject, leftWallObject, rightWallObject, backWallObject, frontWallObject, tallBlockObject, shortBlockObject{-, lightObject-}]
diff --git a/app/src/Distribution.hs b/app/src/Distribution.hs
--- a/app/src/Distribution.hs
+++ b/app/src/Distribution.hs
@@ -1,9 +1,14 @@
 -- Module for generating sample patterns for distributed ray tracing
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 
-module Distribution (generatePointsOnSphere, generatePointsOnQuad, generatePointsOnHemisphere) where
+module Distribution (generatePointsOnSphere, 
+                     generatePointsOnQuad, 
+                     generatePointsOnHemisphere,
+                     generatePointOnHemisphere,
+                     generateRandomUVs,
+                     randomUV) where
 
+import PolymorphicNum
 import Vector
 import System.Random.Mersenne.Pure64
 import Control.Monad.State
@@ -22,34 +27,49 @@
 generateRandomUVs :: Int -> GeneratorState [(Double, Double)]
 generateRandomUVs n = replicateM n randomUV
 
+uvToSphere :: Double -> (Double, Double) -> Position
+uvToSphere r (!u, !v) = Vector (r * x) (r * y) (r * z) 1
+    where
+      !z = 2 * u - 1
+      !t = 2 * pi * v
+      !w = sqrt (1 - z * z)
+      !x = w * cos t
+      !y = w * sin t
+
+uvToHemisphere :: Double -> (Double, Double) -> Position
+uvToHemisphere r (!u, !v) = Vector (r * x) (r * y) (r * z) 1
+    where
+      !z = v
+      !t = 2 * pi * u
+      !w = sqrt (1 - v * v)
+      !x = w * cos t
+      !y = w * sin t
+
 -- Generate a list of random points on a sphere
 generatePointsOnSphere :: Int -> Double -> Int -> [Position]
-generatePointsOnSphere numPoints r seed = map uvToPosition randomUVs
+generatePointsOnSphere numPoints r seed 
+    | numPoints <= 1 = [Vector 0 0 0 1]
+    | otherwise = map (uvToSphere r) randomUVs
     where
       randomUVs = evalState (generateRandomUVs numPoints) (pureMT (fromIntegral seed))
-      uvToPosition (!u, !v) = Vector (r * x) (r * y) (r * z) 1
-          where
-            !z = 2 * u - 1
-            !t = 2 * pi * v
-            !w = sqrt (1 - z * z)
-            !x = w * cos t
-            !y = w * sin t
 
 -- Generate a list of random points on a hemisphere (z > 0)
 generatePointsOnHemisphere :: Int -> Double -> Int -> [Position]
-generatePointsOnHemisphere numPoints r seed = map uvToPosition randomUVs
+generatePointsOnHemisphere numPoints r seed
+    | numPoints <= 1 = [Vector 0 0 0 1]
+    | otherwise = map (uvToHemisphere r) randomUVs
     where
       randomUVs = evalState (generateRandomUVs numPoints) (pureMT (fromIntegral seed))
-      uvToPosition (!u, !v) = Vector (r * x) (r * y) (r * z) 1
-          where
-            !y = 2 * u - 1
-            !t = pi * v
-            !w = sqrt (1 - y * y)
-            !x = w * cos t
-            !z = w * sin t
 
 generatePointsOnQuad :: Position -> Direction -> Direction -> Int -> Int -> [Position]
-generatePointsOnQuad pos deltaU deltaV numPoints seed = map uvToPosition randomUVs
+generatePointsOnQuad pos deltaU deltaV numPoints seed 
+    | numPoints <= 1 = [Vector 0 0 0 1]
+    | otherwise = map (\(u, v) -> pos <+> deltaU <*> u <+> deltaV <*> v) randomUVs
     where
       randomUVs = evalState (generateRandomUVs numPoints) (pureMT (fromIntegral seed))
-      uvToPosition (u, v) = pos + (deltaU <*> u) + (deltaV <*> v)
+
+-- Generate a single random point on a hemisphere
+generatePointOnHemisphere :: PureMT -> Double -> (Position, PureMT)
+generatePointOnHemisphere rndGen r  = (uvToHemisphere r uv, rndGen')
+    where
+      (uv, rndGen') = runState randomUV rndGen
diff --git a/app/src/IrradianceCache.hs b/app/src/IrradianceCache.hs
--- a/app/src/IrradianceCache.hs
+++ b/app/src/IrradianceCache.hs
@@ -1,9 +1,9 @@
 -- The irradiance cache
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 
 module IrradianceCache (IrradianceCache, query, initialiseCache) where
 
+import PolymorphicNum
 import Vector
 import Colour
 import BoundingBox
@@ -33,7 +33,11 @@
 -- The bigger the number, the better the estimate
 errorWeight :: (Position, Direction) -> (Position, CacheSample) -> Double
 {-# SPECIALIZE INLINE errorWeight :: (Position, Direction) -> (Position, CacheSample) -> Double #-}
-errorWeight (pos', dir') (pos, CacheSample (dir, _, r)) = 1 / ((pos `distance` pos') / r + sqrt (1 + (dir `sdot3` dir')))
+errorWeight (pos', dir') (pos, CacheSample (dir, _, r)) 
+    | dot <= 0 = 0
+    | otherwise = 1 / ((pos `distance` pos') / r + sqrt (1 + dot))
+    where
+      !dot = dir `dot3` dir'
 
 -- This slightly convoluted version is written to be tail recursive. I effectively have to maintain a software stack of the
 -- nodes remaining to be traversed
@@ -47,15 +51,15 @@
     where
       !weight = errorWeight posDir (samplePos, sample)
       (CacheSample (_, _, !sampleR)) = sample
-      minimumWeight = 0.4 -- This is approximately the lower bound of the weight at the radius of the sample
+      minimumWeight = 1.5 -- The bigger this weight, the less it will reuse samples and the higher the quality
 findSamples posDir (OctTreeDummy _ : xs) !acc = findSamples posDir xs acc
 findSamples _ [] !acc = acc
 
 -- Sum together a list of samples and error weights
 sumSamples :: [(Vector, CacheSample, Double)] -> Colour
-sumSamples !samples = colourSum Colour.</> weightSum
+sumSamples !samples = colourSum </> weightSum
     where
-      sumSamples' !(!colAcc, !weightAcc) ((_, CacheSample (_, !col, _), !weight):xs) = sumSamples' (colAcc + col Colour.<*> weight, weightAcc + weight) xs
+      sumSamples' !(!colAcc, !weightAcc) ((_, CacheSample (_, !col, _), !weight):xs) = sumSamples' (colAcc <+> col <*> weight, weightAcc + weight) xs
       sumSamples' !(!colAcc, !weightAcc) [] = (colAcc, weightAcc)
       !(!colourSum, !weightSum) = sumSamples' (colBlack, 0) samples
 
diff --git a/app/src/KDTree.hs b/app/src/KDTree.hs
--- a/app/src/KDTree.hs
+++ b/app/src/KDTree.hs
@@ -1,7 +1,8 @@
 -- This is a module for constructing bounding volume hierarchies using a kdtree
-{-# LANGUAGE MagicHash #-}
 
 module KDTree(generateSceneGraphUsingKDTree, makeSplittingPlane, degenerateSplitList, findSplittingPlane) where
+
+import PolymorphicNum
 import Vector
 import Primitive
 import Data.List
@@ -14,19 +15,19 @@
 onPositiveSide (planeNormal, planeDist) obj = planeDist + (planeNormal `dot3` objBoxCentre) > 0.01
     where
       Just (boxMin, boxMax) = primitiveBoundingBox (primitive obj) obj
-      objBoxCentre = (boxMin + boxMax) <*> 0.5
+      objBoxCentre = (boxMin <+> boxMax) <*> (0.5 :: Double)
 
 -- This stuff is generic
 
 -- Generate a plane to split the objects along
 makeSplittingPlane :: AABB -> Int -> (Vector, Double)
-makeSplittingPlane (boxMin, boxMax) buildCycle = case nthLargestAxis (boxMax - boxMin) buildCycle of
+makeSplittingPlane (boxMin, boxMax) buildCycle = case nthLargestAxis (boxMax <-> boxMin) buildCycle of
                                                    0 -> (xaxis, -(vecX midPoint))
                                                    1 -> (yaxis, -(vecY midPoint))
                                                    2 -> (zaxis, -(vecZ midPoint))
                                                    _ -> error "Undefined value"
     where
-      midPoint = (boxMin + boxMax) <*> 0.5
+      midPoint = (boxMin <+> boxMax) <*> (0.5 :: Double)
 
 -- Find a working splitting plane
 findSplittingPlane :: AABB -> Int -> [t] -> ((Vector, Double) -> t -> Bool) -> Maybe (Vector, Double) 
diff --git a/app/src/Light.hs b/app/src/Light.hs
--- a/app/src/Light.hs
+++ b/app/src/Light.hs
@@ -1,6 +1,5 @@
 -- Module for lights
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 
 module Light (applyLight, 
               surfaceEpsilon, 
@@ -15,6 +14,7 @@
               addToPhotonMap, 
               common) where
 
+import PolymorphicNum
 import Vector
 import Colour
 import Ray
@@ -48,60 +48,60 @@
 phongLighting (!shadePos, !tanSpace) (PointLight (CommonLightData !lightColour !inPhotonMap') !lightPos !lightRange) objMaterial sceneGraph !viewDirection 
     | (lightPos `distanceSq` shadePos) < (lightRange * lightRange) && dotProd > 0 = case findAnyIntersection sceneGraph (rayWithPoints intersectionPlusEpsilon lightPos) of
                                                                                       Just _ -> colBlack -- An object is closer to our point of consideration than the light, so occluded
-                                                                                      Nothing -> (lightColour * lightingSum) Colour.<*> attenuation
+                                                                                      Nothing -> (lightColour <*> lightingSum) <*> attenuation
                                                                                           where
-                                                                                            !lightingSum = diffuseLighting + specularLighting
-                                                                                            !attenuation = lightAttenuation lightPos shadePos lightRange
-                                                                                            !specularCorrection = (specularPower objMaterial + 2) / (2 * pi)
-                                                                                            !specularLighting = specular objMaterial Colour.<*> (specularCorrection * saturate (reflection `dot3` Vector.negate viewDirection) ** specularPower objMaterial)
-                                                                                            !reflection = reflect incoming normal
-                                                                                            !diffuseLighting = if inPhotonMap'
-                                                                                                               then colBlack
-                                                                                                               else shaderDiffuse * diffuse objMaterial Colour.<*> saturate dotProd
-                                                                                            !shaderDiffuse = evaluateDiffuse (shader objMaterial) shadePos tanSpace
+                                                                                            lightingSum = diffuseLighting <+> specularLighting
+                                                                                            attenuation = lightAttenuation lightPos shadePos lightRange
+                                                                                            specularCorrection = (specularPower objMaterial + 2) / (2 * pi)
+                                                                                            specularLighting = specular objMaterial <*> (specularCorrection * saturate (reflection `dot3` Vector.negate viewDirection) ** specularPower objMaterial)
+                                                                                            reflection = reflect incoming normal
+                                                                                            diffuseLighting = if inPhotonMap'
+                                                                                                              then colBlack
+                                                                                                              else shaderDiffuse <*> diffuse objMaterial <*> saturate dotProd
+                                                                                            shaderDiffuse = evaluateDiffuse (shader objMaterial) shadePos tanSpace
     | otherwise = colBlack
     where 
-      !intersectionPlusEpsilon = shadePos + (normal Vector.<*> surfaceEpsilon)
-      !incoming = normalise (lightPos - shadePos)
-      !dotProd = normal `dot3` incoming
-      !normal = thr tanSpace
+      intersectionPlusEpsilon = shadePos <+> normal <*> surfaceEpsilon
+      incoming = normalise (lightPos <-> shadePos)
+      dotProd = normal `dot3` incoming
+      normal = thr tanSpace
 phongLighting _ (AmbientLight (CommonLightData _ _)) _ _ _ = error "phongLighting: Do not know how to handle AmbientLight"
 phongLighting (!shadePos, !tanSpace) (QuadLight (CommonLightData !lightColour !inPhotonMap') !lightPos !lightRange !du !dv) objMaterial sceneGraph !viewDirection 
     | (lightCentre `distanceSq` shadePos) < (lightRange * lightRange) && dotProd > 0 = case findAnyIntersection sceneGraph (rayWithPoints intersectionPlusEpsilon lightCentre) of
                                                                                       Just _ -> colBlack -- An object is closer to our point of consideration than the light, so occluded
-                                                                                      Nothing -> (lightColour * lightingSum) Colour.<*> attenuation
+                                                                                      Nothing -> lightColour <*> lightingSum <*> attenuation
                                                                                           where
-                                                                                            !lightingSum = diffuseLighting + specularLighting
-                                                                                            !attenuation = lightAttenuation lightCentre shadePos lightRange
-                                                                                            !specularCorrection = (specularPower objMaterial + 2) / (2 * pi)
-                                                                                            !specularLighting = specular objMaterial Colour.<*> (specularCorrection * saturate (reflection `dot3` Vector.negate viewDirection) ** specularPower objMaterial)
-                                                                                            !reflection = reflect incoming normal
-                                                                                            !diffuseLighting = if inPhotonMap'
-                                                                                                               then colBlack
-                                                                                                               else shaderDiffuse * diffuse objMaterial Colour.<*> saturate dotProd
-                                                                                            !shaderDiffuse = evaluateDiffuse (shader objMaterial) shadePos tanSpace
+                                                                                            lightingSum = diffuseLighting <+> specularLighting
+                                                                                            attenuation = lightAttenuation lightCentre shadePos lightRange
+                                                                                            specularCorrection = (specularPower objMaterial + 2) / (2 * pi)
+                                                                                            specularLighting = specular objMaterial <*> (specularCorrection * saturate (reflection `dot3` Vector.negate viewDirection) ** specularPower objMaterial)
+                                                                                            reflection = reflect incoming normal
+                                                                                            diffuseLighting = if inPhotonMap'
+                                                                                                              then colBlack
+                                                                                                              else shaderDiffuse <*> diffuse objMaterial <*> saturate dotProd
+                                                                                            shaderDiffuse = evaluateDiffuse (shader objMaterial) shadePos tanSpace
     | otherwise = colBlack
     where 
-      !lightCentre = lightPos + du Vector.<*> 0.5 + dv Vector.<*> 0.5
-      !intersectionPlusEpsilon = shadePos + (normal Vector.<*> surfaceEpsilon)
-      !incoming = normalise (lightCentre - shadePos)
-      !dotProd = normal `dot3` incoming
-      !normal = thr tanSpace
+      lightCentre = lightPos <+> du <*> (0.5 :: Double) <+> dv <*> (0.5 :: Double)
+      intersectionPlusEpsilon = shadePos <+> normal <*> surfaceEpsilon
+      incoming = normalise (lightCentre <-> shadePos)
+      dotProd = normal `dot3` incoming
+      normal = thr tanSpace
 
 -- For a given surface point, work out the lighting, including occlusion
 applyLight :: SceneGraph -> SurfaceLocation -> Material -> Direction -> Light -> Colour
-applyLight sceneGraph !intersectionPointNormal !objMaterial !viewDirection light@(PointLight (CommonLightData _ _) _ _)
+applyLight sceneGraph intersectionPointNormal objMaterial viewDirection light@(PointLight (CommonLightData _ _) _ _)
     = phongLighting 
       intersectionPointNormal 
       light
       objMaterial 
       sceneGraph
       viewDirection
-applyLight _ (!intersectionPoint, !intersectionTanSpace) !objMaterial _ (AmbientLight (CommonLightData !ambientColour _)) = 
+applyLight _ (intersectionPoint, intersectionTanSpace) objMaterial _ (AmbientLight (CommonLightData ambientColour _)) = 
     let shaderAmbient = evaluateAmbient (shader objMaterial) intersectionPoint intersectionTanSpace
         materialAmbient = ambient objMaterial
-    in ambientColour * shaderAmbient * materialAmbient
-applyLight sceneGraph !intersectionPointNormal !objMaterial !viewDirection light@(QuadLight (CommonLightData _ _) _ _ _ _)
+    in ambientColour <*> shaderAmbient <*> materialAmbient
+applyLight sceneGraph intersectionPointNormal objMaterial viewDirection light@(QuadLight (CommonLightData _ _) _ _ _ _)
     = phongLighting 
       intersectionPointNormal 
       light
diff --git a/app/src/Main.hs b/app/src/Main.hs
--- a/app/src/Main.hs
+++ b/app/src/Main.hs
@@ -24,6 +24,7 @@
     | DirectPhotonMapVisualisation -- -v
     | DistributedRayTracing -- d
     | IrradianceCaching -- c
+    | PathTrace -- P
       deriving (Eq, Ord, Enum, Show, Bounded)
 
 options :: [OptDescr Option]
@@ -32,7 +33,8 @@
     Option ['p'] [] (NoArg PhotonMap) "Photon map",
     Option ['v'] [] (NoArg DirectPhotonMapVisualisation) "Direct photon map visualisation",
     Option ['d'] [] (NoArg DistributedRayTracing) "Distributed ray tracing",
-    Option ['c'] [] (NoArg IrradianceCaching) "Irradiance caching"
+    Option ['c'] [] (NoArg IrradianceCaching) "Irradiance caching",
+    Option ['P'] [] (NoArg PathTrace) "Path tracing"
     ]
 
 parsedOptions :: [String] -> [Option]
@@ -51,8 +53,10 @@
 renderImage :: Int -> RenderContext -> Maybe PhotonMap -> [Colour]
 renderImage mipLevel renderSettings photonMap = finalImage
     where
-      rawImageOutput = rayTraceImage renderSettings cornellBoxCamera (renderWidth mipLevel) (renderHeight mipLevel) photonMap
-      exposedImage = exposeImage imageAverageLogLuminance rawImageOutput 4
+      rawImageOutput = case renderMode renderSettings of
+                         PathTracer -> pathTraceImage renderSettings cornellBoxCamera (renderWidth mipLevel) (renderHeight mipLevel)
+                         _ -> rayTraceImage renderSettings cornellBoxCamera (renderWidth mipLevel) (renderHeight mipLevel) photonMap
+      exposedImage = exposeImage imageAverageLuminance rawImageOutput 4
       toneMappedImage = toneMapImage toneMapHejlBurgessDawson exposedImage
       finalImage = map (clamp . invGammaCorrect) toneMappedImage
 
@@ -116,6 +120,7 @@
          depthOfFieldFocalDistance' = 400
          renderMode'
              | PhotonMap `Prelude.elem` opts = PhotonMapper
+             | PathTrace `Prelude.elem` opts = PathTracer
              | otherwise = RayTrace
          directPhotonMapVisualisation = DirectPhotonMapVisualisation `Prelude.elem` opts
          enableIrradianceCache = IrradianceCaching `Prelude.elem` opts
@@ -140,6 +145,9 @@
 
   -- Display message about irradiance cache
   Prelude.putStrLn (if useIrradianceCache renderSettings then "Irradiance caching enabled" else "Irrradiance caching disabled")
+
+  -- Display message about path tracing
+  Prelude.putStrLn (if PathTrace `Prelude.elem` opts then "Path tracer enabled" else "Path tracer disabled")
 
   -- Render the image
   let renderSettings' = renderSettings { lights = lights' }
diff --git a/app/src/Material.hs b/app/src/Material.hs
--- a/app/src/Material.hs
+++ b/app/src/Material.hs
@@ -8,6 +8,7 @@
 data Material = Material { ambient :: {-# UNPACK #-} !Colour, 
                            diffuse :: {-# UNPACK #-} !Colour, 
                            specular :: {-# UNPACK #-} !Colour, 
+                           emission :: {-# UNPACK #-} !Colour, 
                            specularPower :: {-# UNPACK #-} !Double,
                            reflectivity :: {-# UNPACK #-} !Double,
                            transmit :: {-# UNPACK #-} !Double,
diff --git a/app/src/Matrix.hs b/app/src/Matrix.hs
--- a/app/src/Matrix.hs
+++ b/app/src/Matrix.hs
@@ -1,6 +1,5 @@
 -- 4D Matrix Library
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 
 module Matrix where
 
diff --git a/app/src/Misc.hs b/app/src/Misc.hs
--- a/app/src/Misc.hs
+++ b/app/src/Misc.hs
@@ -6,6 +6,8 @@
 import GHC.Prim
 import GHC.Types
 import Data.List
+import Control.Parallel.Strategies
+import Control.Monad.State
 
 degreesToRadians :: Double -> Double
 degreesToRadians x = x * pi / 180
@@ -27,5 +29,34 @@
       !(D# !value) = Prelude.max 0 (Prelude.min (D# x) 1)
 
 harmonicMean :: (Num t, Fractional t) => [t] -> t
-harmonicMean (x:xs) = fromIntegral (length (x:xs)) / foldl' (\a b -> b + 1 / a) 0 (x:xs)
+harmonicMean array@(_:_) = fromIntegral (length array) / foldl' (\a b -> b + 1 / a) 0 array
 harmonicMean [] = 0
+
+-- This performs a map, and passes through the state of the completed operation to the next recursion
+-- Couldn't work out the equivalent using the state monad etc
+mapS :: (a -> s -> (b, s)) -> [a] -> s -> ([b], s)
+mapS f z s = mapS' z s []
+    where
+      mapS' !(x:xs) !st !acc = seq (result, st') $ mapS' xs st' (result : acc)
+          where (!result, !st') = f x st `using` rseq
+      mapS' [] !st !acc = (acc, st)
+
+-- Map over a list, passing state from one to the next with the state monad
+mapWithState :: [a] -> s -> (a -> State s b) -> ([b], s)
+mapWithState arr s f = mapWithState' arr s []
+    where
+      mapWithState' (x:xs) st acc = mapWithState' xs st' (result : acc)
+          where
+            (result, st') = runState (f x) st
+      mapWithState' [] st acc = (acc, st)
+
+-- Zip over two lists, passing state from one to the next with the state monad
+zipWithState :: (a -> b -> State s c) -> [a] -> [b] -> s -> ([c], s)
+zipWithState f arr1 arr2 s = zipWithState' arr1 arr2 s []
+    where
+      zipWithState' (x:xs) (y:ys) st acc = zipWithState' xs ys st' (result : acc)
+          where
+            (result, st') = runState (f x y) st
+      zipWithState' (_:_) [] _ _ = error "Lists are of a different size - unhandled case!"
+      zipWithState' [] (_:_) _ _ = error "Lists are of a different size - unhandled case!"
+      zipWithState' [] [] st acc = (acc, st)
diff --git a/app/src/Octree.hs b/app/src/Octree.hs
--- a/app/src/Octree.hs
+++ b/app/src/Octree.hs
@@ -1,13 +1,12 @@
 -- This is a module for constructing bounding volume hierarchies using an octree approach
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 
 module Octree(generateSceneGraphUsingOctree, generateOctreeBoxList, OctTree(OctTreeNode, OctTreeLeaf, OctTreeDummy), create, Octree.insert, gather) where
 
 import Vector
 import Primitive
 import BoundingBox
-import Control.Parallel.Strategies
+import Misc
 
 data OctTree a = OctTreeDummy !AABB
                | OctTreeNode !AABB [OctTree a]
@@ -26,16 +25,6 @@
 
 create :: AABB -> OctTree a
 create box = OctTreeNode box $ map OctTreeDummy (generateOctreeBoxList box)
-
--- This performs a map, and passes through the state of the completed operation to the next recursion
--- Couldn't work out the equivalent using the state monad etc
-mapS :: (a -> s -> (b, s)) -> [a] -> s -> ([b], s)
-mapS f xs state = mapS' f xs state []
-
-mapS' :: (a -> s -> (b, s)) -> [a] -> s -> [b] -> ([b], s)
-mapS' !f !(x:xs) !state !acc = seq (result, state') $ mapS' f xs state' (result : acc)
-    where (!result, !state') = f x state `using` rseq
-mapS' _ [] !state !acc = (acc, state)
 
 -- Insert into an octree
 insert :: Vector -> a -> OctTree a -> OctTree a
diff --git a/app/src/PhotonMap.hs b/app/src/PhotonMap.hs
--- a/app/src/PhotonMap.hs
+++ b/app/src/PhotonMap.hs
@@ -3,6 +3,7 @@
 
 module PhotonMap(buildPhotonMap, PhotonMap(photonList), irradiance, PhotonMapContext(PhotonMapContext)) where
 
+import PolymorphicNum
 import {-# SOURCE #-} Light hiding (position)
 import Vector
 import Distribution
@@ -22,6 +23,7 @@
 import System.Random.Mersenne.Pure64
 import Data.List hiding (union, insert)
 import Primitive
+import RussianRoulette
 
 type GeneratorState = State PureMT
 
@@ -39,11 +41,10 @@
 data PhotonMap = PhotonMap { photonList :: [Photon],
                              photonMapTree :: PhotonMapTree } deriving(Show, Eq)
 
-data PhotonChoice = DiffuseReflect | SpecularReflect | Absorb
-
 instance NFData Photon where
     rnf (Photon power' posDir') = rnf power' `seq` rnf posDir'
 
+-- TODO - Sort this out!
 seedToRefactor :: Int
 seedToRefactor = 12345
 
@@ -54,27 +55,18 @@
 emitPhotons :: Light -> Int -> [(Position, Direction, PureMT, Colour)]
 emitPhotons (PointLight (CommonLightData lightPower True) pos _) numPhotons = zipWith (\dir num -> (pos, dir, pureMT (fromIntegral num), flux)) (generatePointsOnSphere numPhotons 1 seedToRefactor) [1..numPhotons]
     where
-      flux = lightPower Colour.<*> (1.0 / fromIntegral numPhotons)
+      flux = lightPower <*> ((1.0 / fromIntegral numPhotons) :: Double)
 emitPhotons (QuadLight (CommonLightData lightPower True) corner _ du dv) numPhotons = zipWith3 (\pos dir num -> (pos, transformDir dir tanSpace, pureMT (fromIntegral num), flux)) randomPoints randomDirs [1..numPhotons]
     where
       randomPoints = generatePointsOnQuad corner du dv numPhotons seedToRefactor
       randomDirs = generatePointsOnHemisphere numPhotons 1 (seedToRefactor * 10)
       area =  Vector.magnitude (du `cross` dv)
-      flux = lightPower Colour.<*> (area / fromIntegral numPhotons)
+      flux = lightPower <*> (area / fromIntegral numPhotons)
       tanSpace = (normalise du, normalise dv, normalise (du `cross` dv))
 emitPhotons _ _ = []
 
--- Compute russian roulette coefficients
-russianRouletteCoefficients :: Material -> (Double, Double)
-russianRouletteCoefficients mat = (diffuseP, specularP)
-    where
-      (Colour diffuseR diffuseG diffuseB _) = Material.diffuse mat
-      (Colour specularR specularG specularB _) = Material.specular mat
-      diffuseP = (diffuseR + diffuseG + diffuseB) / 3
-      specularP = (specularR + specularG + specularB) / 3
-
 -- Decide what to do with a photon
-choosePhotonFate :: (Double, Double) -> GeneratorState PhotonChoice
+choosePhotonFate :: (Double, Double) -> GeneratorState RussianRouletteChoice
 choosePhotonFate (diffuseP, specularP) = do
   generator <- get
   let (p, generator') = randomDouble generator
@@ -85,26 +77,18 @@
   return $! result
 
 -- Compute new power for a photon
-computeNewPhotonPower :: PhotonChoice -> (Double, Double) -> Colour -> Material -> Colour
+computeNewPhotonPower :: RussianRouletteChoice -> (Double, Double) -> Colour -> Material -> Colour
 computeNewPhotonPower fate (diffuseP, specularP) photonPower mat = case fate of
-                                                                     DiffuseReflect -> photonPower * diffuse mat Colour.</> diffuseP
-                                                                     SpecularReflect -> photonPower * specular mat Colour.</> specularP
+                                                                     DiffuseReflect -> photonPower <*> diffuse mat </> diffuseP
+                                                                     SpecularReflect -> photonPower <*> specular mat </> specularP
                                                                      Absorb -> colBlack
 
--- Compute a new diffuse reflection in spherical co-ordinates
-generateUV :: GeneratorState (Double, Double)
-generateUV = do generator <- get
-                let (u, generator') = randomDouble generator
-                let (v, generator'') = randomDouble generator'
-                put generator''
-                return $! (u, v)
-
 -- Find a diffuse reflection direction in the hemisphere of the normal
 -- Realistic Image Synthesis Using Photon Mapping - Eq 2.24
 diffuseReflectionDirection :: PureMT -> TangentSpace -> (Direction, PureMT)
 diffuseReflectionDirection stdGen tanSpace = (transformDir dir tanSpace, stdGen')
     where
-      ((u, v), stdGen') = runState generateUV stdGen
+      ((u, v), stdGen') = runState randomUV stdGen
       theta = acos (sqrt u)
       phi = 2 * pi * v
       dir = sphericalToDirection theta phi
@@ -143,7 +127,7 @@
             tanSpace = primitiveTangentSpace (primitive obj) subId hitPosition obj
             normal = thr tanSpace
             hitPosition = pointAlongRay ray t
-            surfacePos = hitPosition + (normal Vector.<*> surfaceEpsilon)
+            surfacePos = hitPosition <+> normal <*> surfaceEpsilon
             brightnessEpsilon = 0.1
             storedPhoton = Photon photonPower (surfacePos, snd photonPosDir)
     where
@@ -175,9 +159,9 @@
 buildKDTree (x:[]) = PhotonMapLeaf x
 buildKDTree [] = error "buildKDTree [] should never get called"
 buildKDTree photons = let (boxMin, boxMax) = photonsBoundingBox photons
-                          axis = largestAxis (boxMax - boxMin)
-                          numPhotons = fromIntegral (length photons)
-                          photonsMedian = foldl' (\box photon -> box + (fst . posDir $ photon)) zeroVector photons Vector.</> numPhotons
+                          axis = largestAxis (boxMax <-> boxMin)
+                          numPhotons = (fromIntegral (length photons)) :: Double
+                          photonsMedian = foldl' (\box photon -> box <+> (fst . posDir $ photon)) zeroVector photons </> numPhotons
                           value = component photonsMedian axis
                           photonsGT = Prelude.filter (\p -> component ((fst . posDir) p) axis > value) photons
                           photonsLE = Prelude.filter (\p -> component ((fst . posDir) p) axis <= value) photons
@@ -238,19 +222,19 @@
 -- Return the contribution of a given photon, including a simple cos term to emulate BRDF plus the cone filter
 -- Cone filter is from Realistic Image Synthesis Using Photon Mapping p81
 photonContribution :: Double -> SurfaceLocation -> Photon -> Colour
-photonContribution kr (pos, (_, _, normal)) photon = power photon Colour.<*> ((Vector.negate normal `sdot3` (snd . posDir) photon) * weight)
+photonContribution kr (pos, (_, _, normal)) photon = power photon <*> ((Vector.negate normal `sdot3` (snd . posDir) photon) * weight)
     where
       weight = 1 - (pos `distance` (fst . posDir) photon) / (kr + 0.000000001) -- Add on an epsilon to prevent div0 in cone filter
 
 -- Find the overall contribution of a list of photons
 -- Radiance estimate algorithm from Realistic Image Synthesis Using Photon Mapping p81
 sumPhotonContribution :: Double -> Double -> SurfaceLocation -> [Photon] -> Colour
-sumPhotonContribution r k posTanSpace photons = foldl' (\y x -> y + photonContribution (k * r) posTanSpace x) colBlack photons Colour.<*> (1.0 / ((1.0 - 2.0 / (3.0 * k)) * pi * r * r))
+sumPhotonContribution r k posTanSpace photons = foldl' (\y x -> y <+> photonContribution (k * r) posTanSpace x) colBlack photons <*> (1.0 / ((1.0 - 2.0 / (3.0 * k)) * pi * r * r))
 
 -- Look up the resulting irradiance from the photon map at a given point
 -- Realistic Image Synthesis Using Photon Mapping, e7.6
-irradiance :: PhotonMap -> PhotonMapContext -> Material -> SurfaceLocation -> Colour
-irradiance photonMap photonMapContext mat posTanSpace = sumPhotonContribution r k posTanSpace gatheredPhotons * diffuse mat
+irradiance :: PhotonMap -> PhotonMapContext -> Material -> SurfaceLocation -> (Colour, Double)
+irradiance photonMap photonMapContext mat posTanSpace = (sumPhotonContribution r k posTanSpace gatheredPhotons <*> diffuse mat, harmonicMean $ map (\(GatheredPhoton dist _) -> sqrt dist) nearestPhotons)
     where
       r = photonGatherDistance photonMapContext
       maxPhotons
@@ -258,4 +242,5 @@
           | otherwise = maxGatherPhotons photonMapContext
       k = coneFilterK photonMapContext
       photonHeap = gatherPhotons (photonMapTree photonMap) (fst posTanSpace) (r * r) Data.Heap.empty maxPhotons
-      gatheredPhotons = map (\(GatheredPhoton _ photon) -> photon) (Data.Heap.take maxPhotons photonHeap)
+      nearestPhotons = Data.Heap.take maxPhotons photonHeap
+      gatheredPhotons = map (\(GatheredPhoton _ photon) -> photon) nearestPhotons
diff --git a/app/src/PhotonMap.hs-boot b/app/src/PhotonMap.hs-boot
--- a/app/src/PhotonMap.hs-boot
+++ b/app/src/PhotonMap.hs-boot
@@ -24,5 +24,5 @@
 data PhotonMap = PhotonMap { photonList :: [Photon],
                              photonMapTree :: PhotonMapTree }
 
-irradiance :: PhotonMap -> PhotonMapContext -> Material -> (Position, TangentSpace) -> Colour
+irradiance :: PhotonMap -> PhotonMapContext -> Material -> (Position, TangentSpace) -> (Colour, Double)
 buildPhotonMap :: SceneGraph -> [Light] -> Int -> (PhotonMap, [Light])
diff --git a/app/src/PolymorphicNum.hs b/app/src/PolymorphicNum.hs
new file mode 100644
--- /dev/null
+++ b/app/src/PolymorphicNum.hs
@@ -0,0 +1,17 @@
+-- Module for a generic typeclass to bind together my linear algebra maths - vectors, matrices - rather than instancing off Num
+
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+module PolymorphicNum where
+
+class PolymorphicNum a b c | a b -> c where
+    (<*>) :: a -> b -> c
+    (</>) :: a -> b -> c
+    (<->) :: a -> b -> c
+    (<+>) :: a -> b -> c
+    infixl 7 <*>
+    infixl 7 </>
+    infixl 6 <+>
+    infixl 6 <->
+    
diff --git a/app/src/Primitive.hs b/app/src/Primitive.hs
--- a/app/src/Primitive.hs
+++ b/app/src/Primitive.hs
@@ -1,6 +1,5 @@
 -- Module for general primitives and intersections
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 {-# LANGUAGE UnboxedTuples #-}
 
 module Primitive (primitiveBoundingRadius, 
@@ -31,6 +30,7 @@
                   TangentSpace,
                   Vertex) where
 
+import PolymorphicNum
 import Ray
 import Vector
 import Material
@@ -62,15 +62,15 @@
 
 -- Surface normal for 3 points
 surfaceNormal :: Position -> Position -> Position -> Direction
-surfaceNormal !v1 !v2 !v3 = (v2 - v1) `cross` (v3 - v1)
+surfaceNormal !v1 !v2 !v3 = (v2 <-> v1) `cross` (v3 <-> v1)
 
 -- Make a plane
 makePlane :: Position -> Position -> Position -> Primitive
 makePlane !v1 !v2 !v3 = Plane (tangent, binormal, normal) (-(v1 `dot3` normal))
     where 
       !normal = normalise (surfaceNormal v1 v2 v3)
-      !tangent = normalise (v2 - v1)
-      !binormal = normalise (v3 - v1)
+      !tangent = normalise (v2 <-> v1)
+      !binormal = normalise (v3 <-> v1)
 
 -- -------------------------------------------------------------------------------------------------------------------------------------------------------------------
 -- Triangle base functionality
@@ -82,7 +82,7 @@
           newTanSpace = planeTangentSpace newPlane
           verts = map (\v -> Vertex v zeroVector newTanSpace) [v1, v2, v3]
           edgeVertices = [v1, v2, v3]
-          edges = map normalise [v2 - v1, v3 - v2, v1 - v3]
+          edges = map normalise [v2 <-> v1, v3 <-> v2, v1 <-> v3]
           edgeNormals = map (\edge -> normalise $ thr newTanSpace `cross` edge) edges
           -- TODO - The two vectors passed here are just dummies but they can fairly easily be derived
           newHalfPlanes = zipWith (\edgeNormal edgeVertex -> Plane (Vector 1 0 0 1, Vector 0 1 0 1, edgeNormal) (-(edgeNormal `dot3` edgeVertex))) edgeNormals edgeVertices
@@ -151,9 +151,9 @@
 interpolatedTangentSpace :: Triangle -> Double -> Double -> Double -> TangentSpace
 interpolatedTangentSpace !triangle !triAlpha !triBeta !triGamma = (tangent, binormal, normal)
     where [!(tan1, bi1, norm1), !(tan2, bi2, norm2), !(tan3, bi3, norm3)] = map vertTangentSpace (vertices triangle)
-          tangent = normalise $ (tan1 <*> triAlpha) + (tan2 <*> triBeta) + (tan3 <*> triGamma)
-          binormal = normalise $ (bi1 <*> triAlpha) + (bi2 <*> triBeta) + (bi3 <*> triGamma)
-          normal = normalise $ (norm1 <*> triAlpha) + (norm2 <*> triBeta) + (norm3 <*> triGamma)
+          tangent = normalise $ tan1 <*> triAlpha <+> tan2 <*> triBeta <+> tan3 <*> triGamma
+          binormal = normalise $ bi1 <*> triAlpha <+> bi2 <*> triBeta <+> bi3 <*> triGamma
+          normal = normalise $ norm1 <*> triAlpha <+> norm2 <*> triBeta <+> norm3 <*> triGamma
 
 -- -------------------------------------------------------------------------------------------------------------------------------------------------------------------
 -- Family of intersection functions
@@ -167,7 +167,7 @@
     | root2 >= 0 && root2 <= rayLen = Just (root2, 0)
     | otherwise = Nothing
     where 
-      !delta = rayOrg - getCentre obj
+      !delta = rayOrg <-> getCentre obj
       !b = 2 * (delta `dot3` rayDir)
       !c = (delta `dot3` delta) - sphereRadius ** 2
       !discriminant = b ** 2 - 4 * c -- A is 1 because the ray direction is normalised
@@ -198,7 +198,7 @@
     where
       tangent = Vector 1 0 0 0 -- This is clearly incorrect - fix this later!
       binormal = Vector 0 1 0 0
-      normal = (intersectionPoint - getCentre obj) <*> (1 / sphereRadius)
+      normal = (intersectionPoint <-> getCentre obj) <*> (1 / sphereRadius)
 
 primitiveTangentSpace (Plane !planeNormal _) _ _ _ = planeNormal
 primitiveTangentSpace (TriangleMesh !tris) !triId !intersectionPoint _ = interpolatedTangentSpace triangle triAlpha triBeta triGamma
@@ -227,8 +227,8 @@
 primitiveBoundingBox :: Primitive -> Object -> Maybe AABB
 primitiveBoundingBox (Sphere sphereRadius) obj = Just (boxMin, boxMax)
     where
-      boxMin = getCentre obj - Vector sphereRadius sphereRadius sphereRadius 0
-      boxMax = getCentre obj + Vector sphereRadius sphereRadius sphereRadius 0
+      boxMin = getCentre obj <-> Vector sphereRadius sphereRadius sphereRadius 0
+      boxMax = getCentre obj <+> Vector sphereRadius sphereRadius sphereRadius 0
 primitiveBoundingBox (Plane _ _) _ = Nothing
 primitiveBoundingBox (TriangleMesh tris) obj = Just $ triangleListBoundingBox initialInvalidBox (transform obj) tris
 
@@ -290,7 +290,7 @@
     | root2 >= 0 && root2 <= rayLen = Just root2
     | otherwise = Nothing
     where 
-      !delta = rayOrg - centre
+      !delta = rayOrg <-> centre
       !b = 2.0 * (delta `dot3` rayDir)
       !c = (delta `dot3` delta) - rad**2
       !discriminant = b**2 - 4 * c -- A is 1 because the ray direction is normalised
diff --git a/app/src/Ray.hs b/app/src/Ray.hs
--- a/app/src/Ray.hs
+++ b/app/src/Ray.hs
@@ -1,9 +1,9 @@
 -- Module for handling rays in a raytracer
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 
 module Ray where
 
+import PolymorphicNum
 import Vector
 
 -- For now, we're sticking to Doubles
@@ -11,7 +11,7 @@
 
 -- Make a ray given the start and end position
 rayWithPoints :: Position -> Position -> Ray
-rayWithPoints !start !end = Ray start (normalise (end - start)) (end `distance` start)
+rayWithPoints !start !end = Ray start (normalise (end <-> start)) (end `distance` start)
 
 rayWithDirection :: Position -> Direction -> Double -> Ray
 rayWithDirection !start !dir !rayLen = Ray start dir rayLen
diff --git a/app/src/RayTrace.hs b/app/src/RayTrace.hs
--- a/app/src/RayTrace.hs
+++ b/app/src/RayTrace.hs
@@ -1,9 +1,9 @@
 -- The module where all the tracing actually happens
 {-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE MagicHash #-}
 
-module RayTrace (rayTraceImage, findNearestIntersection, findAnyIntersection, GlobalIlluminationFunc) where
+module RayTrace (rayTraceImage, pathTraceImage, findNearestIntersection, findAnyIntersection, GlobalIlluminationFunc) where
 
+import PolymorphicNum
 import Vector
 import {-# SOURCE #-} Light
 import Primitive
@@ -20,6 +20,9 @@
 import IrradianceCache
 import Control.Monad.State
 import RenderContext
+import System.Random.Mersenne.Pure64
+import Data.List
+import RussianRoulette
 
 -- Intersect a ray against a sphere tree
 intersectSphereTree :: [SphereTreeNode] -> Ray -> Maybe (Object, Double, Int) -> Maybe (Object, Double, Int)
@@ -102,14 +105,11 @@
 -- Accumulate the contributions of the lights
 lightSurface :: [Light] -> Colour -> RenderContext -> SurfaceLocation -> Material -> Vector -> Colour
 lightSurface (x:xs) !acc renderContext !posTanSpace !objMaterial !viewDirection 
-    = let result = acc + applyLight (sceneGraph renderContext) posTanSpace objMaterial viewDirection x
+    = let result = acc <+> emissive <+> applyLight (sceneGraph renderContext) posTanSpace objMaterial viewDirection x
+          emissive = emission objMaterial
       in seq result (lightSurface xs result renderContext posTanSpace objMaterial viewDirection)
 lightSurface [] !acc _ _ _ _ = acc
 
--- Magic number for the usable radius of an irradaiance cache sample
-irrCacheSampleRadius :: Double
-irrCacheSampleRadius = 10
-
 -- Abstraction to permit different GI calculations
 type GlobalIlluminationFunc = (SurfaceLocation -> IrradianceCache -> Object -> RenderContext -> (Colour, IrradianceCache))
 
@@ -119,10 +119,10 @@
     case renderMode renderContext of
       PhotonMapper -> if useIrradianceCache renderContext
                       then query irrCache surfaceLocation irradiance'
-                      else (irradiance photonMap (photonMapContext renderContext) (material obj) surfaceLocation, irrCache)
+                      else (fst $ irradiance photonMap (photonMapContext renderContext) (material obj) surfaceLocation, irrCache)
       _ -> undefined -- Shouldn't hit this path...
     where
-      irradiance' x = (irradiance photonMap (photonMapContext renderContext) (material obj) x, irrCacheSampleRadius)
+      irradiance' x = irradiance photonMap (photonMapContext renderContext) (material obj) x
 photonMapGlobalIllumination _ _ irrCache _ _ = (colBlack, irrCache)
 
 -- Null GI
@@ -151,7 +151,7 @@
           let !intersectionPoint = pointAlongRay ray intersectionDistance
           let !tanSpace = primitiveTangentSpace (primitive obj) hitId intersectionPoint obj
           let !(!surfaceIrradiance, !newIrrCache) = calculateGI renderContext photonMap (intersectionPoint, tanSpace) irrCache obj renderContext
-          -- TODO - Need to plug irradiance values into surface shading more correctly
+          -- TODO - Need to plug irradiance values into shader model correctly
           let resultColour = lightSurface (lights renderContext) surfaceIrradiance renderContext (intersectionPoint, tanSpace) (material obj) viewDir
           put newIrrCache
           return $! resultColour
@@ -166,7 +166,7 @@
           let !tanSpace = primitiveTangentSpace (primitive obj) hitId intersectionPoint obj
           let !normal = thr tanSpace
           let !incoming = Vector.negate $ direction ray
-          -- TODO - Need to plug irradiance values into surface shading more correctly
+          -- TODO - Need to plug irradiance values into shader model correctly
 
           -- Evaluate result from irradiance cache
           irrCache <- get
@@ -199,58 +199,152 @@
           put irrCache'''
 
           -- Final colour combine
-          return $! (surfaceShading + (reflection Colour.<*> shine) + (refraction Colour.<*> transmittance))
+          return $! (surfaceShading <+> reflection <*> shine <+> refraction <*> transmittance)
               where
                 enteringObject !incoming !normal = incoming `dot3` normal > 0
 
--- This function converts a pixel co-ordinate to a direction of the ray
-makeRayDirection :: Int -> Int -> Camera -> (Int, Int) -> Vector
-makeRayDirection !renderWidth !renderHeight !camera (x, y) =
-    let x' = (fromIntegral x / fromIntegral renderWidth) * 2.0 - 1.0
-        y' = (fromIntegral y / fromIntegral renderHeight) * 2.0 - 1.0
-        fov = 0.5 * fieldOfView camera
-        fovX = tan (degreesToRadians fov)
-        fovY = -tan (degreesToRadians fov)
-        aspectRatio = fromIntegral renderWidth / fromIntegral renderHeight
-        !dirX = fovX * x'
-        !dirY = fovY * (-y') / aspectRatio
-        rayDir = normalise (Vector dirX dirY 1 0)
-    in normalise $ transformVector (worldToCamera camera) rayDir
-
 -- Trace a list of distributed samples with tail recursion
-traceDistributedSample :: RenderContext -> Colour -> [Position] -> Maybe PhotonMap -> (Position, Direction) -> Double -> RayTraceState
-traceDistributedSample renderContext !acc (x:xs) photonMap !eyeViewDir !sampleWeighting = 
+rayTracePixelSample :: RenderContext -> Colour -> [Position] -> Maybe PhotonMap -> (Position, Direction) -> Double -> RayTraceState
+rayTracePixelSample renderContext !acc (x:xs) photonMap !eyeViewDir !sampleWeighting = 
     do
       irrCache <- get
-      let !dofFocalDistance = depthOfFieldFocalDistance renderContext
-      let jitteredRayPosition jitter = fst eyeViewDir + jitter
+      let dofFocalDistance = depthOfFieldFocalDistance renderContext
+      let jitteredRayPosition jitter = fst eyeViewDir <+> jitter
       let jitteredRayDirection jitter = normalise $ madd jitter (snd eyeViewDir) dofFocalDistance
       let (sampleColour, irrCache') = runState (traceRay renderContext photonMap (rayWithDirection (jitteredRayPosition x) (jitteredRayDirection x) 100000.0) (maximumRayDepth renderContext) (snd eyeViewDir) 1 1) irrCache
-      let result = (sampleColour Colour.<*> sampleWeighting) + acc
+      let result = sampleColour <*> sampleWeighting <+> acc
       put irrCache'
-      let (col, irrCache'') = runState (traceDistributedSample renderContext result xs photonMap eyeViewDir sampleWeighting) irrCache'
+      let (col, irrCache'') = runState (rayTracePixelSample renderContext result xs photonMap eyeViewDir sampleWeighting) irrCache'
       put irrCache''
       return $! col
-traceDistributedSample _ !acc [] _ _ _ = return $! acc
+rayTracePixelSample _ !acc [] _ _ _ = return $! acc
 
--- Need to remove hard coded constants of 8 here
 -- This traces for a given pixel (x, y)
-tracePixel :: RenderContext -> Position -> Maybe PhotonMap -> Direction -> RayTraceState
-tracePixel renderContext !eye photonMap !viewDirection = do
+rayTracePixel :: RenderContext -> Position -> Maybe PhotonMap -> Direction -> RayTraceState
+rayTracePixel renderContext eye photonMap viewDirection = do
   irrCache <- get
-  let !distributedPositions = generatePointsOnSphere (numDistribSamples renderContext) (rayOriginDistribution renderContext) 12345
-  let (!pixelColour, irrCache') = runState (traceDistributedSample renderContext colBlack distributedPositions photonMap (eye, viewDirection) (1.0 / (fromIntegral . numDistribSamples $ renderContext))) irrCache
+  let distributedPositions = generatePointsOnSphere (numDistribSamples renderContext) (rayOriginDistribution renderContext) 12345 -- TODO Fix this magic number!
+  let (pixelColour, irrCache') = runState (rayTracePixelSample renderContext colBlack distributedPositions photonMap (eye, viewDirection) (1.0 / (fromIntegral . numDistribSamples $ renderContext))) irrCache
   put irrCache'
   return $! pixelColour
 
+-- This function converts a pixel co-ordinate to a direction of the ray
+makeRayDirection :: Int -> Int -> Camera -> (Double, Double) -> Vector
+makeRayDirection !renderWidth !renderHeight camera (x, y) =
+    let !x' = (x / fromIntegral renderWidth) * 2.0 - 1.0
+        !y' = (y / fromIntegral renderHeight) * 2.0 - 1.0
+        !fov = 0.5 * fieldOfView camera
+        !fovX = tan (degreesToRadians fov)
+        !fovY = -tan (degreesToRadians fov)
+        !aspectRatio = fromIntegral renderWidth / fromIntegral renderHeight
+        !dirX = fovX * x'
+        !dirY = fovY * (-y') / aspectRatio
+        !rayDir = normalise (Vector dirX dirY 1 0)
+    in normalise $ transformVector (worldToCamera camera) rayDir
+
 -- Generate a list of colours which contains a raytraced image. In parallel
 rayTraceImage :: RenderContext -> Camera -> Int -> Int -> Maybe PhotonMap -> [Colour]
 rayTraceImage renderContext camera renderWidth renderHeight photonMap = tracePixelPassingState rayDirections irrCache `using` parListChunk 256 rdeepseq
-    where !rayDirections = [makeRayDirection renderWidth renderHeight camera (x, y) | y <- [0..(renderHeight - 1)], x <- [0..(renderWidth - 1)]]
-          !eyePosition = Camera.position camera
+    where rayDirections = [makeRayDirection renderWidth renderHeight camera (fromIntegral x, fromIntegral y) | y <- [0..(renderHeight - 1)], x <- [0..(renderWidth - 1)]]
+          eyePosition = Camera.position camera
           irrCache = initialiseCache (sceneGraph renderContext)
           -- This function is the equivalent to map, but it passes the ending state of one invocation to the next invocation
-          tracePixelPassingState (x:xs) st = result : tracePixelPassingState xs st'
+          -- I'm using this rather than my mapWithState routine because I'm concerned to do so may break the parallelism of parListChunk
+          tracePixelPassingState !(x:xs) !st = result : tracePixelPassingState xs st'
               where
-                (!result, !st') = runState (tracePixel renderContext eyePosition photonMap x) st
+                (!result, !st') = runState (rayTracePixel renderContext eyePosition photonMap x) st
           tracePixelPassingState [] _ = []
+
+type PathTraceState = State PureMT Colour
+pathTrace :: RenderContext -> Ray -> Int -> Direction -> Double -> Colour -> PathTraceState
+
+-- General case
+pathTrace renderContext !ray depth !viewDir !currentIOR !weight =
+    case findNearestIntersection (sceneGraph renderContext) ray of
+        Nothing -> return $! colBlack
+        Just (obj, intersectionDistance, hitId) ->           
+            -- Evaluate surface-location specific things such as shader results
+            let intersectionPoint = pointAlongRay ray intersectionDistance
+                tanSpace = primitiveTangentSpace (primitive obj) hitId intersectionPoint obj
+                normal = thr tanSpace
+                incoming = Vector.negate $ direction ray
+
+                -- Thunk for emitted light
+                emittedLight = (emission . material) obj
+            
+                -- Compute radiance at this point
+                radiance = lightSurface (lights renderContext) colBlack renderContext (intersectionPoint, tanSpace) (material obj) viewDir
+
+                (diffuseP, specularP) = russianRouletteCoefficients (material obj)
+            in do
+              -- TODO - Need to evaluate the shader model here!
+
+              -- Randomly decide the fate at this intersection
+              gen <- get
+              let (p, gen') = randomDouble gen
+              put gen'
+              let interaction | p < diffuseP = DiffuseReflect
+                              | p < (diffuseP + specularP) = SpecularReflect
+                              | otherwise = Absorb
+
+              -- Thunk for reflected light 
+              let (randomDir, gen'') = generatePointOnHemisphere gen' 1
+              put gen''
+
+              -- Set up expression for reflected light
+              let reflectedDir 
+                      | interaction == DiffuseReflect = transformDir randomDir tanSpace
+                      | otherwise = incoming `reflect` normal
+              let ray' = rayWithDirection intersectionPoint reflectedDir (rayLength ray)
+              let weight' = (diffuse . material) obj <*> weight <*> (normal `sdot3` reflectedDir)
+              let (tracedPathColour, gen''') = runState (pathTrace renderContext ray' (depth + 1) viewDir currentIOR weight') gen''
+              let reflectedLight = tracedPathColour <*> weight'
+              put gen'''
+
+              -- Have to divide by probability to correctly account for that relative proportion of the domain
+              return $! case interaction of DiffuseReflect -> (emittedLight <+> radiance <+> reflectedLight) </> diffuseP
+                                            SpecularReflect -> (emittedLight <+> radiance <+> reflectedLight) </> (diffuseP + specularP)
+                                            Absorb -> (emittedLight <+> radiance) </> (1 - diffuseP - specularP)
+
+-- Path-trace a sub-sample
+pathTracePixelSample :: RenderContext -> Camera -> (Int, Int) -> (Int, Int) -> (Double, Double) -> (Double, Double) -> PathTraceState
+pathTracePixelSample renderContext camera xy (width, height) jitterUV stratUV = pathTrace renderContext ray 0 rayDirection 1 colWhite
+    where
+      jitteredX = (fromIntegral . fst) xy + fst stratUV + fst jitterUV
+      jitteredY = (fromIntegral . snd) xy + snd stratUV + snd jitterUV
+      rayDirection = makeRayDirection width height camera (jitteredX, jitteredY)
+      ray = rayWithDirection (Camera.position camera) rayDirection (farClip camera)
+
+-- Path-trace a pixel. Do stratified sub-sampling
+pathTracePixel :: RenderContext -> Camera -> (Int, Int) -> (Int, Int) -> PathTraceState
+pathTracePixel renderContext camera pixelCoords renderTargetSize =
+    do
+      -- Work out the jittered UV offsets
+      gen <- get
+      let (offsetUVs, gen') = runState (generateRandomUVs numPathTraceSamples) gen
+      let offsetUVs' = map (\(u, v) -> (u * du, v * dv)) offsetUVs
+      put gen'
+
+      -- Mung it all together
+      let (pixelSamples, gen'') = zipWithState (pathTracePixelSample renderContext camera pixelCoords renderTargetSize) offsetUVs' stratifiedCentres gen'
+      put gen''
+      return $! foldl' (\x y -> x <*> weight <+> y) colBlack pixelSamples
+    where
+      -- Total number of samples to take
+      numPathTraceSamplesRoot = 6 :: Int
+      numPathTraceSamples = numPathTraceSamplesRoot * numPathTraceSamplesRoot
+      weight = (1.0 :: Double) / fromIntegral numPathTraceSamples
+
+      -- Work out a set of stratified centres to jitter from
+      du = (1.0 :: Double) / fromIntegral numPathTraceSamplesRoot
+      dv = du
+      stratifiedCentres = [(fromIntegral x * du, fromIntegral y * dv) | y <- [0..(numPathTraceSamplesRoot - 1)], x <- [0..(numPathTraceSamplesRoot - 1)]]
+
+pathTraceImage :: RenderContext -> Camera -> Int -> Int -> [Colour]
+pathTraceImage renderContext camera renderWidth renderHeight = zipWith
+                                                               (\x y -> evalState (pathTracePixel renderContext camera x (renderWidth, renderHeight)) (pureMT y))
+                                                               [(x, y) | y <- [0..(renderHeight - 1)], x <- [0..(renderWidth - 1)]]
+                                                               [1..] 
+                                                               `using` parListChunk 256 rdeepseq
+
+-- TODO Re-unify the path tracer and ray tracer/photon map code paths as much as practical to ease maintenance (?)
diff --git a/app/src/RenderContext.hs b/app/src/RenderContext.hs
--- a/app/src/RenderContext.hs
+++ b/app/src/RenderContext.hs
@@ -8,6 +8,10 @@
 
 data RenderMode = RayTrace | PhotonMapper | PathTracer deriving (Show)
 
+data RenderMethodConfiguration = RenderMethodRayTrace
+                               | RenderMethodPhotonMap 
+                               | RenderMethodPathTracer deriving (Show)
+
 data RenderContext = RenderContext {
       numDistribSamples :: Int,
       sceneGraph :: SceneGraph,
diff --git a/app/src/RussianRoulette.hs b/app/src/RussianRoulette.hs
new file mode 100644
--- /dev/null
+++ b/app/src/RussianRoulette.hs
@@ -0,0 +1,15 @@
+-- Shared module for russian roulette across path tracer and photon mapper
+
+module RussianRoulette where
+
+import Material
+import Colour
+
+data RussianRouletteChoice = DiffuseReflect | SpecularReflect | Absorb deriving Eq
+
+-- Compute russian roulette coefficients
+russianRouletteCoefficients :: Material -> (Double, Double)
+russianRouletteCoefficients mat = (diffuseP, specularP)
+    where
+      diffuseP = (magnitude . Material.diffuse) mat
+      specularP = (magnitude . Material.specular) mat
diff --git a/app/src/SceneGraph.hs b/app/src/SceneGraph.hs
--- a/app/src/SceneGraph.hs
+++ b/app/src/SceneGraph.hs
@@ -1,8 +1,8 @@
 -- The graph structure holding the scene
-{-# LANGUAGE MagicHash #-}
 
 module SceneGraph (buildSceneGraph, SphereTreeNode(boundingRadius, boundingCentre, object, children), SceneGraph(root, infiniteObjects, finiteBox)) where
 
+import PolymorphicNum
 import Primitive
 import Vector
 import BoundingBox
@@ -13,13 +13,13 @@
 
 -- Find the mean of a collection of objects
 calculateMeanPosition' :: [Object] -> Vector -> Vector
-calculateMeanPosition' (obj : objects) acc = calculateMeanPosition' objects acc + getCentre obj
+calculateMeanPosition' (obj : objects) acc = calculateMeanPosition' objects acc <+> getCentre obj
 calculateMeanPosition' [] acc = acc
 
 calculateMeanPosition :: [Object] -> Vector
 calculateMeanPosition objects = setWTo1 (calculateMeanPosition' objects zeroVector </> len')
     where
-      !len' = fromIntegral (length objects)
+      !len' = (fromIntegral (length objects)) :: Double
 
 -- Find the overall bounding radius of a list of objects
 calculateBoundingRadius :: [Object] -> Vector -> Double
diff --git a/app/src/SceneGraphTest.hs b/app/src/SceneGraphTest.hs
deleted file mode 100644
--- a/app/src/SceneGraphTest.hs
+++ /dev/null
@@ -1,30 +0,0 @@
-module SceneGraphTest where
-
-import SceneGraph
-import Primitive
-import Material
-import Matrix
-import Shader
-import Colour
-
-testObjs1 = [
-    Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100)
-    ]
-
-testObjs2 = [
-    Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100),
-    Object (Sphere 50) (Material colGreen colGreen colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 0 100)
-    ]
-
-testObjs3 = [
-    Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100),
-    Object (Sphere 50) (Material colGreen colGreen colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 0 100),
-    Object (Sphere 50) (Material colBlue colBlue colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 200 100)
-    ]
-
-testObjs4 = [
-    Object (Sphere 50) (Material colRed colRed colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 0 100),
-    Object (Sphere 50) (Material colGreen colGreen colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 0 100),
-    Object (Sphere 50) (Material colBlue colBlue colWhite 50 0.5 0 iorAir NullShader) (translationMatrix (-200) 200 100),
-    Object (Sphere 50) (Material colYellow colYellow colWhite 50 0.5 0 iorAir NullShader) (translationMatrix 200 200 100)
-    ]
diff --git a/app/src/Shader.hs b/app/src/Shader.hs
--- a/app/src/Shader.hs
+++ b/app/src/Shader.hs
@@ -1,8 +1,8 @@
 -- Generic shaders to return colour and texture information
-{-# LANGUAGE MagicHash #-}
 
 module Shader where
 
+import PolymorphicNum
 import Vector
 import Colour
 import Misc
@@ -19,7 +19,7 @@
 
 -- Checked shaders
 evaluateDiffuse (CheckedShader checkScale checkColour1 checkColour2) position _ = 
-    let scaledPosition = checkScale * position
+    let scaledPosition = checkScale <*> position
         scaledX = round (vecX scaledPosition) :: Int
         scaledY = round (vecY scaledPosition) :: Int
         scaledZ = round (vecZ scaledPosition) :: Int
@@ -36,9 +36,9 @@
 evaluateAmbient = evaluateDiffuse
 
 -- New style shader interface
-shadePoint (CheckedShader checkScale checkColour1 checkColour2) (position, _) (ambient, diffuse, specular) = (ambient + diffuse + specular) * checkColour
+shadePoint (CheckedShader checkScale checkColour1 checkColour2) (position, _) (ambient, diffuse, specular) = (ambient <+> diffuse <+> specular) <*> checkColour
     where
-      scaledPosition = checkScale * position
+      scaledPosition = checkScale <*> position
       scaledX = round (vecX scaledPosition) :: Int
       scaledY = round (vecY scaledPosition) :: Int
       scaledZ = round (vecZ scaledPosition) :: Int
@@ -47,4 +47,4 @@
 shadePoint ShowNormalShader (_, norm) (_, _, _) = encodeNormal norm
 
 -- Default fallback
-shadePoint _ (_, _) (ambient, diffuse, specular) = ambient + diffuse + specular
+shadePoint _ (_, _) (ambient, diffuse, specular) = ambient <+> diffuse <+> specular
diff --git a/app/src/Tests/ColourTest.hs b/app/src/Tests/ColourTest.hs
--- a/app/src/Tests/ColourTest.hs
+++ b/app/src/Tests/ColourTest.hs
@@ -1,42 +1,28 @@
 module Tests.ColourTest where
 
+import PolymorphicNum
 import Colour
 import Test.HUnit
 
-test_Add = TestCase (assertEqual "Colour addition" expected (v1 + v2))
+test_Add = TestCase (assertEqual "Colour addition" expected (v1 <+> v2))
     where
       v1 = Colour 1 2 3 4
       v2 = Colour 10 20 30 40
       expected = Colour 11 22 33 44
 
-test_Sub = TestCase (assertEqual "Colour subtraction" expected (v1 - v2))
+test_Sub = TestCase (assertEqual "Colour subtraction" expected (v1 <-> v2))
     where
       v1 = Colour 10 20 30 40
       v2 = Colour 5 10 15 20
       expected = Colour 5 10 15 20
 
-test_Mul = TestCase (assertEqual "Colour multiplication" expected (v1 * v2))
+test_Mul = TestCase (assertEqual "Colour multiplication" expected (v1 <*> v2))
     where
       v1 = Colour 10 20 30 40
       v2 = Colour 5 10 15 20
       expected = Colour 50 200 450 800
 
-test_Abs = TestCase (assertEqual "Colour abs" expected (abs v1))
-    where
-      v1 = Colour 10 (-20) 0 (-40)
-      expected = Colour 10 20 0 40
-
-test_Signum = TestCase (assertEqual "Colour signum" expected (signum v1))
-    where
-      v1 = Colour 10 (-20) 0 (-40)
-      expected = Colour 1 (-1) 0 (-1)
-
-test_fromInteger = TestCase (assertEqual "Colour fromInteger" expected (signum v1))
-    where
-      v1 = fromInteger 5
-      expected = Colour 5 5 5 5
-
-test_Div = TestCase (assertEqual "Colour division" expected (v1 / v2))
+test_Div = TestCase (assertEqual "Colour division" expected (v1 </> v2))
     where
       v1 = Colour 10 20 30 40
       v2 = Colour 5 10 15 20
@@ -45,13 +31,13 @@
 test_AddScalar = TestCase (assertEqual "Colour add scalar" expected (v1 <+> k))
     where
       v1 = Colour 10 20 30 40
-      k = 1
+      k = 1 :: Double
       expected = Colour 11 21 31 41
 
 test_SubScalar = TestCase (assertEqual "Colour sub scalar" expected (v1 <-> k))
     where
       v1 = Colour 10 20 30 40
-      k = 1
+      k = 1 :: Double
       expected = Colour 9 19 29 39
 
 test_Clamp = TestCase (assertEqual "Colour clamp" expected (clamp v1))
@@ -73,8 +59,6 @@
                 TestLabel "Addition" test_Add,
                 TestLabel "Subtraction" test_Sub,
                 TestLabel "Multiplication" test_Mul,
-                TestLabel "Abs" test_Abs,
-                TestLabel "Signum" test_Signum,
                 TestLabel "Div" test_Div,
                 TestLabel "Add scalar" test_AddScalar,
                 TestLabel "Sub scalar" test_SubScalar,
diff --git a/app/src/Tests/VectorTest.hs b/app/src/Tests/VectorTest.hs
--- a/app/src/Tests/VectorTest.hs
+++ b/app/src/Tests/VectorTest.hs
@@ -1,41 +1,27 @@
 module Tests.VectorTest where
 
+import PolymorphicNum
 import Vector
 import Test.HUnit
 
-test_Add = TestCase (assertEqual "Vector addition" expectedResult (v1 + v2))
+test_Add = TestCase (assertEqual "Vector addition" expectedResult (v1 <+> v2))
     where
       v1 = Vector 1 2 3 4
       v2 = Vector 10 20 30 40
       expectedResult = Vector 11 22 33 44
 
-test_Sub = TestCase (assertEqual "Vector subtraction" expectedResult (v1 - v2))
+test_Sub = TestCase (assertEqual "Vector subtraction" expectedResult (v1 <-> v2))
     where
       v1 = Vector 10 20 30 40
       v2 = Vector 1 2 3 4
       expectedResult = Vector 9 18 27 36
 
-test_Mul = TestCase (assertEqual "Vector multiplication" expectedResult (v1 * v2))
+test_Mul = TestCase (assertEqual "Vector multiplication" expectedResult (v1 <*> v2))
     where
       v1 = Vector 1 0 2 3
       v2 = Vector 1 10 (-2) 3
       expectedResult = Vector 1 0 (-4) 9
 
-test_Negate = TestCase (assertEqual "Vector negation" expectedResult (Vector.negate v1))
-    where
-      v1 = Vector 1 (-2) 3 (-4)
-      expectedResult = Vector (-1) 2 (-3) 4
-
-test_Abs = TestCase (assertEqual "Vector abs" expectedResult (abs v1))
-    where
-      v1 = Vector 0 1 (-2) 3
-      expectedResult = Vector 0 1 2 3
-
-test_Signum = TestCase (assertEqual "Vector signum" expectedResult (signum v1))
-    where
-      v1 = Vector 1 (-1) 0 (-2)
-      expectedResult = Vector 1 (-1) 0 (-1)
-
 test_Madd = TestCase (assertEqual "Vector madd" expectedResult (madd pos dir k))
     where
       pos = Vector 1 2 3 1
@@ -43,16 +29,16 @@
       k = 10
       expectedResult = Vector 6 2 13 1
 
-test_ScalarMul = TestCase (assertEqual "Vector-scalar mul" expectedResult (vectorScalarMul vec k))
+test_ScalarMul = TestCase (assertEqual "Vector-scalar mul" expectedResult (vec <*> k))
     where
       vec = Vector 1 2 (-3) 1
-      k = 2
+      k = 2 :: Double
       expectedResult = Vector 2 4 (-6) 2
 
 test_ScalarDiv = TestCase (assertEqual "Vector-scalar div" expectedResult (vec </> k))
     where
       vec = Vector 10 20 (-30) 40
-      k = 2
+      k = 2 :: Double
       expectedResult = Vector 5 10 (-15) 20
 
 test_ScalarDot3 = TestCase (assertEqual "dot3" expectedResult (v1 `dot3` v2))
@@ -124,9 +110,6 @@
                 TestLabel "Addition" test_Add, 
                 TestLabel "Subtraction" test_Sub,
                 TestLabel "Multiplication" test_Mul,
-                TestLabel "Negation" test_Negate,
-                TestLabel "Abs" test_Abs,
-                TestLabel "Signum" test_Signum,
                 TestLabel "Madd" test_Madd,
                 TestLabel "Vector * scalar" test_ScalarMul,
                 TestLabel "Vector / scalar" test_ScalarDiv,
diff --git a/app/src/ToneMap.hs b/app/src/ToneMap.hs
--- a/app/src/ToneMap.hs
+++ b/app/src/ToneMap.hs
@@ -10,6 +10,7 @@
                imageAverageLogLuminance,
                imageAverageLuminance) where
 
+import PolymorphicNum
 import Colour
 
 -- x = x
@@ -30,9 +31,9 @@
 toneMapHejlBurgessDawson :: [Colour] -> [Colour]
 toneMapHejlBurgessDawson = map f
     where
-      f colour = (x * (x <*> 6.2 <+> 0.5)) / (x * (x <*> 6.2 <+> 1.7) <+> 0.06)
+      f colour = (x <*> (x <*> (6.2 :: Double) <+> (0.5 :: Double))) </> (x <*> (x <*> (6.2 :: Double) <+> (1.7 :: Double)) <+> (0.06:: Double))
           where
-            x = (\x' -> fold max x' 0) (colour <-> 0.004)
+            x = (\x' -> fold max x' 0) (colour <-> (0.004 :: Double))
 
 -- Apply a tone map operator
 toneMapImage :: ([Colour] -> [Colour]) -> [Colour] -> [Colour]
diff --git a/app/src/Vector.hs b/app/src/Vector.hs
--- a/app/src/Vector.hs
+++ b/app/src/Vector.hs
@@ -2,8 +2,12 @@
 
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE MagicHash #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
 
 module Vector where
+
+import Prelude
+import PolymorphicNum as L
 import Data.List
 import Misc
 import GHC.Prim
@@ -16,41 +20,54 @@
                        vecW :: {-# UNPACK #-} !Double } deriving (Ord, Eq)
 type Position = Vector
 type Direction = Vector
-type Normal = Direction
-type TangentSpace = (Normal, Normal, Normal)
+type Normal = Vector
+type TangentSpace = (Vector, Vector, Vector)
 type SurfaceLocation = (Position, TangentSpace)
 
-instance Num Vector where
-    {-# SPECIALIZE INLINE (+) :: Vector -> Vector -> Vector #-}
-    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) + (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x +## x') (D# $ y +## y') (D# $ z +## z') (D# $ w +## w')
-    {-# SPECIALIZE INLINE (-) :: Vector -> Vector -> Vector #-}
-    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) - (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x -## x') (D# $ y -## y') (D# $ z -## z') (D# $ w -## w')
-    {-# SPECIALIZE INLINE (*) :: Vector -> Vector -> Vector #-}
-    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) * (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x *## x') (D# $ y *## y') (D# $ z *## z') (D# $ w *## w')
-    abs (Vector !x !y !z !w) = Vector absX absY absZ absW
-        where
-          !absX = abs x
-          !absY = abs y
-          !absZ = abs z
-          !absW = abs w
-    signum (Vector !x !y !z !w) = Vector signumX signumY signumZ signumW
-        where
-          !signumX = signum x
-          !signumY = signum y
-          !signumZ = signum z
-          !signumW = signum w
+instance PolymorphicNum Vector Vector Vector where
+    {-# SPECIALIZE INLINE (<+>) :: Vector -> Vector -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) <+> (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x +## x') (D# $ y +## y') (D# $ z +## z') (D# $ w +## w')
+    {-# SPECIALIZE INLINE (<->) :: Vector -> Vector -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) <-> (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x -## x') (D# $ y -## y') (D# $ z -## z') (D# $ w -## w')
+    {-# SPECIALIZE INLINE (<*>) :: Vector -> Vector -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) <*> (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x *## x') (D# $ y *## y') (D# $ z *## z') (D# $ w *## w')
+    {-# SPECIALIZE INLINE (</>) :: Vector -> Vector -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) </> (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x /## x') (D# $ y /## y') (D# $ z /## z') (D# $ w /## w')
 
-    fromInteger x = Vector x' x' x' x'
-        where
-          !x' = fromInteger x
+instance PolymorphicNum Vector Double Vector where
+    {-# SPECIALIZE INLINE (<+>) :: Vector -> Double -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) <+> (!(D# k)) = Vector (D# $ x +## k) (D# $ y +## k) (D# $ z +## k) (D# $ w +## k)
+    {-# SPECIALIZE INLINE (<->) :: Vector -> Double -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) <-> (!(D# k)) = Vector (D# $ x -## k) (D# $ y -## k) (D# $ z -## k) (D# $ w -## k)
+    {-# SPECIALIZE INLINE (<*>) :: Vector -> Double -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) <*> (!(D# k)) = Vector (D# $ x *## k) (D# $ y *## k) (D# $ z *## k) (D# $ w *## k)
+    {-# SPECIALIZE INLINE (</>) :: Vector -> Double -> Vector #-}
+    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) </> (!(D# k)) = Vector (D# $ x /## k) (D# $ y /## k) (D# $ z /## k) (D# $ w /## k)
 
-instance Fractional Vector where
-    {-# SPECIALIZE INLINE (/) :: Vector -> Vector -> Vector #-}
-    (Vector !(D# x) !(D# y) !(D# z) !(D# w)) / (Vector !(D# x') !(D# y') !(D# z') !(D# w')) = Vector (D# $ x /## x') (D# $ y /## y') (D# $ z /## z') (D# $ w /## w')
-    fromRational x = Vector x' x' x' x'
-        where
-          !x' = fromRational x
+instance PolymorphicNum Double Vector Vector where
+    {-# SPECIALIZE INLINE (<+>) :: Double -> Vector -> Vector #-}
+    (!(D# k)) <+> (Vector !(D# x) !(D# y) !(D# z) !(D# w)) = Vector (D# $ k +## x) (D# $ k +## y) (D# $ k +## z) (D# $ k +## w)
+    {-# SPECIALIZE INLINE (<->) :: Double -> Vector -> Vector #-}
+    (!(D# k)) <-> (Vector !(D# x) !(D# y) !(D# z) !(D# w)) = Vector (D# $ k -## x) (D# $ k -## y) (D# $ k -## z) (D# $ k -## w)
+    {-# SPECIALIZE INLINE (<*>) :: Double -> Vector -> Vector #-}
+    (!(D# k)) <*> (Vector !(D# x) !(D# y) !(D# z) !(D# w)) = Vector (D# $ k *## x) (D# $ k *## y) (D# $ k *## z) (D# $ k *## w)
+    {-# SPECIALIZE INLINE (</>) :: Double -> Vector -> Vector #-}
+    (!(D# k)) </> (Vector !(D# x) !(D# y) !(D# z) !(D# w)) = Vector (D# $ k /## x) (D# $ k /## y) (D# $ k /## z) (D# $ k /## w)
 
+{-
+instance PolymorphicNum Vector (Num a) Vector where
+    (Vector x y z w) <+> k = Vector (x + k) (y + k) (z + k) (w + k)
+    (Vector x y z w) <-> k = Vector (x - k) (y - k) (z - k) (w - k)
+    (Vector x y z w) <*> k = Vector (x * k) (y * k) (z * k) (w * k)
+    (Vector x y z w) </> k = Vector (x / k) (y / k) (z / k) (w / k)
+
+instance PolymorphicNum (Num a) Vector Vector where
+    k <+> (Vector x y z w) = Vector (k + x) (k + y) (k + z) (k + w)
+    k <-> (Vector x y z w) = Vector (k - x) (k - y) (k - z) (k - w)
+    k <*> (Vector x y z w) = Vector (k * x) (k * y) (k * z) (k * w)
+    k </> (Vector x y z w) = Vector (k / x) (k / y) (k / z) (k / w)
+-}
+
 instance Show Vector where
     show (Vector !x !y !z !w) = "(" ++ show x ++ ", " ++ show y ++ ", " ++ show z ++ ", " ++ show w ++ ")"
 
@@ -106,23 +123,6 @@
 {-# SPECIALIZE INLINE Vector.negate :: Vector -> Vector #-}
 negate (Vector !x !y !z !w) = Vector (-x) (-y) (-z) (-w)
 
-vectorScalarMul :: Vector -> Double -> Vector
-{-# SPECIALIZE INLINE vectorScalarMul :: Vector -> Double -> Vector #-}
-(Vector !x !y !z !w) `vectorScalarMul` k = Vector (x * k) (y * k) (z * k) (w * k)
-
-infixl 7 <*>
-infixl 7 </>
---infixl 6 <+>
---infixl 6 <->
-
-(</>) :: Vector -> Double -> Vector
-{-# SPECIALIZE INLINE (</>) :: Vector -> Double -> Vector #-}
-a </> b = a `vectorScalarMul` (1 / b)
-
-(<*>) :: Vector -> Double -> Vector
-{-# SPECIALIZE INLINE (<*>) :: Vector -> Double -> Vector #-}
-a <*> b = a `vectorScalarMul` b
-
 dot3 :: Vector -> Vector -> Double
 {-# SPECIALIZE INLINE dot3 :: Vector -> Vector -> Double #-}
 (Vector !(D# x) !(D# y) !(D# z) _) `dot3` (Vector !(D# x') !(D# y') !(D# z') _) = D# $ (x *## x') +## (y *## y') +## (z *## z')
@@ -157,25 +157,25 @@
 
 normalise :: Direction -> Direction
 {-# SPECIALIZE INLINE normalise :: Direction -> Direction #-}
-normalise !a = setWTo0 (a `vectorScalarMul` (1 / magnitude a))
+normalise !a = setWTo0 (a </> magnitude a)
 
 distance :: Position -> Position -> Double
 {-# SPECIALIZE INLINE distance :: Position -> Position -> Double #-}
-distance !a !b = magnitude (a - b)
+distance !a !b = magnitude (a <-> b)
 
 distanceSq :: Position -> Position -> Double
 {-# SPECIALIZE INLINE distanceSq :: Position -> Position -> Double #-}
-distanceSq !a !b = magnitudeSq (a - b)
+distanceSq !a !b = magnitudeSq (a <-> b)
 
 reflect :: Direction -> Direction -> Direction
 {-# SPECIALIZE INLINE reflect :: Direction -> Direction -> Direction #-}
-reflect !incoming !normal = restoreOriginalW incoming $ (normal `vectorScalarMul` (2 * (normal `dot3` incoming))) - incoming
+reflect !incoming !normal = setWTo0 $ (normal <*> (2 * (normal `dot3` incoming))) <-> incoming
 
 refract :: Direction -> Direction -> Double -> Direction
 {-# SPECIALIZE INLINE refract :: Vector -> Vector -> Double -> Vector #-}
 refract !incoming !normal !eta
-    | cosTheta1 >## 0.0## = setWTo0 $ (l `vectorScalarMul` eta) + (normal `vectorScalarMul` D# (eta# *## cosTheta1 -## cosTheta2))
-    | otherwise = setWTo0 $ (l `vectorScalarMul` eta) + (normal `vectorScalarMul` D# (eta# *## cosTheta1 +## cosTheta2))
+    | cosTheta1 >## 0.0## = setWTo0 $ (l <*> eta) <+> (normal <*> D# (eta# *## cosTheta1 -## cosTheta2))
+    | otherwise = setWTo0 $ (l <*> eta) <+> (normal <*> D# (eta# *## cosTheta1 +## cosTheta2))
     where !(D# cosTheta1) = normal `dot3` incoming
           !cosTheta2 = sqrtDouble# (1.0## -## eta# **## 2.0## *## (1.0## -## cosTheta1 **## 2.0##))
           !l = Vector.negate incoming
@@ -230,4 +230,4 @@
 
 transformDir :: Direction -> TangentSpace -> Direction
 {-# SPECIALIZE INLINE transformDir :: Direction -> TangentSpace -> Direction #-}
-transformDir (Vector !x !y !z _) !(tangent, binormal, normal) = setWTo0 ((tangent <*> x) + (binormal <*> y) + (normal <*> z))
+transformDir (Vector !x !y !z _) !(tangent, binormal, normal) = setWTo0 $ tangent <*> x <+> binormal <*> y <+> normal <*> z
diff --git a/app/src/scripts/unit_test b/app/src/scripts/unit_test
--- a/app/src/scripts/unit_test
+++ b/app/src/scripts/unit_test
@@ -1,2 +1,2 @@
 #!/bin/sh
-runhaskell Tests/UnitTests.hs
+runhaskell -- -fobject-code Tests/UnitTests.hs
diff --git a/crocodile.cabal b/crocodile.cabal
--- a/crocodile.cabal
+++ b/crocodile.cabal
@@ -8,7 +8,7 @@
 -- The package version. See the Haskell package versioning policy
 -- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for
 -- standards guiding when and how versions should be incremented.
-Version:             0.1.1
+Version:             0.1.2
 
 stability:	     Experimental