{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
Maintainer  :  numericprelude@henning-thielemann.de
Stability   :  provisional
Portability :  portable (?)

Quaternions
-}

module Number.Quaternion
        (
        -- * Cartesian form
        T(real,imag),
        fromReal,
        (+::),

        -- * Conversions
        toRotationMatrix,
        fromRotationMatrix,
        fromRotationMatrixDenorm,
        toComplexMatrix,
        fromComplexMatrix,

        -- * Operations
        scalarProduct,
        crossProduct,
        conjugate,
        scale,
        norm,
        normSqr,
        normalize,
        similarity,
        slerp,
        )  where

import qualified Algebra.NormedSpace.Euclidean as NormedEuc
import qualified Algebra.VectorSpace  as VectorSpace
import qualified Algebra.Module       as Module
import qualified Algebra.Vector       as Vector
import qualified Algebra.Transcendental as Trans
import qualified Algebra.Algebraic    as Algebraic
import qualified Algebra.Field        as Field
import qualified Algebra.Ring         as Ring
import qualified Algebra.Additive     as Additive
import qualified Algebra.ZeroTestable as ZeroTestable

import Algebra.Module((<*>.*>), )

import qualified Number.Complex as Complex

import Number.Complex ((+:))

import qualified NumericPrelude.Elementwise as Elem
import Algebra.Additive ((<*>.+), (<*>.-), (<*>.-$), )

import Data.Array (Array, (!))
import qualified Data.Array as Array

import NumericPrelude.Base
import NumericPrelude.Numeric hiding (signum)
import Text.Show.HT (showsInfixPrec, )
import Text.Read.HT (readsInfixPrec, )


{- TODO:
conversion to and from complex matrix
-}


infix  6  +::, `Cons`

{- |
Quaternions could be defined based on Complex numbers.
However quaternions are often considered as real part and three imaginary parts.
-}
data T a
  = Cons {T a -> a
real :: !a           -- ^ real part
         ,T a -> (a, a, a)
imag :: !(a, a, a)   -- ^ imaginary parts
         }
  deriving (T a -> T a -> Bool
(T a -> T a -> Bool) -> (T a -> T a -> Bool) -> Eq (T a)
forall a. Eq a => T a -> T a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: T a -> T a -> Bool
$c/= :: forall a. Eq a => T a -> T a -> Bool
== :: T a -> T a -> Bool
$c== :: forall a. Eq a => T a -> T a -> Bool
Eq)

fromReal :: Additive.C a => a -> T a
fromReal :: a -> T a
fromReal a
x = a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons a
x (a, a, a)
forall a. C a => a
zero


plusPrec :: Int
plusPrec :: Int
plusPrec = Int
6

instance (Show a) => Show (T a) where
   showsPrec :: Int -> T a -> ShowS
showsPrec Int
prec (a
x `Cons` (a, a, a)
y) = String -> Int -> Int -> a -> (a, a, a) -> ShowS
forall a b.
(Show a, Show b) =>
String -> Int -> Int -> a -> b -> ShowS
showsInfixPrec String
"+::" Int
plusPrec Int
prec a
x (a, a, a)
y

instance (Read a) => Read (T a) where
   readsPrec :: Int -> ReadS (T a)
readsPrec Int
prec = String -> Int -> Int -> (a -> (a, a, a) -> T a) -> ReadS (T a)
forall a b c.
(Read a, Read b) =>
String -> Int -> Int -> (a -> b -> c) -> ReadS c
readsInfixPrec String
"+::" Int
plusPrec Int
prec a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
(+::)


-- | Construct a quaternion from real and imaginary part.
(+::) :: a -> (a,a,a) -> T a
+:: :: a -> (a, a, a) -> T a
(+::) = a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons

-- | The conjugate of a quaternion.
{-# SPECIALISE conjugate :: T Double -> T Double #-}
conjugate        :: (Additive.C a) => T a -> T a
conjugate :: T a -> T a
conjugate (Cons a
r (a, a, a)
i) =  a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons a
r ((a, a, a) -> (a, a, a)
forall a. C a => a -> a
negate (a, a, a)
i)

-- | Scale a quaternion by a real number.
{-# SPECIALISE scale :: Double -> T Double -> T Double #-}
scale            :: (Ring.C a) => a -> T a -> T a
scale :: a -> T a -> T a
scale a
r (Cons a
xr (a, a, a)
xi) =  a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons (a
r a -> a -> a
forall a. C a => a -> a -> a
* a
xr) (a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag a
r (a, a, a)
xi)

-- | like Module.*> but without additional class dependency
scaleImag        :: (Ring.C a) => a -> (a,a,a) -> (a,a,a)
scaleImag :: a -> (a, a, a) -> (a, a, a)
scaleImag a
r (a
xi,a
xj,a
xk) =  (a
r a -> a -> a
forall a. C a => a -> a -> a
* a
xi, a
r a -> a -> a
forall a. C a => a -> a -> a
* a
xj, a
r a -> a -> a
forall a. C a => a -> a -> a
* a
xk)

-- | the same as NormedEuc.normSqr but with a simpler type class constraint
normSqr          :: (Ring.C a) => T a -> a
normSqr :: T a -> a
normSqr (Cons a
xr (a, a, a)
xi) = a
xra -> a -> a
forall a. C a => a -> a -> a
*a
xr a -> a -> a
forall a. C a => a -> a -> a
+ (a, a, a) -> (a, a, a) -> a
forall a. C a => (a, a, a) -> (a, a, a) -> a
scalarProduct (a, a, a)
xi (a, a, a)
xi

norm             :: (Algebraic.C a) => T a -> a
norm :: T a -> a
norm T a
x = a -> a
forall a. C a => a -> a
sqrt (T a -> a
forall a. C a => T a -> a
normSqr T a
x)

-- | scale a quaternion into a unit quaternion
normalize        :: (Algebraic.C a) => T a -> T a
normalize :: T a -> T a
normalize T a
x = a -> T a -> T a
forall a. C a => a -> T a -> T a
scale (a -> a
forall a. C a => a -> a
recip (T a -> a
forall a. C a => T a -> a
norm T a
x)) T a
x

scalarProduct    :: (Ring.C a) => (a,a,a) -> (a,a,a) -> a
scalarProduct :: (a, a, a) -> (a, a, a) -> a
scalarProduct (a
xi,a
xj,a
xk) (a
yi,a
yj,a
yk) =
   a
xia -> a -> a
forall a. C a => a -> a -> a
*a
yi a -> a -> a
forall a. C a => a -> a -> a
+ a
xja -> a -> a
forall a. C a => a -> a -> a
*a
yj a -> a -> a
forall a. C a => a -> a -> a
+ a
xka -> a -> a
forall a. C a => a -> a -> a
*a
yk

crossProduct     :: (Ring.C a) => (a,a,a) -> (a,a,a) -> (a,a,a)
crossProduct :: (a, a, a) -> (a, a, a) -> (a, a, a)
crossProduct (a
xi,a
xj,a
xk) (a
yi,a
yj,a
yk) =
   (a
xja -> a -> a
forall a. C a => a -> a -> a
*a
yk a -> a -> a
forall a. C a => a -> a -> a
- a
xka -> a -> a
forall a. C a => a -> a -> a
*a
yj, a
xka -> a -> a
forall a. C a => a -> a -> a
*a
yi a -> a -> a
forall a. C a => a -> a -> a
- a
xia -> a -> a
forall a. C a => a -> a -> a
*a
yk, a
xia -> a -> a
forall a. C a => a -> a -> a
*a
yj a -> a -> a
forall a. C a => a -> a -> a
- a
xja -> a -> a
forall a. C a => a -> a -> a
*a
yi)

{- | similarity mapping as needed for rotating 3D vectors

It holds
@similarity (cos(a\/2) +:: scaleImag (sin(a\/2)) v) (0 +:: x) == (0 +:: y)@
where @y@ results from rotating @x@ around the axis @v@ by the angle @a@.
-}
similarity       :: (Field.C a) => T a -> T a -> T a
similarity :: T a -> T a -> T a
similarity T a
c T a
x = T a
cT a -> T a -> T a
forall a. C a => a -> a -> a
*T a
xT a -> T a -> T a
forall a. C a => a -> a -> a
/T a
c

{-
rotate   :: (Field.C a) =>
      (a,a,a)  {- ^ rotation axis, must be normalized -}
   -> T a
   -> T a
rotate c x = c*x/c
-}

{- |
Let @c@ be a unit quaternion, then it holds
@similarity c (0+::x) == toRotationMatrix c * x@
-}
toRotationMatrix :: (Ring.C a) => T a -> Array (Int,Int) a
toRotationMatrix :: T a -> Array (Int, Int) a
toRotationMatrix (Cons a
r (a
i,a
j,a
k)) =
   let r2 :: a
r2 = a
ra -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2
       i2 :: a
i2 = a
ia -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2;   j2 :: a
j2 = a
ja -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2;   k2 :: a
k2 = a
ka -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2
       ri :: a
ri = a
2a -> a -> a
forall a. C a => a -> a -> a
*a
ra -> a -> a
forall a. C a => a -> a -> a
*a
i; rj :: a
rj = a
2a -> a -> a
forall a. C a => a -> a -> a
*a
ra -> a -> a
forall a. C a => a -> a -> a
*a
j; rk :: a
rk = a
2a -> a -> a
forall a. C a => a -> a -> a
*a
ra -> a -> a
forall a. C a => a -> a -> a
*a
k
       jk :: a
jk = a
2a -> a -> a
forall a. C a => a -> a -> a
*a
ja -> a -> a
forall a. C a => a -> a -> a
*a
k; ki :: a
ki = a
2a -> a -> a
forall a. C a => a -> a -> a
*a
ka -> a -> a
forall a. C a => a -> a -> a
*a
i; ij :: a
ij = a
2a -> a -> a
forall a. C a => a -> a -> a
*a
ia -> a -> a
forall a. C a => a -> a -> a
*a
j
   in  ((Int, Int), (Int, Int)) -> [a] -> Array (Int, Int) a
forall i e. Ix i => (i, i) -> [e] -> Array i e
Array.listArray ((Int
0,Int
0),(Int
2,Int
2)) ([a] -> Array (Int, Int) a) -> [a] -> Array (Int, Int) a
forall a b. (a -> b) -> a -> b
$ [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[a]] -> [a]) -> [[a]] -> [a]
forall a b. (a -> b) -> a -> b
$
          [a
r2a -> a -> a
forall a. C a => a -> a -> a
+a
i2a -> a -> a
forall a. C a => a -> a -> a
-a
j2a -> a -> a
forall a. C a => a -> a -> a
-a
k2, a
ija -> a -> a
forall a. C a => a -> a -> a
-a
rk,       a
kia -> a -> a
forall a. C a => a -> a -> a
+a
rj      ] [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
:
          [a
ija -> a -> a
forall a. C a => a -> a -> a
+a
rk,       a
r2a -> a -> a
forall a. C a => a -> a -> a
-a
i2a -> a -> a
forall a. C a => a -> a -> a
+a
j2a -> a -> a
forall a. C a => a -> a -> a
-a
k2, a
jka -> a -> a
forall a. C a => a -> a -> a
-a
ri      ] [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
:
          [a
kia -> a -> a
forall a. C a => a -> a -> a
-a
rj,       a
jka -> a -> a
forall a. C a => a -> a -> a
+a
ri,       a
r2a -> a -> a
forall a. C a => a -> a -> a
-a
i2a -> a -> a
forall a. C a => a -> a -> a
-a
j2a -> a -> a
forall a. C a => a -> a -> a
+a
k2] [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
:
          []

fromRotationMatrix :: (Algebraic.C a) => Array (Int,Int) a -> T a
fromRotationMatrix :: Array (Int, Int) a -> T a
fromRotationMatrix =
   T a -> T a
forall a. C a => T a -> T a
normalize (T a -> T a)
-> (Array (Int, Int) a -> T a) -> Array (Int, Int) a -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Array (Int, Int) a -> T a
forall a. C a => Array (Int, Int) a -> T a
fromRotationMatrixDenorm


checkBounds :: (Int,Int) -> Array (Int,Int) a -> Array (Int,Int) a
checkBounds :: (Int, Int) -> Array (Int, Int) a -> Array (Int, Int) a
checkBounds (Int
c,Int
r) Array (Int, Int) a
arr =
   let bnds :: ((Int, Int), (Int, Int))
bnds@((Int
c0,Int
r0), (Int
c1,Int
r1)) = Array (Int, Int) a -> ((Int, Int), (Int, Int))
forall i e. Array i e -> (i, i)
Array.bounds Array (Int, Int) a
arr
   in  if Int
c1Int -> Int -> Int
forall a. C a => a -> a -> a
-Int
c0Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
c Bool -> Bool -> Bool
&& Int
r1Int -> Int -> Int
forall a. C a => a -> a -> a
-Int
r0Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
r
         then ((Int, Int), (Int, Int)) -> [a] -> Array (Int, Int) a
forall i e. Ix i => (i, i) -> [e] -> Array i e
Array.listArray ((Int
0,Int
0), (Int
c1Int -> Int -> Int
forall a. C a => a -> a -> a
-Int
c0, Int
r1Int -> Int -> Int
forall a. C a => a -> a -> a
-Int
r0))
                              (Array (Int, Int) a -> [a]
forall i e. Array i e -> [e]
Array.elems Array (Int, Int) a
arr)
         else String -> Array (Int, Int) a
forall a. HasCallStack => String -> a
error (String
"Quaternion.checkBounds: invalid matrix size "
                         String -> ShowS
forall a. [a] -> [a] -> [a]
++ ((Int, Int), (Int, Int)) -> String
forall a. Show a => a -> String
show ((Int, Int), (Int, Int))
bnds)


{- |
The rotation matrix must be normalized.
(I.e. no rotation with scaling)
The computed quaternion is not normalized.
-}
fromRotationMatrixDenorm :: (Ring.C a) => Array (Int,Int) a -> T a
fromRotationMatrixDenorm :: Array (Int, Int) a -> T a
fromRotationMatrixDenorm Array (Int, Int) a
mat' =
   let mat :: Array (Int, Int) a
mat = (Int, Int) -> Array (Int, Int) a -> Array (Int, Int) a
forall a. (Int, Int) -> Array (Int, Int) a -> Array (Int, Int) a
checkBounds (Int
2,Int
2) Array (Int, Int) a
mat'
       trace :: a
trace = [a] -> a
forall a. C a => [a] -> a
sum ((Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
i -> Array (Int, Int) a
mat Array (Int, Int) a -> (Int, Int) -> a
forall i e. Ix i => Array i e -> i -> e
! (Int
i,Int
i)) [Int
0..Int
2])
       dif :: (Int, Int) -> a
dif (Int
i,Int
j) = Array (Int, Int) a
matArray (Int, Int) a -> (Int, Int) -> a
forall i e. Ix i => Array i e -> i -> e
!(Int
i,Int
j) a -> a -> a
forall a. C a => a -> a -> a
- Array (Int, Int) a
matArray (Int, Int) a -> (Int, Int) -> a
forall i e. Ix i => Array i e -> i -> e
!(Int
j,Int
i)
   in  a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons (a
tracea -> a -> a
forall a. C a => a -> a -> a
+a
1) ((Int, Int) -> a
dif (Int
2,Int
1), (Int, Int) -> a
dif (Int
0,Int
2), (Int, Int) -> a
dif (Int
1,Int
0))

{- |
Map a quaternion to complex valued 2x2 matrix,
such that quaternion addition and multiplication
is mapped to matrix addition and multiplication.
The determinant of the matrix equals the squared quaternion norm ('normSqr').
Since complex numbers can be turned into real (orthogonal) matrices,
a quaternion could also be converted into a real matrix.
-}
toComplexMatrix :: (Additive.C a) =>
   T a -> Array (Int,Int) (Complex.T a)
toComplexMatrix :: T a -> Array (Int, Int) (T a)
toComplexMatrix (Cons a
r (a
i,a
j,a
k)) =
   ((Int, Int), (Int, Int)) -> [T a] -> Array (Int, Int) (T a)
forall i e. Ix i => (i, i) -> [e] -> Array i e
Array.listArray ((Int
0,Int
0), (Int
1,Int
1))
      [a
ra -> a -> T a
forall a. a -> a -> T a
+:a
i, (-a
j)a -> a -> T a
forall a. a -> a -> T a
+:(-a
k), a
ja -> a -> T a
forall a. a -> a -> T a
+:(-a
k), a
ra -> a -> T a
forall a. a -> a -> T a
+:(-a
i)]


{- |
Revert 'toComplexMatrix'.
-}
fromComplexMatrix :: (Field.C a) =>
   Array (Int,Int) (Complex.T a) -> T a
fromComplexMatrix :: Array (Int, Int) (T a) -> T a
fromComplexMatrix Array (Int, Int) (T a)
mat =
   let xs :: [T a]
xs = Array (Int, Int) (T a) -> [T a]
forall i e. Array i e -> [e]
Array.elems ((Int, Int) -> Array (Int, Int) (T a) -> Array (Int, Int) (T a)
forall a. (Int, Int) -> Array (Int, Int) a -> Array (Int, Int) a
checkBounds (Int
1,Int
1) Array (Int, Int) (T a)
mat)
       [a
ar,a
br,a
cr,a
dr] = (T a -> a) -> [T a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map T a -> a
forall a. T a -> a
Complex.real [T a]
xs
       [a
ai,a
bi,a
ci,a
di] = (T a -> a) -> [T a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map T a -> a
forall a. T a -> a
Complex.imag [T a]
xs
   in  a -> T a -> T a
forall a. C a => a -> T a -> T a
scale (a
1a -> a -> a
forall a. C a => a -> a -> a
/a
2) (a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons (a
ara -> a -> a
forall a. C a => a -> a -> a
+a
dr) (a
aia -> a -> a
forall a. C a => a -> a -> a
-a
di, a
cra -> a -> a
forall a. C a => a -> a -> a
-a
br, -a
cia -> a -> a
forall a. C a => a -> a -> a
-a
bi))


{- |
Spherical Linear Interpolation

Can be generalized to any transcendent Hilbert space.
In fact, we should also include the real part in the interpolation.
-}
slerp :: (Trans.C a) =>
      a   {- ^ For @0@ return vector @v@,
               for @1@ return vector @w@ -}
   -> (a,a,a)  {- ^ vector @v@, must be normalized -}
   -> (a,a,a)  {- ^ vector @w@, must be normalized -}
   -> (a,a,a)
slerp :: a -> (a, a, a) -> (a, a, a) -> (a, a, a)
slerp a
c (a, a, a)
v (a, a, a)
w =
   let scal :: a
scal  = (a, a, a) -> (a, a, a) -> a
forall a. C a => (a, a, a) -> (a, a, a) -> a
scalarProduct (a, a, a)
v (a, a, a)
w a -> a -> a
forall a. C a => a -> a -> a
/
                  a -> a
forall a. C a => a -> a
sqrt ((a, a, a) -> (a, a, a) -> a
forall a. C a => (a, a, a) -> (a, a, a) -> a
scalarProduct (a, a, a)
v (a, a, a)
v a -> a -> a
forall a. C a => a -> a -> a
* (a, a, a) -> (a, a, a) -> a
forall a. C a => (a, a, a) -> (a, a, a) -> a
scalarProduct (a, a, a)
w (a, a, a)
w)
       angle :: a
angle = a -> a
forall a. C a => a -> a
Trans.acos a
scal
   in  a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag (a -> a
forall a. C a => a -> a
recip (a -> a
forall a. C a => a -> a
Algebraic.sqrt (a
1a -> a -> a
forall a. C a => a -> a -> a
-a
scala -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2)))
         (a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag (a -> a
forall a. C a => a -> a
Trans.sin ((a
1a -> a -> a
forall a. C a => a -> a -> a
-a
c)a -> a -> a
forall a. C a => a -> a -> a
*a
angle)) (a, a, a)
v (a, a, a) -> (a, a, a) -> (a, a, a)
forall a. C a => a -> a -> a
+
          a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag (a -> a
forall a. C a => a -> a
Trans.sin (   a
c a -> a -> a
forall a. C a => a -> a -> a
*a
angle)) (a, a, a)
w)



instance (NormedEuc.Sqr a b) => NormedEuc.Sqr a (T b) where
   normSqr :: T b -> a
normSqr (Cons b
r (b, b, b)
i) = b -> a
forall a v. Sqr a v => v -> a
NormedEuc.normSqr b
r a -> a -> a
forall a. C a => a -> a -> a
+ (b, b, b) -> a
forall a v. Sqr a v => v -> a
NormedEuc.normSqr (b, b, b)
i

instance (Algebraic.C a, NormedEuc.Sqr a b) => NormedEuc.C a (T b) where
   norm :: T b -> a
norm = T b -> a
forall a v. (C a, Sqr a v) => v -> a
NormedEuc.defltNorm



instance (ZeroTestable.C a) => ZeroTestable.C (T a)  where
   isZero :: T a -> Bool
isZero (Cons a
r (a, a, a)
i)  = a -> Bool
forall a. C a => a -> Bool
isZero a
r Bool -> Bool -> Bool
&& (a, a, a) -> Bool
forall a. C a => a -> Bool
isZero (a, a, a)
i

instance (Additive.C a) => Additive.C (T a)  where
   {-# SPECIALISE instance Additive.C (T Float) #-}
   {-# SPECIALISE instance Additive.C (T Double) #-}
   zero :: T a
zero   = a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons a
forall a. C a => a
zero (a, a, a)
forall a. C a => a
zero
   + :: T a -> T a -> T a
(+)    = T (T a, T a) (T a) -> T a -> T a -> T a
forall x y a. T (x, y) a -> x -> y -> a
Elem.run2 (T (T a, T a) (T a) -> T a -> T a -> T a)
-> T (T a, T a) (T a) -> T a -> T a -> T a
forall a b. (a -> b) -> a -> b
$ (a -> (a, a, a) -> T a) -> T (T a, T a) (a -> (a, a, a) -> T a)
forall a v. a -> T v a
Elem.with a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons T (T a, T a) (a -> (a, a, a) -> T a)
-> (T a -> a) -> T (T a, T a) ((a, a, a) -> T a)
forall x v a. C x => T (v, v) (x -> a) -> (v -> x) -> T (v, v) a
<*>.+  T a -> a
forall a. T a -> a
real T (T a, T a) ((a, a, a) -> T a)
-> (T a -> (a, a, a)) -> T (T a, T a) (T a)
forall x v a. C x => T (v, v) (x -> a) -> (v -> x) -> T (v, v) a
<*>.+  T a -> (a, a, a)
forall a. T a -> (a, a, a)
imag
   (-)    = T (T a, T a) (T a) -> T a -> T a -> T a
forall x y a. T (x, y) a -> x -> y -> a
Elem.run2 (T (T a, T a) (T a) -> T a -> T a -> T a)
-> T (T a, T a) (T a) -> T a -> T a -> T a
forall a b. (a -> b) -> a -> b
$ (a -> (a, a, a) -> T a) -> T (T a, T a) (a -> (a, a, a) -> T a)
forall a v. a -> T v a
Elem.with a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons T (T a, T a) (a -> (a, a, a) -> T a)
-> (T a -> a) -> T (T a, T a) ((a, a, a) -> T a)
forall x v a. C x => T (v, v) (x -> a) -> (v -> x) -> T (v, v) a
<*>.-  T a -> a
forall a. T a -> a
real T (T a, T a) ((a, a, a) -> T a)
-> (T a -> (a, a, a)) -> T (T a, T a) (T a)
forall x v a. C x => T (v, v) (x -> a) -> (v -> x) -> T (v, v) a
<*>.-  T a -> (a, a, a)
forall a. T a -> (a, a, a)
imag
   negate :: T a -> T a
negate = T (T a) (T a) -> T a -> T a
forall v a. T v a -> v -> a
Elem.run  (T (T a) (T a) -> T a -> T a) -> T (T a) (T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$ (a -> (a, a, a) -> T a) -> T (T a) (a -> (a, a, a) -> T a)
forall a v. a -> T v a
Elem.with a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons T (T a) (a -> (a, a, a) -> T a)
-> (T a -> a) -> T (T a) ((a, a, a) -> T a)
forall x v a. C x => T v (x -> a) -> (v -> x) -> T v a
<*>.-$ T a -> a
forall a. T a -> a
real T (T a) ((a, a, a) -> T a) -> (T a -> (a, a, a)) -> T (T a) (T a)
forall x v a. C x => T v (x -> a) -> (v -> x) -> T v a
<*>.-$ T a -> (a, a, a)
forall a. T a -> (a, a, a)
imag

instance (Ring.C a) => Ring.C (T a)  where
   {-# SPECIALISE instance Ring.C (T Float) #-}
   {-# SPECIALISE instance Ring.C (T Double) #-}
   one :: T a
one                          =  a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons a
forall a. C a => a
one (a, a, a)
forall a. C a => a
zero
   fromInteger :: Integer -> T a
fromInteger                  =  a -> T a
forall a. C a => a -> T a
fromReal (a -> T a) -> (Integer -> a) -> Integer -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> a
forall a. C a => Integer -> a
fromInteger
   (Cons a
xr (a, a, a)
xi) * :: T a -> T a -> T a
* (Cons a
yr (a, a, a)
yi)  =
       a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons (a
xra -> a -> a
forall a. C a => a -> a -> a
*a
yr a -> a -> a
forall a. C a => a -> a -> a
- (a, a, a) -> (a, a, a) -> a
forall a. C a => (a, a, a) -> (a, a, a) -> a
scalarProduct (a, a, a)
xi (a, a, a)
yi)
            (a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag a
xr (a, a, a)
yi (a, a, a) -> (a, a, a) -> (a, a, a)
forall a. C a => a -> a -> a
+ a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag a
yr (a, a, a)
xi (a, a, a) -> (a, a, a) -> (a, a, a)
forall a. C a => a -> a -> a
+
             (a, a, a) -> (a, a, a) -> (a, a, a)
forall a. C a => (a, a, a) -> (a, a, a) -> (a, a, a)
crossProduct (a, a, a)
xi (a, a, a)
yi)

instance (Field.C a) => Field.C (T a)  where
   {-# SPECIALISE instance Field.C (T Float) #-}
   {-# SPECIALISE instance Field.C (T Double) #-}
   recip :: T a -> T a
recip T a
x = a -> T a -> T a
forall a. C a => a -> T a -> T a
scale (a -> a
forall a. C a => a -> a
recip (T a -> a
forall a. C a => T a -> a
normSqr T a
x)) (T a -> T a
forall a. C a => T a -> T a
conjugate T a
x)
   (Cons a
xr (a, a, a)
xi) / :: T a -> T a -> T a
/ y :: T a
y@(Cons a
yr (a, a, a)
yi) =
       a -> T a -> T a
forall a. C a => a -> T a -> T a
scale (a -> a
forall a. C a => a -> a
recip (T a -> a
forall a. C a => T a -> a
normSqr T a
y))
          (a -> (a, a, a) -> T a
forall a. a -> (a, a, a) -> T a
Cons (a
xra -> a -> a
forall a. C a => a -> a -> a
*a
yr a -> a -> a
forall a. C a => a -> a -> a
+ (a, a, a) -> (a, a, a) -> a
forall a. C a => (a, a, a) -> (a, a, a) -> a
scalarProduct (a, a, a)
xi (a, a, a)
yi)
                (a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag a
yr (a, a, a)
xi (a, a, a) -> (a, a, a) -> (a, a, a)
forall a. C a => a -> a -> a
- a -> (a, a, a) -> (a, a, a)
forall a. C a => a -> (a, a, a) -> (a, a, a)
scaleImag a
xr (a, a, a)
yi (a, a, a) -> (a, a, a) -> (a, a, a)
forall a. C a => a -> a -> a
- (a, a, a) -> (a, a, a) -> (a, a, a)
forall a. C a => (a, a, a) -> (a, a, a) -> (a, a, a)
crossProduct (a, a, a)
xi (a, a, a)
yi))

instance Vector.C T where
   zero :: T a
zero  = T a
forall a. C a => a
zero
   <+> :: T a -> T a -> T a
(<+>) = T a -> T a -> T a
forall a. C a => a -> a -> a
(+)
   *> :: a -> T a -> T a
(*>)  = a -> T a -> T a
forall a. C a => a -> T a -> T a
scale

-- | The '(*>)' method can't replace 'scale'
--   because it requires the Algebra.Module constraint
instance (Module.C a b) => Module.C a (T b) where
   *> :: a -> T b -> T b
(*>) = T (a, T b) (T b) -> a -> T b -> T b
forall x y a. T (x, y) a -> x -> y -> a
Elem.run2 (T (a, T b) (T b) -> a -> T b -> T b)
-> T (a, T b) (T b) -> a -> T b -> T b
forall a b. (a -> b) -> a -> b
$ (b -> (b, b, b) -> T b) -> T (a, T b) (b -> (b, b, b) -> T b)
forall a v. a -> T v a
Elem.with b -> (b, b, b) -> T b
forall a. a -> (a, a, a) -> T a
Cons T (a, T b) (b -> (b, b, b) -> T b)
-> (T b -> b) -> T (a, T b) ((b, b, b) -> T b)
forall a x v c.
C a x =>
T (a, v) (x -> c) -> (v -> x) -> T (a, v) c
<*>.*> T b -> b
forall a. T a -> a
real T (a, T b) ((b, b, b) -> T b)
-> (T b -> (b, b, b)) -> T (a, T b) (T b)
forall a x v c.
C a x =>
T (a, v) (x -> c) -> (v -> x) -> T (a, v) c
<*>.*> T b -> (b, b, b)
forall a. T a -> (a, a, a)
imag

instance (VectorSpace.C a b) => VectorSpace.C a (T b)