{-# 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))