import Math.CayleyDickson
import Data.Proxy (Proxy(Proxy))
import Data.Ratio (Ratio, (%))
import Control.Monad (replicateM, replicateM_, forM_, liftM)
import System.Random (Random, randomRIO)
----------------------------------------------------------
-- alternate formulas
pureDir :: (Tag n, Conjugable a, Floating a) => Nion n a -> Nion n a
pureDir x = p /. (norm p) where p = purePart x
cos' :: (Tag n, Conjugable a, RealFloat a) => Nion n a -> Nion n a
cos' x = (exp (u * x) + exp (- u * x)) / 2 where u = pureDir x
sin' :: (Tag n, Conjugable a, RealFloat a) => Nion n a -> Nion n a
sin' x = ((exp (u * x) - exp (- u * x)) * recip u) / 2 where u = pureDir x
cosh' :: (Tag n, Conjugable a, RealFloat a) => Nion n a -> Nion n a
cosh' x = (exp x + exp (- x)) / 2
sinh' :: (Tag n, Conjugable a, RealFloat a) => Nion n a -> Nion n a
sinh' x = (exp x - exp (- x)) / 2
dot' :: (Tag n, Conjugable a, Fractional a) => Nion n a -> Nion n a -> a
dot' x y = scalarPart $ (y * conj x + x * conj y) / 2
cross' :: (Tag n, Conjugable a, Fractional a) =>
Nion n a -> Nion n a -> Nion n a
cross' x y = (y * conj x - x * conj y) / 2
qmul :: Num a => a -> a -> a -> a ->
a -> a -> a -> a -> (a, a, a, a)
qmul aw ax ay az bw bx by bz = (aw*bw - ax*bx - ay*by - az*bz,
aw*bx + ax*bw + ay*bz - az*by,
aw*by - ax*bz + ay*bw + az*bx,
aw*bz + ax*by - ay*bx + az*bw)
----------------------------------------------------------
-- test utils
epsilon :: Double
epsilon = 1e-7
assert :: Bool -> IO ()
assert True = putChar '.'
assert False = error "assertion failed"
close :: Tag n => Nion n Double -> Nion n Double -> Bool
close x y = norm (x - y) < epsilon
closeReal :: Double -> Double -> Bool
closeReal x y = abs (x - y) < epsilon
assertClose :: Tag n => Nion n Double -> Nion n Double -> IO ()
assertClose x y = assert $ close x y
assertCloseReal :: Double -> Double -> IO ()
assertCloseReal x y = assert $ closeReal x y
boundsI :: (Integer, Integer)
boundsI = (-100000, 100000)
boundsD :: (Double, Double)
boundsD = (-1, 1)
randomElt' :: (Tag n, Conjugable a, Random a) =>
Proxy n -> (a, a) -> IO (Nion n a)
randomElt' n' bounds = liftM nion $ replicateM (2^n) (randomRIO bounds) where
n = tagVal n'
randomElt :: (Tag n, Conjugable a, Random a) => (a, a) -> IO (Nion n a)
randomElt = randomElt' Proxy
randomEltD :: Tag n => IO (Nion n Double)
randomEltD = randomElt boundsD
randomEltI :: Tag n => IO (Nion n Integer)
randomEltI = randomElt boundsI
randomEltI' :: (Tag n1, Tag n2) => Integer -> IO (Nion n1 (Nion n2 Integer))
randomEltI' n = liftM nion $ replicateM (fromIntegral n) randomEltI
randomEltI2 :: Tag n => IO (Complex (Nion n Integer))
randomEltI2 = randomEltI' 2
randomEltI4 :: Tag n => IO (Quaternion (Nion n Integer))
randomEltI4 = randomEltI' 4
----------------------------------------------------------
-- checks
checkFloating1' :: Tag n => Proxy n -> Nion n Double -> IO ()
checkFloating1' n' x = do
if sqnorm (purePart x) /= 0
then do assertClose (cos x) (cos' x)
assertClose (sin x) (sin' x)
else return ()
assertClose (cosh x) (cosh' x)
assertClose (sinh x) (sinh' x)
if sqnorm x /= 0
then assertCloseReal 1 (scalarPart $ x * recip x)
else return ()
forM_ (zip (coords x) ([0..] :: [Integer])) $ \(e, i) -> do
assert $ e == coord x i
assert $ 999 == coord (setCoord x i 999) i
if n /= 0
then do let (s, t, u) = polar x
assertClose x $ s .* exp (t .* u)
else return ()
where
n = tagVal n'
checkFloating1 :: Tag n => Nion n Double -> IO ()
checkFloating1 = checkFloating1' Proxy
checkFloating2' :: Tag n => Proxy n -> Nion n Double -> Nion n Double -> IO ()
checkFloating2' _ x y = do
assertCloseReal (x `dot` y) (x `dot'` y)
assertCloseReal (5 `dot` y) (5 `dot'` y)
assertCloseReal (x `dot` 5) (x `dot'` 5)
assertClose (x `cross` y) (x `cross'` y)
assertClose (5 `cross` y) (5 `cross'` y)
assertClose (x `cross` 5) (x `cross'` 5)
checkFloating2 :: Tag n => Nion n Double -> Nion n Double -> IO ()
checkFloating2 = checkFloating2' Proxy
checkFloating3 :: Tag n => Nion n Double -> IO ()
checkFloating3 x' = do
assertCloseReal (scalarPart $ exp x) (exp $ scalarPart x)
assertClose (purePart $ exp x) 0
assertCloseReal (scalarPart $ cos x) (cos $ scalarPart x)
assertClose (purePart $ cos x) 0
assertCloseReal (scalarPart $ sin x) (sin $ scalarPart x)
assertClose (purePart $ sin x) 0
assertCloseReal (scalarPart $ cosh x) (cosh $ scalarPart x)
assertClose (purePart $ cosh x) 0
assertCloseReal (scalarPart $ sinh x) (sinh $ scalarPart x)
assertClose (purePart $ sinh x) 0
where
x = scalarPart x' .+ (x' - x')
checkFloating' :: Tag n => IO (Nion n Double) -> IO (Nion n Double) -> IO ()
checkFloating' x y = do
x' <- x
y' <- y
checkFloating1 x'
checkFloating2 x' y'
checkFloating3 x'
checkFloating :: IO ()
checkFloating = do
checkFloating' (randomEltD :: IO (Nion Tag0 Double))
(randomEltD :: IO (Nion Tag0 Double))
checkFloating' (randomEltD :: IO (Complex Double))
(randomEltD :: IO (Complex Double))
checkFloating' (randomEltD :: IO (Quaternion Double))
(randomEltD :: IO (Quaternion Double))
checkFloating' (randomEltD :: IO (Sedenion Double))
(randomEltD :: IO (Sedenion Double))
checkScalar :: IO ()
checkScalar = do
let x = nion [3] :: Nion Tag0 Integer
y = nion [4] :: Nion Tag0 Integer
z = nion [12] :: Nion Tag0 Integer
assert $ x * y == z
checkComplex :: IO ()
checkComplex = do
let x = complex 1 0 :: Complex Integer
y = complex 0 1 :: Complex Integer
z = complex 1 1 :: Complex Integer
assert $ x + y == z
assert $ x == basisElement (0 :: Integer)
assert $ y == basisElement (1 :: Integer)
assert $ (complex 1 2 :: Complex Integer) == nion [1..]
checkQuaternion :: IO ()
checkQuaternion = do
let x = quaternion 0 1 0 0 :: Quaternion Integer
y = quaternion 0 0 1 0 :: Quaternion Integer
z = quaternion 0 0 0 1 :: Quaternion Integer
assert $ x == basisElement (1 :: Integer)
assert $ y == basisElement (2 :: Integer)
assert $ z == basisElement (3 :: Integer)
assert $ x `cross` y == z
assert $ (quaternion 1 2 3 4 :: Quaternion Integer) == nion [1..]
a <- randomEltI :: IO (Quaternion Integer)
b <- randomEltI :: IO (Quaternion Integer)
let (a0:a1:a2:a3:[]) = coords a
(b0:b1:b2:b3:[]) = coords b
(cw, cx, cy, cz) = qmul a0 a1 a2 a3 b0 b1 b2 b3
assert $ quaternion cw cx cy cz == a * b
checkOctonion :: IO ()
checkOctonion = do
let x = octonion 0 1 0 0 0 0 0 0 :: Octonion Integer
y = octonion 0 0 1 0 0 0 0 0 :: Octonion Integer
z = octonion 0 0 0 1 0 0 0 0 :: Octonion Integer
assert $ x == basisElement (1 :: Integer)
assert $ y == basisElement (2 :: Integer)
assert $ z == basisElement (3 :: Integer)
assert $ x `cross` y == z
assert $ (octonion 1 2 3 4 5 6 7 8 :: Octonion Integer) == nion [1..]
checkSedenion :: IO ()
checkSedenion = do
let x = sedenion 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 :: Sedenion Integer
y = sedenion 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 :: Sedenion Integer
z = sedenion 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 :: Sedenion Integer
e = basisElement :: Integer -> Sedenion Integer
assert $ x == e 1
assert $ y == e 2
assert $ z == e 3
assert $ x `cross` y == z
assert $ (e 3 + e 10) * (e 6 - e 15) == 0
assert $ (sedenion 1 2 3 4 5 6 7 8
9 10 11 12 13 14 15 16 :: Sedenion Integer) == nion [1..]
checkBig :: IO ()
checkBig = do
let x = nion [0,1,0,0] :: Nion Tag6 Integer
y = nion [0,0,1,0] :: Nion Tag6 Integer
z = nion [0,0,0,1] :: Nion Tag6 Integer
assert $ x `cross` y == z
checkMixed :: IO ()
checkMixed = do
let x = quaternion 1 2 3 4 :: Quaternion Integer
assert $ 10 .+ x == quaternion 11 2 3 4
assert $ x +. 10 == quaternion 11 2 3 4
assert $ 10 .- x == nion [9, -2, -3, -4]
assert $ x -. 10 == nion [-9, 2, 3, 4]
assert $ 10 .* x == quaternion 10 20 30 40
assert $ x *. 10 == quaternion 10 20 30 40
assert $ (quaternion 1 2 3 4 :: Quaternion (Ratio Integer)) /. 2 ==
nion [1 % 2, 1, 3 % 2, 2]
forM_ [-2, -1, -0.5, 0, 0.5, 1, 2] $ \r -> do
let y = fromScalar r :: Quaternion Double
(s, t, u) = polar y
assert $ s >= 0
assertClose (fromScalar r) $ s .* exp (t .* u)
assertClose (sqrt y) $ (sqrt s) .* exp ((t .* u) /. 2)
checkInverses :: IO ()
checkInverses = do
f $ (nion [0.1] :: Nion Tag0 Double)
f $ complex 0.1 0.2
f $ quaternion 0.1 0.2 0.3 0.4
f $ octonion 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
return ()
where
f :: Tag n => Nion n Double -> IO ()
f x = do
assertClose x $ (cos . acos) x
assertClose x $ (acos . cos) x
assertClose x $ (sin . asin) x
assertClose x $ (asin . sin) x
assertClose x $ (tan . atan) x
assertClose x $ (atan . tan) x
checkBasic :: IO ()
checkBasic = do
let x = quaternion 3 4 5 6 :: Quaternion Integer
assert $ negate x == quaternion (-3) (-4) (-5) (-6)
assert $ sqnorm x == 3^(2::Integer) + 4^(2::Integer) +
5^(2::Integer) + 6^(2::Integer)
assert $ x + 99 == quaternion 102 4 5 6
assert $ 99 + x == quaternion 102 4 5 6
assert $ x - 1 == quaternion 2 4 5 6
assert $ 1 - x == quaternion (-2) (-4) (-5) (-6)
assert $ x * 2 == quaternion 6 8 10 12
assert $ 2 * x == quaternion 6 8 10 12
assert $ x `dot` 7 == 21
assert $ 7 `dot` x == 21
assert $ coord (nion [5] :: Nion Tag0 Integer) (0::Integer) == 5
assert $ coord (nion [5] :: Nion Tag4 Integer) (0::Integer) == 5
assert $ coord (nion [5] :: Nion Tag4 Integer) (1::Integer) == 0
assert $ setCoord (nion [5] :: Nion Tag0 Integer)
(0::Integer) (9::Integer) == nion [9]
assert $ setCoord (nion [5] :: Nion Tag4 Integer) (1::Integer) (9::Integer) ==
nion [5, 9]
assert $ (fromScalar 5 :: Nion Tag0 Integer) == nion [5]
assert $ (fromScalar 5 :: Nion Tag4 Integer) == nion [5]
assert $ nion [5] == (fromScalar 5 :: Nion Tag4 Integer)
let y = quaternion 1 2 3 4 :: Quaternion (Ratio Integer)
assert $ y / 2 == quaternion (1 % 2) 1 (3 % 2) 2
assert $ recip y == quaternion (1 % 30) (-1 % 15) (-1 % 10) (-2 % 15)
assert $ y * recip y == 1
checkPower :: IO ()
checkPower = do
let x = quaternion 1 2 3 4 :: Quaternion Integer
assert $ x ^. (0 :: Integer) == 1
assert $ x ^. (1 :: Integer) == x
assert $ x ^. (2 :: Integer) == x * x
assert $ x ^. (3 :: Integer) == x * x * x
assert $ x ^. (4 :: Integer) == x * x * x * x
let y = quaternion 1 2 3 4 :: Quaternion (Ratio Integer)
assert $ y ^^. (0 :: Integer) == 1
assert $ y ^^. (1 :: Integer) == y
assert $ y ^^. (2 :: Integer) == y * y
assert $ y ^^. (3 :: Integer) == y * y * y
assert $ y ^^. (4 :: Integer) == y * y * y * y
assert $ y ^^. (-1 :: Integer) == recip y
assert $ y ^^. (-2 :: Integer) == recip (y * y)
checkZeroAndOne :: (Show a, Eq a, Conjugable a) => Nion n1 (Nion n2 a) -> IO ()
checkZeroAndOne x = do
assert $ 0 + x == x
assert $ 1 * x == x
checkDistributive :: (Show a, Eq a, Conjugable a) =>
Nion n1 (Nion n2 a) -> Nion n1 (Nion n2 a) ->
Nion n2 a -> Nion n2 a ->
IO ()
checkDistributive x y r s = do
assert $ (r + s) .* x == r .* x + s .* x
assert $ r .* (x + y) == r .* x + r .* y
checkModule :: (Show a, Eq a, Conjugable a) =>
Nion n1 (Nion n2 a) -> Nion n1 (Nion n2 a) ->
Nion n2 a -> Nion n2 a -> IO ()
checkModule x y r s = do
checkDistributive x y r s
checkZeroAndOne x
assert $ (r * s) .* x == r .* (s .* x)
checkIsomorphism :: (Conjugable a, Show a, Eq a) =>
((Nion n1 (Nion n2 a)) -> Nion n3 a) ->
(Nion n1 (Nion n2 a)) -> (Nion n1 (Nion n2 a)) ->
IO ()
checkIsomorphism f x y = do
assert $ f 0 == 0
assert $ f 1 == 1
assert $ f (conj x) == conj (f x)
assert $ f (negate x) == negate (f x)
assert $ f (x + y) == f x + f y
assert $ f (x - y) == f x - f y
assert $ f (x * y) == f x * f y
assert $ scalarPart (sqnorm x) == sqnorm (f x)
phi :: (Tag n1, Tag n2, Tag n3, Conjugable a) =>
(Nion n1 (Nion n2 a)) -> Nion n3 a
phi = nion . concatMap coords . coords
checkProperties1 :: IO ()
checkProperties1 = do
let f = phi :: Complex (Complex Integer) -> Quaternion Integer
r <- randomEltI :: IO (Complex Integer)
s <- randomEltI :: IO (Complex Integer)
x <- randomEltI2 :: IO (Complex (Complex Integer))
y <- randomEltI2 :: IO (Complex (Complex Integer))
checkIsomorphism f x y
checkModule x y r s
checkProperties2 :: IO ()
checkProperties2 = do
let f = phi :: Complex (Quaternion Integer) -> Octonion Integer
r <- randomEltI :: IO (Quaternion Integer)
s <- randomEltI :: IO (Quaternion Integer)
x <- randomEltI2 :: IO (Complex (Quaternion Integer))
y <- randomEltI2 :: IO (Complex (Quaternion Integer))
checkIsomorphism f x y
checkModule x y r s
checkProperties3 :: IO ()
checkProperties3 = do
let f = phi :: Complex (Octonion Integer) -> Sedenion Integer
r <- randomEltI :: IO (Octonion Integer)
s <- randomEltI :: IO (Octonion Integer)
x <- randomEltI2 :: IO (Complex (Octonion Integer))
y <- randomEltI2 :: IO (Complex (Octonion Integer))
checkIsomorphism f x y
checkDistributive x y r s
checkZeroAndOne x
checkProperties4 :: IO ()
checkProperties4 = do
let f = phi :: Quaternion (Complex Integer) -> Octonion Integer
r <- randomEltI :: IO (Complex Integer)
s <- randomEltI :: IO (Complex Integer)
x <- randomEltI4 :: IO (Quaternion (Complex Integer))
y <- randomEltI4 :: IO (Quaternion (Complex Integer))
checkIsomorphism f x y
checkModule x y r s
main :: IO ()
main = do
checkBasic
checkScalar
checkComplex
checkQuaternion
checkOctonion
checkSedenion
checkBig
checkPower
checkInverses
checkMixed
replicateM_ 20 $ do
checkFloating
checkProperties1
checkProperties2
checkProperties3
checkProperties4
putStrLn "\nAll tests passed."