{-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances #-}

module Vector where

type R = Double

data V1 = V1 !R deriving (Show, Eq, Ord)
data V2 = V2 !R !R deriving (Show, Eq, Ord)
data V3 = V3 !R !R !R deriving (Show, Eq, Ord)
data V4 = V4 !R !R !R !R deriving (Show, Eq, Ord)

class V a where
  o       :: a
  (^+^)   :: a -> a -> a
  (^-^)   :: a -> a -> a
  (*^)    :: R -> a -> a
  (^*)    :: a -> R -> a
  (^/)    :: a -> R -> a
  dot     :: a -> a -> R

(|-|) :: V a => a -> a -> R
u |-| v = let d = u ^-^ v in d `dot` d

norm :: V a => a -> a
norm v = v ^/ sqrt (v`dot`v)

instance V V1 where
  o = V1 0
  V1 a ^+^ V1 x = V1 (a + x)
  V1 a ^-^ V1 x = V1 (a - x)
  k *^ V1 x = V1 (k * x)
  V1 a ^* k = V1 (a * k)
  V1 a ^/ k = V1 (a / k)
  V1 a `dot` V1 x = a * x

instance V V2 where
  o = V2 0 0
  V2 a b ^+^ V2 x y = V2 (a + x) (b + y)
  V2 a b ^-^ V2 x y = V2 (a - x) (b - y)
  k *^ V2 x y = V2 (k * x) (k * y)
  V2 a b ^* k = V2 (a * k) (b * k)
  V2 a b ^/ k = V2 (a / k) (b / k)
  V2 a b `dot` V2 x y = a * x + b * y

instance V V3 where
  o = V3 0 0 0
  V3 a b c ^+^ V3 x y z = V3 (a + x) (b + y) (c + z)
  V3 a b c ^-^ V3 x y z = V3 (a - x) (b - y) (c - z)
  k *^ V3 x y z = V3 (k * x) (k * y) (k * z)
  V3 a b c ^* k = V3 (a * k) (b * k) (c * k)
  V3 a b c ^/ k = V3 (a / k) (b / k) (c / k)
  V3 a b c `dot` V3 x y z = a * x + b * y + c * z

cross3 :: V3 -> V3 -> V3
cross3 (V3 a1 a2 a3) (V3 b1 b2 b3) =
  V3 (a2*b3-a3*b2) (a3*b1-a1*b3) (a1*b2-a2*b1)

instance V V4 where
  o = V4 0 0 0 0
  V4 a b c d ^+^ V4 x y z w = V4 (a + x) (b + y) (c + z) (d + w)
  V4 a b c d ^-^ V4 x y z w = V4 (a - x) (b - y) (c - z) (d - w)
  k *^ V4 x y z w = V4 (k * x) (k * y) (k * z) (k * w)
  V4 a b c d ^* k = V4 (a * k) (b * k) (c * k) (d * k)
  V4 a b c d ^/ k = V4 (a / k) (b / k) (c / k) (d / k)
  V4 a b c d `dot` V4 x y z w = a * x + b * y + c * z + d * w

data M1 = M1 !R deriving (Show, Eq, Ord)
data M2 = M2 !R !R  !R !R deriving (Show, Eq, Ord)
data M3 = M3 !R !R !R  !R !R !R  !R !R !R deriving (Show, Eq, Ord)
data M4 = M4 !R !R !R !R  !R !R !R !R  !R !R !R !R  !R !R !R !R deriving (Show, Eq, Ord)

class M a where
  i       :: a
  (^^+^^) :: a -> a -> a
  (^^-^^) :: a -> a -> a
  (^^*^^) :: a -> a -> a
  (*^^)   :: R -> a -> a
  (^^*)   :: a -> R -> a
  (^^/)   :: a -> R -> a
  mdot    :: a -> a -> R
  det     :: a -> R
  inv     :: a -> a

(||-||) :: M a => a -> a -> R
u ||-|| v = let d = u ^^-^^ v in d `mdot` d

instance M M1 where
  i = M1 1
  M1 a11 ^^+^^ M1 b11 = M1 (a11 + b11)
  M1 a11 ^^-^^ M1 b11 = M1 (a11 - b11)
  M1 a11 ^^*^^ M1 b11 = M1 (a11 * b11)
  a *^^ M1 b11 = M1 (a * b11)
  M1 a11 ^^* b = M1 (a11 * b)
  M1 a11 ^^/ b = M1 (a11 / b)
  M1 a11 `mdot` M1 b11 = (a11 * b11)
  det (M1 a11) = a11
  inv (M1 a11) = M1 (1/a11)

instance M M2 where
  i = M2 1 0 0 1
  M2 a11 a12 a21 a22 ^^+^^ M2 b11 b12 b21 b22 =
    M2 (a11 + b11) (a12 + b12) (a21 + b21) (a22 + b22)
  M2 a11 a12 a21 a22 ^^-^^ M2 b11 b12 b21 b22 =
    M2 (a11 - b11) (a12 - b12) (a21 - b21) (a22 - b22)
  M2 a11 a12 a21 a22 ^^*^^ M2 b11 b12 b21 b22 =
    M2 (a11*b11 + a12*b21) (a11*b12 + a12*b22) (a21*b11 + a22*b21) (a21*b12 + a22*b22)
  a *^^ M2 b11 b12 b21 b22 = M2 (a * b11) (a * b12) (a * b21) (a * b22)
  M2 a11 a12 a21 a22 ^^* b = M2 (a11 * b) (a12 * b) (a21 * b) (a22 * b)
  M2 a11 a12 a21 a22 ^^/ b = M2 (a11 / b) (a12 / b) (a21 / b) (a22 / b)
  M2 a11 a12 a21 a22 `mdot` M2 b11 b12 b21 b22 =
    a11*b11 + a12*b12 + a21*b21 + a22*b22
  det (M2 a11 a12 a21 a22) = a11 * a22 - a12 * a21
  inv a@(M2 a11 a12 a21 a22) = (M2 a22 (-a12) (-a21) a11) ^^/ det a

instance M M3 where
  i = M3 1 0 0 0 1 0 0 0 1
  M3 a11 a12 a13 a21 a22 a23 a31 a32 a33 ^^+^^ M3 b11 b12 b13 b21 b22 b23 b31 b32 b33 =
    M3 (a11 + b11) (a12 + b12) (a13 + b13) (a21 + b21) (a22 + b22) (a23 + b23) (a31 + b31) (a32 + b32) (a33 + b33)
  M3 a11 a12 a13 a21 a22 a23 a31 a32 a33 ^^-^^ M3 b11 b12 b13 b21 b22 b23 b31 b32 b33 =
    M3 (a11 - b11) (a12 - b12) (a13 - b13) (a21 - b21) (a22 - b22) (a23 - b23) (a31 - b31) (a32 - b32) (a33 - b33)
  M3 a11 a12 a13 a21 a22 a23 a31 a32 a33 ^^*^^ M3 b11 b12 b13 b21 b22 b23 b31 b32 b33 =
    M3 (a11*b11 + a12*b21 + a13*b31) (a11*b12 + a12*b22 + a13*b32) (a11*b13 + a12*b23 + a13*b33)
       (a21*b11 + a22*b21 + a23*b31) (a21*b12 + a22*b22 + a23*b32) (a21*b13 + a22*b23 + a23*b33)
       (a31*b11 + a32*b21 + a33*b31) (a31*b12 + a32*b22 + a33*b32) (a31*b13 + a32*b23 + a33*b33)
  a *^^ M3 b11 b12 b13 b21 b22 b23 b31 b32 b33 =
    M3 (a * b11) (a * b12) (a * b13) (a * b21) (a * b22) (a * b23) (a * b31) (a * b32) (a * b33)
  M3 a11 a12 a13 a21 a22 a23 a31 a32 a33 ^^* b =
    M3 (a11 * b) (a12 * b) (a13 * b) (a21 * b) (a22 * b) (a23 * b) (a31 * b) (a32 * b) (a33 * b)
  M3 a11 a12 a13 a21 a22 a23 a31 a32 a33 ^^/ b =
    M3 (a11 / b) (a12 / b) (a13 / b) (a21 / b) (a22 / b) (a23 / b) (a31 / b) (a32 / b) (a33 / b)
  M3 a11 a12 a13 a21 a22 a23 a31 a32 a33 `mdot` M3 b11 b12 b13 b21 b22 b23 b31 b32 b33 =
    a11*b11 + a12*b12 + a13*b13 +
    a21*b21 + a22*b22 + a23*b23 +
    a31*b31 + a32*b32 + a33*b33
  det (M3 a11 a12 a13 a21 a22 a23 a31 a32 a33) =
    let m11 = M2 a22 a23 a32 a33
        m12 = M2 a21 a23 a31 a33
        m13 = M2 a21 a22 a31 a32
    in a11 * det m11 - a12 * det m12 + a13 * det m13
  inv a@(M3 a11 a12 a13 a21 a22 a23 a31 a32 a33) =
    let m11 =   a33 * a22 - a32 * a23
        m12 = -(a33 * a12 - a32 * a13)
        m13 =   a23 * a12 - a22 * a13
        m21 = -(a33 * a21 - a31 * a23)
        m22 =   a33 * a11 - a31 * a13
        m23 = -(a23 * a11 - a21 * a13)
        m31 =   a32 * a21 - a31 * a22
        m32 = -(a32 * a11 - a31 * a12)
        m33 =   a22 * a11 - a21 * a12
    in (M3 m11 m12 m13 m21 m22 m23 m31 m32 m33) ^^/ det a

instance M M4 where
  i = M4 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
  M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44 ^^+^^ M4 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 b41 b42 b43 b44 =
    M4 (a11 + b11) (a12 + b12) (a13 + b13) (a14 + b14)
       (a21 + b21) (a22 + b22) (a23 + b23) (a24 + b24)
       (a31 + b31) (a32 + b32) (a33 + b33) (a34 + b34)
       (a41 + b41) (a42 + b42) (a43 + b43) (a44 + b44)
  M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44 ^^-^^ M4 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 b41 b42 b43 b44 =
    M4 (a11 - b11) (a12 - b12) (a13 - b13) (a14 - b14)
       (a21 - b21) (a22 - b22) (a23 - b23) (a24 - b24)
       (a31 - b31) (a32 - b32) (a33 - b33) (a34 - b34)
       (a41 - b41) (a42 - b42) (a43 - b43) (a44 - b44)
  M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44 ^^*^^ M4 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 b41 b42 b43 b44 =
    M4 (a11*b11 + a12*b21 + a13*b31 + a14*b41) (a11*b12 + a12*b22 + a13*b32 + a14*b42) (a11*b13 + a12*b23 + a13*b33 + a14*b43) (a11*b14 + a12*b24 + a13*b34 + a14*b44)
       (a21*b11 + a22*b21 + a23*b31 + a24*b41) (a21*b12 + a22*b22 + a23*b32 + a24*b42) (a21*b13 + a22*b23 + a23*b33 + a24*b43) (a21*b14 + a22*b24 + a23*b34 + a24*b44)
       (a31*b11 + a32*b21 + a33*b31 + a34*b41) (a31*b12 + a32*b22 + a33*b32 + a34*b42) (a31*b13 + a32*b23 + a33*b33 + a34*b43) (a31*b14 + a32*b24 + a33*b34 + a34*b44)
       (a41*b11 + a42*b21 + a43*b31 + a44*b41) (a41*b12 + a42*b22 + a43*b32 + a44*b42) (a41*b13 + a42*b23 + a43*b33 + a44*b43) (a41*b14 + a42*b24 + a43*b34 + a44*b44)
  a *^^ M4 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 b41 b42 b43 b44 =
    M4 (a * b11) (a * b12) (a * b13) (a * b14)
       (a * b21) (a * b22) (a * b23) (a * b24)
       (a * b31) (a * b32) (a * b33) (a * b34)
       (a * b41) (a * b42) (a * b43) (a * b44)
  M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44 ^^* b =
    M4 (a11 * b) (a12 * b) (a13 * b) (a14 * b)
       (a21 * b) (a22 * b) (a23 * b) (a24 * b)
       (a31 * b) (a32 * b) (a33 * b) (a34 * b)
       (a41 * b) (a42 * b) (a43 * b) (a44 * b)
  M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44 ^^/ b =
    M4 (a11 / b) (a12 / b) (a13 / b) (a14 / b)
       (a21 / b) (a22 / b) (a23 / b) (a24 / b)
       (a31 / b) (a32 / b) (a33 / b) (a34 / b)
       (a41 / b) (a42 / b) (a43 / b) (a44 / b)
  M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44 `mdot` M4 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 b41 b42 b43 b44 =
    a11*b11 + a12*b12 + a13*b13 + a14*b14 +
    a21*b21 + a22*b22 + a23*b23 + a24*b24 +
    a31*b31 + a32*b32 + a33*b33 + a34*b34 +
    a41*b41 + a42*b42 + a43*b43 + a44*b44
  det (M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44) =
    let m11 = M3 a22 a23 a24 a32 a33 a34 a42 a43 a44
        m12 = M3 a21 a23 a24 a31 a33 a34 a41 a43 a44
        m13 = M3 a21 a22 a24 a31 a32 a34 a41 a42 a44
        m14 = M3 a21 a22 a23 a31 a32 a33 a41 a42 a43
    in a11 * det m11 - a12 * det m12 + a13 * det m13 - a14 * det m14
  inv (M4 a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44) =
    let a = M2 a11 a12 a21 a22
        b = M2 a13 a14 a23 a24
        c = M2 a31 a32 a41 a42
        d = M2 a33 a34 a43 a44
        a1 = inv a
        ca1 = c ^^*^^ a1
        dcab1 = inv (d ^^-^^ (ca1 ^^*^^ b))
        a1bdcab1 = a1 ^^*^^ (b ^^*^^ dcab1)
        M2 m11 m12 m21 m22 = a1 ^^+^^ (a1bdcab1 ^^*^^ ca1)
        M2 m13 m14 m23 m24 = (-1) *^^ a1bdcab1
        M2 m31 m32 m41 m42 = (-1) *^^ (dcab1 ^^*^^ ca1)
        M2 m33 m34 m43 m44 = dcab1
    in M4 m11 m12 m13 m14 m21 m22 m23 m24 m31 m32 m33 m34 m41 m42 m43 m44

class MV m v where
  (^^*^) :: m -> v -> v
  (^*^!) :: v -> v -> m

reflector :: (M m, V v, MV m v) => v -> m
reflector v = i ^^-^^ (2 *^^ (v ^*^! v))

instance MV M1 V1 where
  M1 m ^^*^ V1 v = V1 (m * v)
  V1 a ^*^! V1 b = M1 (a * b)

instance MV M2 V2 where
  M2 m11 m12 m21 m22 ^^*^ V2 v1 v2 =
    V2 (m11*v1 + m12*v2)
       (m21*v1 + m22*v2)
  V2 a1 a2 ^*^! V2 b1 b2 =
    M2 (a1 * b1) (a1 * b2)
       (a2 * b1) (a2 * b2)

instance MV M3 V3 where
  M3 m11 m12 m13 m21 m22 m23 m31 m32 m33 ^^*^ V3 v1 v2 v3 =
    V3 (m11*v1 + m12*v2 + m13*v3)
       (m21*v1 + m22*v2 + m23*v3)
       (m31*v1 + m32*v2 + m33*v3)
  V3 a1 a2 a3 ^*^! V3 b1 b2 b3 =
    M3 (a1 * b1) (a1 * b2) (a1 * b3)
       (a2 * b1) (a2 * b2) (a2 * b3)
       (a3 * b1) (a3 * b2) (a3 * b3)

instance MV M4 V4 where
  M4 m11 m12 m13 m14 m21 m22 m23 m24 m31 m32 m33 m34 m41 m42 m43 m44 ^^*^ V4 v1 v2 v3 v4 =
    V4 (m11*v1 + m12*v2 + m13*v3 + m14*v4)
       (m21*v1 + m22*v2 + m23*v3 + m24*v4)
       (m31*v1 + m32*v2 + m33*v3 + m34*v4)
       (m41*v1 + m42*v2 + m43*v3 + m44*v4)
  V4 a1 a2 a3 a4 ^*^! V4 b1 b2 b3 b4 =
    M4 (a1 * b1) (a1 * b2) (a1 * b3) (a1 * b4)
       (a2 * b1) (a2 * b2) (a2 * b3) (a2 * b4)
       (a3 * b1) (a3 * b2) (a3 * b3) (a3 * b4)
       (a4 * b1) (a4 * b2) (a4 * b3) (a4 * b4)

cross4 :: V4 -> V4 -> V4 -> V4
cross4 (V4 u0 u1 u2 u3) (V4 v0 v1 v2 v3) (V4 w0 w1 w2 w3) =
{-
  let vw01 = v0 * w1 - v1 * w0
      vw02 = v0 * w2 - v2 * w0
      vw03 = v0 * w3 - v3 * w0
      vw12 = v1 * w2 - v2 * w1
      vw13 = v1 * w3 - v3 * w1
      vw23 = v2 * w3 - v3 * w2
      r0 =   u1 * vw23 - u2 * vw13 + u3 * vw12
      r1 = - u0 * vw23 + u2 * vw03 - u3 * vw02
      r2 =   u0 * vw13 - u1 * vw03 + u3 * vw02
      r3 = - u0 * vw12 + u1 * vw02 - u2 * vw01
  in V4 r0 r1 r2 r3
-}
  let m0 = M3 u1 u2 u3 v1 v2 v3 w1 w2 w3
      m1 = M3 u0 u2 u3 v0 v2 v3 w0 w2 w3
      m2 = M3 u0 u1 u3 v0 v1 v3 w0 w1 w3
      m3 = M3 u0 u1 u2 v0 v1 v2 w0 w1 w2
  in V4 (det m0) (-(det m1)) (det m2) (-(det m3))