This is a type safe API for https://github.com/guibou/streamray which hides the `V3` from `Linear` under a `newtype` qualified with a type phantom describing if it is a position, direction (normalized or not) or color. The typeclass `Add`, `Mul`, `Sub` describes the allowed operations.

Algebra.hs

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE TypeFamilies #-}

module Algebra
  ( pattern P,
    pattern N,
    pattern D,
    pattern C,
    Space (..),
    DirectionKind (..),
    V3,
    Add (..),
    Mul (..),
    Sub (..),
    normalize,
    dot,
  )
where

import qualified Linear

data DirectionKind = Normalized | 'NotNormalized

data Space = Position | Direction DirectionKind | Color

newtype V3 (k :: Space) = V3 (Linear.V3 Float)
    deriving (Show)

-- * Add
class Add a b where
  type AddResult a b
  (.+.) :: a -> b -> AddResult a b

infixl 6 .+.

type family AddResultV3 (a :: Space) (b :: Space) where
  AddResultV3 'Color 'Color = 'Color
  AddResultV3 'Position ('Direction k) = 'Position
  AddResultV3 ('Direction k) 'Position = 'Position
  AddResultV3 ('Direction k) ('Direction k') = 'Direction 'NotNormalized

unsafeLiftLinear :: (Linear.V3 Float -> Linear.V3 Float -> Linear.V3 Float) -> V3 k1 -> V3 k2 -> V3 k3
unsafeLiftLinear f (V3 x) (V3 y) = V3 (f x y)

instance Add (V3 a) (V3 b) where
  type AddResult (V3 a) (V3 b) = V3 (AddResultV3 a b)
  (.+.) = unsafeLiftLinear (+)

-- * Sub

class Sub a b where
  type SubResult a b
  (.-.) :: a -> b -> SubResult a b

infixl 6 .-.

type family SubResultV3 a b where
  SubResultV3 'Position ('Direction k) = 'Position
  SubResultV3 'Position 'Position = 'Direction 'NotNormalized
  SubResultV3 ('Direction k) ('Direction k') = 'Direction 'NotNormalized

instance Sub (V3 a) (V3 b) where
  type SubResult (V3 a) (V3 b) = V3 (SubResultV3 a b)
  (.-.) = unsafeLiftLinear (-)

-- * Mul

class Mul a b where
  type MulResult a b
  (.*.) :: a -> b -> MulResult a b

infixl 7 .*.

type family MulResultV3 a b where
  MulResultV3 'Color 'Color = 'Color
  MulResultV3 'Position ('Direction k) = 'Position
  MulResultV3 ('Direction k) ('Direction k') = 'Direction 'NotNormalized

instance Mul (V3 a) (V3 b) where
  type MulResult (V3 a) (V3 b) = V3 (MulResultV3 a b)
  (.*.) = unsafeLiftLinear (*)

instance Mul (V3 k) Float where
  type MulResult (V3 k) Float = V3 (MulResultScalarV3 k)
  V3 v .*. f = V3 (v Linear.^* f)

type family MulResultScalarV3 a where
  MulResultScalarV3 ('Direction k) = 'Direction 'NotNormalized
  MulResultScalarV3 'Color = 'Color

instance Mul Float (V3 k) where
  type MulResult Float (V3 k) = V3 (MulResultScalarV3 k)
  f .*. V3 v = V3 (f Linear.*^ v)

-- * Ctor and aliases

-- 'Position
pattern P :: Float -> Float -> Float -> V3 'Position
pattern P x y z = V3 (Linear.V3 x y z)

-- 'Direction
pattern D :: Float -> Float -> Float -> V3 ('Direction 'NotNormalized)
pattern D x y z = V3 (Linear.V3 x y z)

pattern N :: Float -> Float -> Float -> V3 ('Direction 'Normalized)
pattern N x y z <-
  V3 (Linear.V3 x y z)
  where
    N x y z = V3 (Linear.V3 x' y' z')
      where
        n = sqrt (x * x + y * y + z * z)
        x' = x / n
        y' = y / n
        z' = z / n

pattern C :: Float -> Float -> Float -> V3 'Color
pattern C x y z = V3 (Linear.V3 x y z)

-- functions

dot :: V3 ('Direction k) -> V3 ('Direction k') -> Float
dot (V3 x) (V3 y) = Linear.dot x y

normalize :: V3 ('Direction 'NotNormalized) -> V3 ('Direction 'Normalized)
normalize (V3 v) = V3 (Linear.normalize v)

Main.hs.diff

diff --git a/app/Main.hs b/app/Main.hs
index fc8dcf1..5c6ba2c 100644
--- a/app/Main.hs
+++ b/app/Main.hs
@@ -1,26 +1,29 @@
+{-# LANGUAGE DataKinds #-}
 {-# LANGUAGE MultiWayIf #-}
 {-# LANGUAGE NamedFieldPuns #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TupleSections #-}
+{-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE ViewPatterns #-}
 
 module Main where
 
+import Algebra
 import Codec.Picture
 import Data.List
 import Data.Maybe
 import Data.Ord (comparing)
-import Linear
 
 -- | This is a ray
 -- Any point X on the ray can be represented using X = Origin + t * Direction.
 data Ray = Ray
-  { origin :: V3 Float,
-    direction :: V3 Float
+  { origin :: V3 'Position,
+    direction :: V3 ('Direction 'Normalized)
   }
   deriving (Show)
 
 data Sphere = Sphere
-  { center :: V3 Float,
+  { center :: V3 'Position,
     radius :: Float
   }
   deriving (Show)
@@ -65,7 +68,7 @@ rayIntersectSphere Ray {origin, direction} Sphere {radius, center} =
       --
       -- >>> introduces some simplifications
 
-      oc = center - origin
+      oc = center .-. origin
       r2 = radius ** 2
       -- length(t * D - OC) ^ 2 = r2
       --
@@ -109,29 +112,29 @@ sphereRadius = 5000
 scene :: [Object]
 scene =
   [ Object
-      (Material (V3 1 1 0) Diffuse)
-      (Sphere (V3 (sphereRadius + 500) 250 0) sphereRadius), -- Right
+      (Material (C 1 1 0) Diffuse)
+      (Sphere (P (sphereRadius + 500) 250 0) sphereRadius), -- Right
     Object
-      (Material (V3 0 1 1) Diffuse)
-      (Sphere (V3 (- sphereRadius) 250 0) sphereRadius), -- Left
+      (Material (C 0 1 1) Diffuse)
+      (Sphere (P (- sphereRadius) 250 0) sphereRadius), -- Left
     Object
-      (Material (V3 1 1 1) Diffuse)
-      (Sphere (V3 250 (- sphereRadius) 0) sphereRadius), -- Top
+      (Material (C 1 1 1) Diffuse)
+      (Sphere (P 250 (- sphereRadius) 0) sphereRadius), -- Top
     Object
-      (Material (V3 1 1 1) Diffuse)
-      (Sphere (V3 250 (sphereRadius + 500) 0) sphereRadius), -- Bottom
+      (Material (C 1 1 1) Diffuse)
+      (Sphere (P 250 (sphereRadius + 500) 0) sphereRadius), -- Bottom
     Object
-      (Material (V3 1 1 1) Diffuse)
-      (Sphere (V3 250 250 (sphereRadius + 500)) sphereRadius), -- Back
+      (Material (C 1 1 1) Diffuse)
+      (Sphere (P 250 250 (sphereRadius + 500)) sphereRadius), -- Back
       -- Small sphere 1
     Object
-      (Material (V3 1 1 1) Diffuse)
-      (Sphere (V3 150 350 350) 80),
+      (Material (C 1 1 1) Diffuse)
+      (Sphere (P 150 350 350) 80),
     --
     -- Small sphere 1
     Object
-      (Material (V3 1 1 1) Diffuse)
-      (Sphere (V3 350 350 350) 80)
+      (Material (C 1 1 1) Diffuse)
+      (Sphere (P 350 350 350) 80)
   ]
 
 -- | Represents a object, with its shape and material
@@ -140,7 +143,7 @@ data Object = Object Material Sphere
 
 -- | Material, with albedo and behavior
 data Material
-  = Material (V3 Float) MaterialBehavior
+  = Material (V3 'Color) MaterialBehavior
   deriving (Show)
 
 data MaterialBehavior
@@ -152,21 +155,21 @@ data MaterialBehavior
     Mirror
   deriving (Show)
 
-lightPosition :: V3 Float
-lightPosition = V3 250 250 250
+lightPosition :: V3 'Position
+lightPosition = P 250 250 250
 
-lightEmission :: V3 Float
-lightEmission = V3 30000 30000 30000
+lightEmission :: V3 'Color
+lightEmission = C 30000 30000 30000
 
-(-->) :: Num a => a -> a -> a
-x --> y = y - x
+(-->) :: V3 'Position -> V3 'Position -> V3 ('Direction 'NotNormalized)
+x --> y = y .-. x
 
 -- | Returns the pixel color associated with a 'Ray'
 radiance :: Ray -> PixelRGBA8
 radiance ray = case rayIntersectObjets ray scene of
   Nothing -> PixelRGBA8 0 0 0 255
   Just (t, Object (Material albedo behavior) sphere) -> do
-    let x = origin ray + t *^ direction ray
+    let x = origin ray .+. t .*. direction ray
         directionToLight = x --> lightPosition
         normal = normalize (center sphere --> x)
 
@@ -175,16 +178,19 @@ radiance ray = case rayIntersectObjets ray scene of
 
         lightDistance2 = dot directionToLight directionToLight
 
-    tonemap (lightEmission * (coef *^ albedo))
+    tonemap (lightEmission .*. (coef .*. albedo))
 
 -- | Convert a light measure to a pixel value
-tonemap :: V3 Float -> PixelRGBA8
-tonemap v =
-  -- truncate converts to Word8
-  -- min/max clamps to the acceptable range
-  -- pow (1 / 2.2) is doing gamma correction
-  let V3 x y z = truncate . max 0 . min 255 <$> ((v ** (1 / 2.2)) * 255)
-   in PixelRGBA8 x y z 255
+tonemap :: V3 'Color -> PixelRGBA8
+tonemap v = PixelRGBA8 x y z 255
+  where
+    -- truncate converts to Word8
+    -- min/max clamps to the acceptable range
+    -- pow (1 / 2.2) is doing gamma correction
+    C (ftonemap -> x) (ftonemap -> y) (ftonemap -> z) = v
+
+ftonemap :: Float -> Pixel8
+ftonemap = truncate @Float @Pixel8 . max 0 . min 255 . (* 255) . (** (1 / 2.2))
 
 -- | Raytrace a 500x500 image
 -- This function is called for each pixel
@@ -194,11 +200,11 @@ raytrace (fromIntegral -> x) (fromIntegral -> y) = radiance ray
     -- Generate a ray in the XY plane and pointing in the Z direction
     coefOpening = 1.001
 
-    n = V3 x y 0
+    n = P x y 0
 
     -- n' is on the plane [-250:250]
-    n'@(V3 x' y' 0) = n - V3 250 250 0
-    f = V3 (coefOpening * x') (coefOpening * y') 1
+    n'@(P x' y' 0) = n .-. D 250 250 0
+    f = P (coefOpening * x') (coefOpening * y') 1
 
     d = normalize (n' --> f)
     ray = Ray n d