    -- Copyright (c) 2011-2015, David Amos. All rights reserved.

    {-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, TypeSynonymInstances, NoMonomorphismRestriction #-}

    -- |A module defining the (non-associative) algebra of octonions over an arbitrary field.
    --
    -- The octonions are the algebra defined by the basis {1,i0,i1,i2,i3,i4,i5,i6},
    -- where each i_n * i_n = -1, and i_n+1 * i_n+2 = i_n+4 (where the indices are modulo 7).
    module Math.Algebras.Octonions where

    import Prelude hiding ( (*>) )

    import Math.Core.Field
    import Math.Algebras.VectorSpace
    import Math.Algebras.TensorProduct hiding (i1,i2)
    import Math.Algebras.Structures
    import Math.Algebras.Quaternions

    import Math.Combinatorics.FiniteGeometry (ptsAG)

    -- Conway & Smith, On Quaternions and Octonions

    -- OCTONIONS

    data OBasis = O Int deriving (Eq,Ord)
    -- map (return . O) [-1..6] -> [1,i0,i1,i2,i3,i4,i5,i6]

    type Octonion k = Vect k OBasis

    instance Show OBasis where
        show (O n) | n == -1 = "1"
                   | 0 <= n && n <= 6 = "i" ++ show n
                   | otherwise = error "Octonion: invalid basis element"

    i0, i1, i2, i3, i4, i5, i6 :: Octonion Q
    i0 = return (O 0)
    i1 = return (O 1)
    i2 = return (O 2)
    i3 = return (O 3)
    i4 = return (O 4)
    i5 = return (O 5)
    i6 = return (O 6)

    i_ :: Num k => Int -> Octonion k
    i_ n = return (O n)

    instance (Eq k, Num k) => Algebra k OBasis where
        unit x = x *> return (O (-1))
        mult = linear m where
            m (O (-1), O n) = return (O n)
            m (O n, O (-1)) = return (O n)
            m (O a, O b) = case (b-a) `mod` 7 of
                           0 -> -1
                           1 -> i_ ((a+3) `mod` 7)       -- i_n+1 * i_n+2 == i_n+4
                           2 -> i_ ((a+6) `mod` 7)       -- i_n+2 * i_n+4 == i_n+1
                           3 -> -1 *> i_ ((a+1) `mod` 7) -- i_n+1 * i_n+4 == -i_n+2
                           4 -> i_ ((a+5) `mod` 7)       -- i_n+4 * i_n+1 == i_n+2
                           5 -> -1 *> i_ ((a+4) `mod` 7) -- i_n+4 * i_n+2 == -i_n+1
                           6 -> -1 *> i_ ((a+2) `mod` 7) -- i_n+2 * i_n+1 == -i_n+4

    instance (Eq k, Num k) => HasConjugation k OBasis where
        conj = (>>= conj') where
            conj' (O n) = (if n == -1 then 1 else -1) *> return (O n)
        -- ie conj = linear conj', but avoiding unnecessary nf call
        sqnorm x = sum $ map ((^2) . snd) $ terms x
        -- sqnorm x = scalarPart (x * conj x)

    -- Hence, the octonions inherit a Fractional instance

    -- octonions fq = [sum $ zipWith (\x n -> x *> i_ n) xs [-1..6] | xs <- ptsAG 8 fq]