diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,20 @@ +Copyright (c) 2015 James M. Lawrence + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be included +in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,7 @@ +Cayley-Dickson constructions (complex numbers, quaternions, octonions, +sedenions, etc.) over general scalars without limit to the number of +dimensions. + +License: MIT + +Author: James M. Lawrence diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cayley-dickson.cabal b/cayley-dickson.cabal new file mode 100644 --- /dev/null +++ b/cayley-dickson.cabal @@ -0,0 +1,34 @@ +name: cayley-dickson +version: 0.1.0.0 +synopsis: Complex numbers, quaternions, octonions, sedenions, etc. +description: Cayley-Dickson constructions (complex numbers, quaternions, + octonions, sedenions, etc.) over general scalars without + limit to the number of dimensions. +license: MIT +license-file: LICENSE +author: James M. Lawrence +maintainer: James M. Lawrence +copyright: Copyright (c) James M. Lawrence +category: Algebra, Data, Data Structures, Math +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.10 + +source-repository head + type: git + location: git://github.com/lmj/cayley-dickson.git + +library + exposed-modules: Math.CayleyDickson + hs-source-dirs: src + build-depends: base >= 4.8 && < 5 + default-language: Haskell2010 + ghc-options: -Wall -O2 + +Test-Suite test + type: exitcode-stdio-1.0 + hs-source-dirs: test, src + main-is: test.hs + build-depends: base >= 4.8 && < 5, random >= 1 + default-language: Haskell2010 + ghc-options: -Wall diff --git a/src/Math/CayleyDickson.hs b/src/Math/CayleyDickson.hs new file mode 100644 --- /dev/null +++ b/src/Math/CayleyDickson.hs @@ -0,0 +1,613 @@ +-- Copyright (c) 2015 James M. Lawrence +-- +-- Permission is hereby granted, free of charge, to any person obtaining +-- a copy of this software and associated documentation files (the +-- "Software"), to deal in the Software without restriction, including +-- without limitation the rights to use, copy, modify, merge, publish, +-- distribute, sublicense, and/or sell copies of the Software, and to +-- permit persons to whom the Software is furnished to do so, subject to +-- the following conditions: +-- +-- The above copyright notice and this permission notice shall be included +-- in all copies or substantial portions of the Software. +-- +-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +-- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +-- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +-- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +-- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +-- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +-- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +-- +----------------------------------------------------------------------------- +-- | +-- Module : Math.CayleyDickson +-- Copyright : (c) James M. Lawrence +-- License : MIT +-- +-- Maintainer : James M. Lawrence +-- Stability : provisional +-- Portability : portable +-- +-- Cayley-Dickson constructions (complex numbers, quaternions, +-- octonions, sedenions, etc.) over general scalars without limit to +-- the number of dimensions. +-- +-- An element of this structure is composed of an m-dimensional +-- /scalar part/ and an m*(2^n - 1)-dimensional /pure part/ (unrelated +-- to Haskell's uses of "pure"). An element whose scalar part is zero +-- is called a /pure/. Construction with real scalars yields the +-- Cayley-Dickson algebras, in which case the scalar part is also +-- called the /real part/. Other structures may be obtained by +-- considering general scalars, for instance the quaternions over +-- complex scalars. +----------------------------------------------------------------------------- + +module Math.CayleyDickson ( + -- * Types + Nion, Complex, Quaternion, Octonion, Sedenion, + + -- * Construction + nion, fromScalar, complex, quaternion, octonion, sedenion, + + -- * Operations + dot, cross, sqnorm, norm, polar, + -- ** Operations with scalars + -- + -- | The mnemonic is that the period (".") is on the side of the + -- scalar. + (^.), (^^.), (**.), + (.+), (+.), (.-), (-.), (.*), (*.), (/.), + + -- * Accessors + coord, coords, setCoord, scalarPart, purePart, + + -- * Constants + basisElement, + + -- * Classes + Conjugable(conj), + + -- ** Tags + -- + -- | Tags serve to determine a type's dimension, which is 2 raised + -- to `tagVal`. Tag instances are included for convenience only, + -- as you may create your own tag. + Tag(tagVal), + Tag0, Tag1, Tag2, Tag3, Tag4, Tag5, Tag6, Tag7, Tag8, Tag9, + Tag10, Tag11, Tag12, Tag13, Tag14, Tag15, Tag16, Tag17, Tag18, Tag19, + Tag20, Tag21, Tag22, Tag23, Tag24, Tag25, Tag26, Tag27, Tag28, Tag29, + Tag30, + + ) where + +---------------------------------------------------------- +-- import + +import Data.List (genericSplitAt, genericTake, genericReplicate, genericLength) +import Data.Bits (Bits, testBit) +import Data.Proxy (Proxy(Proxy)) +import qualified Data.Int as Z +import qualified Data.Ratio as Q +import qualified Data.Complex as C +import qualified Data.Fixed as F +import qualified Data.Word as W + +---------------------------------------------------------- +-- infix + +infix 6 :@ + +infix 6 .+ +infix 6 +. +infix 6 .- +infix 6 -. + +infix 7 .* +infix 7 *. +infix 7 /. + +infixr 8 ^. +infixr 8 ^^. +infixr 8 **. + +---------------------------------------------------------- +-- Nion + +-- | General Cayley-Dickson construction producing \"N-ions\". The +-- first parameter is a 'Tag' instance that determines the dimension, +-- which is 2 raised to 'tagVal'. The second parameter is the scalar +-- type. +data Nion n a = Scalar a | Nion n a :@ Nion n a + +---------------------------------------------------------- +-- basic operations + +-- | Equivalent to @'coord' x 0@. +scalarPart :: Nion n a -> a +scalarPart (Scalar x) = x +scalarPart (x :@ _) = scalarPart x + +-- | Equivalent to @'setCoord' x 0 0@. +purePart :: Num a => Nion n a -> Nion n a +purePart (Scalar _) = Scalar 0 +purePart (x :@ y) = purePart x :@ y + +-- | Dot product (actually the Hermitian inner product, a +-- generalization of the dot product). +dot :: Conjugable a => Nion n a -> Nion n a -> a +Scalar x `dot` Scalar y = conj x * y -- also defined as x * conj y +x@(Scalar _) `dot` (y1 :@ _) = x `dot` y1 +(x1 :@ _) `dot` y@(Scalar _) = x1 `dot` y +(x1 :@ x2) `dot` (y1 :@ y2) = (x1 `dot` y1) + (x2 `dot` y2) + +-- | Cross product. +cross :: Conjugable a => Nion n a -> Nion n a -> Nion n a +x `cross` y = y * conj x -. x `dot` y + +-- | Squared norm: the dot product of an element with itself. +sqnorm :: Conjugable a => Nion n a -> a +sqnorm x = x `dot` x + +-- | Square root of @sqnorm@. +norm :: (Conjugable a, Floating a) => Nion n a -> a +norm = sqrt . sqnorm + +-- | Promote a scalar, returning an element whose scalar part is the +-- argument and whose pure part is zero. The element behaves as if it +-- were padded with zeros, but no actual padding is done. +fromScalar :: a -> Nion n a +fromScalar = Scalar + +---------------------------------------------------------- +-- power operations + +-- | Raise to a non-negative integral power. +(^.) :: (Conjugable a, Integral b) => Nion n a -> b -> Nion n a +Scalar x ^. y = Scalar $ x ^ y +-- Copied from GHC's (^) with modifications. (c) The University of +-- Glasgow, 1994-2002. +x0 ^. y0 | y0 < 0 = error "(^.): negative exponent" + | y0 == 0 = Scalar 1 + | otherwise = f x0 y0 + where -- f : x0 ^ y0 = x ^ y + f x y | even y = f (x * x) (y `quot` 2) + | y == 1 = x + | otherwise = g (x * x) ((y - 1) `quot` 2) x + -- g : x0 ^ y0 = (x ^ y) * z + g x y z | even y = g (x * x) (y `quot` 2) z + | y == 1 = x * z + | otherwise = g (x * x) ((y - 1) `quot` 2) (x * z) + +-- | Raise to an integral power. +(^^.) :: (Conjugable a, Fractional a, Integral b) => Nion n a -> b -> Nion n a +Scalar x ^^. n = Scalar $ x ^^ n +x ^^. n | n >= 0 = x ^. n + | otherwise = recip $ x ^. negate n + +-- | Raise to a scalar power. +(**.) :: (Tag n, Conjugable a, RealFloat a) => Nion n a -> a -> Nion n a +Scalar x **. y = Scalar $ x ** y +x **. y = exp (Scalar y * log x) + +---------------------------------------------------------- +-- operations with scalars + +leftScalarOp :: (Nion n a -> Nion n a -> Nion n a) -> a -> Nion n a -> Nion n a +leftScalarOp f x y = f (Scalar x) y + +rightScalarOp :: (Nion n a -> Nion n a -> Nion n a) -> Nion n a -> a -> Nion n a +rightScalarOp f x y = f x (Scalar y) + +-- | Equivalent to @'fromScalar' x + y@. +(.+) :: Conjugable a => a -> Nion n a -> Nion n a +(.+) = leftScalarOp (+) + +-- | Equivalent to @'fromScalar' x - y@. +(.-) :: Conjugable a => a -> Nion n a -> Nion n a +(.-) = leftScalarOp (-) + +-- | Equivalent to @'fromScalar' x * y@. +(.*) :: Conjugable a => a -> Nion n a -> Nion n a +(.*) = leftScalarOp (*) + +-- | Equivalent to @x + 'fromScalar' y@. +(+.) :: Conjugable a => Nion n a -> a -> Nion n a +(+.) = rightScalarOp (+) + +-- | Equivalent to @x - 'fromScalar' y@. +(-.) :: Conjugable a => Nion n a -> a -> Nion n a +(-.) = rightScalarOp (-) + +-- | Equivalent to @x * 'fromScalar' y@. +(*.) :: Conjugable a => Nion n a -> a -> Nion n a +(*.) = rightScalarOp (*) + +-- | Equivalent to @x / 'fromScalar' y@. +(/.) :: (Conjugable a, Fractional a) => Nion n a -> a -> Nion n a +(/.) = rightScalarOp (/) + +---------------------------------------------------------- +-- polar form and complex function application + +polarUsing :: (Conjugable a, Floating a, Ord a) => + Nion n a -> Nion n a -> (a, a, Nion n a) +polarUsing sqrtMinus1 x + | sqnormp == 0 = if r >= 0 + then (r, 0, sqrtMinus1) + else (-r, pi, sqrtMinus1) + | otherwise = (absx, acos (r / absx), u) + where + r = scalarPart x + sqnormp = sqnorm x - r*r + u = purePart x /. (sqrt sqnormp) + absx = norm x + +-- | Return @(s, t, u)@ such that (approximately) +-- +-- @x == s .* 'exp' (t .* u)@ +-- +-- where @s@ and @t@ are scalars, @s >= 0@, and @u@ is a unit pure. +-- +-- If @x@ has no pure part then @u@ is arbitrarily chosen to be the +-- first pure basis element. +polar :: (Tag n, Conjugable a, Floating a, Ord a) => + Nion n a -> (a, a, Nion n a) +polar (Scalar _) = error "polar: no polar form for scalars" +polar x = polarUsing basisElement1 x + +applyUsing :: (Conjugable a, RealFloat a) => + Nion n a -> (a -> a) -> (C.Complex a -> C.Complex a) -> + Nion n a -> Nion n a +applyUsing _ f _ (Scalar s) = Scalar $ f s +applyUsing sqrtMinus1 _ f z = x .+ u *. y + where (s, t, u) = polarUsing sqrtMinus1 z + -- handle special cases for a little more accuracy + x C.:+ y | t == 0 = f s' + | t == pi = f $ (-s) C.:+ 0 -- avoid -0.0 + | otherwise = f $ s' * exp (t' * u') + where s' = s C.:+ 0 + t' = t C.:+ 0 + u' = 0 C.:+ 1 + +---------------------------------------------------------- +-- constants + +fill' :: Tag n => Proxy n -> a -> Nion n a +fill' n s = f $ tagVal n where + f 0 = Scalar s + f k = f k' :@ f k' where k' = k - 1 + +fill :: Tag n => a -> Nion n a +fill = fill' Proxy + +paddedZero :: (Tag n, Num a) => Nion n a +paddedZero = fill 0 + +validIndex :: (Tag n, Num b, Ord b) => Proxy n -> b -> Bool +validIndex n i = i >= 0 && i < 2 ^ tagVal n + +basisElement' :: (Tag n, Conjugable a, Bits i, Integral i) => + Proxy n -> i -> Nion n a +basisElement' _ 0 = Scalar 1 +basisElement' n index + | validIndex n index = setCoord paddedZero index 1 + | otherwise = error "basisElement: out of range" + +-- | The nth basis element. +basisElement :: (Tag n, Conjugable a, Bits i, Integral i) => i -> Nion n a +basisElement = basisElement' Proxy + +basisElement1 :: (Tag n, Conjugable a) => Nion n a +basisElement1 = basisElement (1 :: Integer) + +---------------------------------------------------------- +-- accessors + +coords' :: (Tag n, Num a) => Proxy n -> Nion n a -> [a] +coords' n' (Scalar x) = x : replicate (fromInteger $ 2^n - 1) 0 where + n = tagVal n' +coords' _ x = foldr (:) [] x + +-- | List of coordinates for this element. +coords :: (Tag n, Num a) => Nion n a -> [a] +coords = coords' Proxy + +coord' :: (Tag n, Num a, Integral b, Bits b) => Proxy n -> Nion n a -> b -> a +coord' _ (Scalar x) 0 = x +coord' _ (Scalar _) _ = 0 +coord' n elt index + | validIndex n index = f elt $ fromInteger $ tagVal n - 1 + | otherwise = error "coord: out of range" + where + f (Scalar x) _ = x + f (x :@ y) k = case testBit index k of + False -> f x k' + True -> f y k' + where k' = k - 1 + +-- | Get the nth coordinate. +coord :: (Tag n, Num a, Integral b, Bits b) => Nion n a -> b -> a +coord = coord' Proxy + +setCoord' :: (Tag n, Conjugable a, Num b, Bits b) => + Proxy n -> Nion n a -> b -> a -> Nion n a +setCoord' _ (Scalar _) 0 value = Scalar value +setCoord' _ (Scalar x) index value = setCoord (x .+ paddedZero) index value +setCoord' n elt index value = f elt $ fromInteger $ tagVal n - 1 where + f (Scalar _) _ = Scalar value + f (x :@ y) k = case testBit index k of + False -> f x k' :@ y + True -> x :@ f y k' + where k' = k - 1 + +-- | Set the nth coordinate, returning a new element. +setCoord :: (Tag n, Conjugable a, Num b, Bits b) => + Nion n a -> b -> a -> Nion n a +setCoord = setCoord' Proxy + +---------------------------------------------------------- +-- construction + +fromList :: Integer -> [a] -> Nion n a +fromList _ (x:[]) = Scalar x +fromList k xs = fromList k' l :@ fromList k' r where + k' = k `div` 2 + (l, r) = genericSplitAt k' xs + +nion' :: (Tag n, Num a) => Proxy n -> [a] -> Nion n a +nion' n elems = fromList d $ taken ++ padding where + d = 2 ^ tagVal n + taken = genericTake d elems + padding = genericReplicate (d - genericLength taken) 0 + +-- | Construct an element from a list of coordinates. If the list is +-- too small then the remaining coordinates are padded with zeros. If +-- the list is too large then the extra values are ignored. +nion :: (Tag n, Num a) => [a] -> Nion n a +nion = nion' Proxy + +---------------------------------------------------------- +-- instances + +instance (Tag n, Show a, Num a) => Show (Nion n a) where + show x = "nion " ++ show (coords x) + +instance (Conjugable a, Eq a) => Eq (Nion n a) where + Scalar x == Scalar y = x == y + x@(Scalar _) == y1 :@ y2 = x == y1 && y2 == 0 + x1 :@ x2 == y@(Scalar _) = x1 == y && x2 == 0 + x1 :@ x2 == y1 :@ y2 = x1 == y1 && x2 == y2 + +instance Functor (Nion n) where + fmap f (Scalar s) = Scalar $ f s + fmap f (x :@ y) = fmap f x :@ fmap f y + +instance Tag n => Applicative (Nion n) where + pure = fill + + Scalar f <*> Scalar x = Scalar $ f x + Scalar f <*> x@(_ :@ _) = pure f <*> x + f@(_ :@ _) <*> (Scalar x) = f <*> pure x + (f1 :@ f2) <*> (x1 :@ x2) = (f1 <*> x1) :@ (f2 <*> x2) + +instance Foldable (Nion n) where + foldr f acc (Scalar x) = f x acc + foldr f acc (x :@ y) = foldr f (foldr f acc y) x + +instance Traversable (Nion n) where + traverse f (Scalar x) = Scalar <$> (f x) + traverse f (x :@ y) = (:@) <$> traverse f x <*> traverse f y + +instance Conjugable a => Num (Nion n a) where + Scalar x + Scalar y = Scalar $ x + y + x@(Scalar _) + (y1 :@ y2) = (x + y1) :@ y2 + (x1 :@ x2) + y@(Scalar _) = (x1 + y) :@ x2 + (x1 :@ y1) + (x2 :@ y2) = (x1 + x2) :@ (y1 + y2) + + Scalar x - Scalar y = Scalar $ x - y + x@(Scalar _) - (y1 :@ y2) = (x - y1) :@ negate y2 + (x1 :@ x2) - y@(Scalar _) = (x1 - y) :@ x2 + (x1 :@ y1) - (x2 :@ y2) = (x1 - x2) :@ (y1 - y2) + + Scalar x * Scalar y = Scalar $ x * y + x@(Scalar _) * (y1 :@ y2) = (x * y1) :@ (x * y2) + (x1 :@ x2) * y@(Scalar _) = (x1 * y) :@ (x2 * y) + (x1 :@ x2) * (y1 :@ y2) = (x1 * y1 - conj y2 * x2) :@ (y2 * x1 + x2 * conj y1) + + negate = fmap negate + fromInteger = fromScalar . fromInteger + abs = doNotUse + signum = doNotUse + +instance (Conjugable a, Fractional a) => Fractional (Nion n a) where + Scalar x / Scalar y = Scalar $ x / y + x@(Scalar _) / y@(_ :@ _) = x * recip y + x@(_ :@ _) / Scalar y = fmap (/ y) x + x@(_ :@ _) / y@(_ :@ _) = (x * conj y) /. sqnorm y + + recip x = conj x /. sqnorm x + fromRational = fromScalar . fromRational + +-- | The first pure basis element is arbitrarily chosen as sqrt (-1). +instance (Tag n, Conjugable a, RealFloat a) => + Floating (Nion n a) where + pi = Scalar pi + exp = applyUsing basisElement1 exp exp + log = applyUsing basisElement1 log log + sqrt = applyUsing basisElement1 sqrt sqrt + sin = applyUsing basisElement1 sin sin + cos = applyUsing basisElement1 cos cos + tan = applyUsing basisElement1 tan tan + asin = applyUsing basisElement1 asin asin + acos = applyUsing basisElement1 acos acos + atan = applyUsing basisElement1 atan atan + sinh = applyUsing basisElement1 sinh sinh + cosh = applyUsing basisElement1 cosh cosh + tanh = applyUsing basisElement1 tanh tanh + asinh = applyUsing basisElement1 asinh asinh + acosh = applyUsing basisElement1 acosh acosh + atanh = applyUsing basisElement1 atanh atanh + +---------------------------------------------------------- +-- convenience types + +type Complex a = Nion Tag1 a +type Quaternion a = Nion Tag2 a +type Octonion a = Nion Tag3 a +type Sedenion a = Nion Tag4 a + +-- | Construct a complex number. +complex :: a -> a -> Complex a +complex x y = (:@) (Scalar x) (Scalar y) + +-- | Construct a quaternion. +quaternion :: a -> a -> a -> a -> Quaternion a +quaternion w x y z = (:@) ((:@) (Scalar w) (Scalar x)) + ((:@) (Scalar y) (Scalar z)) + +-- | Construct an octonion. +octonion :: a -> a -> a -> a -> + a -> a -> a -> a -> Octonion a +octonion s t u v + w x y z = (:@) ((:@) ((:@) (Scalar s) (Scalar t)) + ((:@) (Scalar u) (Scalar v))) + ((:@) ((:@) (Scalar w) (Scalar x)) + ((:@) (Scalar y) (Scalar z))) + +-- | Construct a sedenion. +sedenion :: a -> a -> a -> a -> + a -> a -> a -> a -> + a -> a -> a -> a -> + a -> a -> a -> a -> Sedenion a +sedenion k l m n + o p q r + s t u v + w x y z = (:@) ((:@) ((:@) ((:@) (Scalar k) (Scalar l)) + ((:@) (Scalar m) (Scalar n))) + ((:@) ((:@) (Scalar o) (Scalar p)) + ((:@) (Scalar q) (Scalar r)))) + ((:@) ((:@) ((:@) (Scalar s) (Scalar t)) + ((:@) (Scalar u) (Scalar v))) + ((:@) ((:@) (Scalar w) (Scalar x)) + ((:@) (Scalar y) (Scalar z)))) + +---------------------------------------------------------- +-- Conjugable + +class Num a => Conjugable a where + -- | Conjugate. + conj :: a -> a + +instance Conjugable a => Conjugable (Nion n a) where + conj (Scalar x) = Scalar $ conj x + conj (x :@ y) = conj x :@ negate y + +instance RealFloat a => Conjugable (C.Complex a) where + conj = C.conjugate + +instance Conjugable Int where conj = id +instance Conjugable Integer where conj = id +instance Conjugable Float where conj = id +instance Conjugable Double where conj = id +instance Conjugable Z.Int8 where conj = id +instance Conjugable Z.Int16 where conj = id +instance Conjugable Z.Int32 where conj = id +instance Conjugable Z.Int64 where conj = id +instance Conjugable W.Word8 where conj = id +instance Conjugable W.Word16 where conj = id +instance Conjugable W.Word32 where conj = id +instance Conjugable W.Word64 where conj = id +instance Integral a => Conjugable (Q.Ratio a) where conj = id +instance F.HasResolution a => Conjugable (F.Fixed a) where conj = id + +----------------------------------------------------------------------------- +-- doNotUse + +rant :: String +rant = unlines $ + ["", + "The Num class is a bit messed up, having tied (+), (-), and (*) to abs", + "and signum. Number systems that have no appropriate definition for abs", + "or signum must either invent their own operators for addition,", + "subtraction, and multiplication, else break the contract with Num by", + "raising an error such as this one when someone uses abs or signum.", + "", + "For some time I resisted hijacking Num, but eventually the replacement", + "operators became too cumbersome and, coupled with the lack of numeric", + "promotion, significantly detracted from the usability of the package.", + "So here we are. Good luck, and stay away from abs and signum, which", + "officially have cooties."] + +doNotUse :: a -> a +doNotUse _ = error rant + +---------------------------------------------------------- +-- Tag + +class Tag n where + tagVal :: Proxy n -> Integer + +data Tag0 +data Tag1 +data Tag2 +data Tag3 +data Tag4 +data Tag5 +data Tag6 +data Tag7 +data Tag8 +data Tag9 +data Tag10 +data Tag11 +data Tag12 +data Tag13 +data Tag14 +data Tag15 +data Tag16 +data Tag17 +data Tag18 +data Tag19 +data Tag20 +data Tag21 +data Tag22 +data Tag23 +data Tag24 +data Tag25 +data Tag26 +data Tag27 +data Tag28 +data Tag29 +data Tag30 + +instance Tag Tag0 where tagVal _ = 0 +instance Tag Tag1 where tagVal _ = 1 +instance Tag Tag2 where tagVal _ = 2 +instance Tag Tag3 where tagVal _ = 3 +instance Tag Tag4 where tagVal _ = 4 +instance Tag Tag5 where tagVal _ = 5 +instance Tag Tag6 where tagVal _ = 6 +instance Tag Tag7 where tagVal _ = 7 +instance Tag Tag8 where tagVal _ = 8 +instance Tag Tag9 where tagVal _ = 9 +instance Tag Tag10 where tagVal _ = 10 +instance Tag Tag11 where tagVal _ = 11 +instance Tag Tag12 where tagVal _ = 12 +instance Tag Tag13 where tagVal _ = 13 +instance Tag Tag14 where tagVal _ = 14 +instance Tag Tag15 where tagVal _ = 15 +instance Tag Tag16 where tagVal _ = 16 +instance Tag Tag17 where tagVal _ = 17 +instance Tag Tag18 where tagVal _ = 18 +instance Tag Tag19 where tagVal _ = 19 +instance Tag Tag20 where tagVal _ = 20 +instance Tag Tag21 where tagVal _ = 21 +instance Tag Tag22 where tagVal _ = 22 +instance Tag Tag23 where tagVal _ = 23 +instance Tag Tag24 where tagVal _ = 24 +instance Tag Tag25 where tagVal _ = 25 +instance Tag Tag26 where tagVal _ = 26 +instance Tag Tag27 where tagVal _ = 27 +instance Tag Tag28 where tagVal _ = 28 +instance Tag Tag29 where tagVal _ = 29 +instance Tag Tag30 where tagVal _ = 30 diff --git a/test/test.hs b/test/test.hs new file mode 100644 --- /dev/null +++ b/test/test.hs @@ -0,0 +1,422 @@ +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] + +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 + +checkApplicative :: IO () +checkApplicative = do + let x = quaternion 1 2 3 4 :: Quaternion Integer + y = quaternion 5 6 7 8 :: Quaternion Integer + r = nion [3] :: Nion Tag0 Integer + s = nion [4] :: Nion Tag0 Integer + assert $ ((+) <$> x <*> y) == x + y + assert $ ((-) <$> x <*> y) == x - y + assert $ ((*) <$> 3 <*> x) == 3 * x + assert $ ((*) <$> x <*> 3) == x * 3 + assert $ ((*) <$> r <*> s) == r * s + +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 + checkApplicative + checkPower + checkInverses + checkMixed + replicateM_ 20 $ do + checkFloating + checkProperties1 + checkProperties2 + checkProperties3 + checkProperties4 + putStrLn "\nAll tests passed."