{-# LANGUAGE ScopedTypeVariables #-}
import qualified Data.Set as S
import Data.Set ((\\), powerSet, Set, cartesianProduct)
import qualified Data.Array as A
type Geo a = Set a -> Float
memo :: forall a.(Bounded a, Enum a, Ord a) => Geo a -> Geo a
memo v =
let upperBound = 2 ^ (1 fromEnum (maxBound :: a)) - 1
basisToPosition = foldr (\x -> ( 2^fromEnum x)) 0 :: [a] -> Int
arr = A.array (0, upperBound)
[ (basisToPosition (S.toList b), v b)
| b <- S.toList (powerSet pseudoscalarBasis)
]
in \b -> arr A.! basisToPosition (S.toList b)
(โณ) :: Ord a => Set a -> Set a -> Set a
a โณ b = S.union (a \\ b) (b \\ a)
pseudoscalarBasis :: (Enum a, Bounded a, Ord a) => Set a
pseudoscalarBasis = S.fromList [minBound .. maxBound]
sign :: (Ord a) => Set a -> Set a -> Float
sign u v =
product $ liftA2
(\a b -> case compare a b of {LT -> 1; EQ -> 0; GT -> -1})
(S.toList u)
(S.toList v)
dual :: (Enum a, Bounded a, Ord a) => Geo a -> Geo a
dual f s = (f * sign s) (pseudoscalarBasis โณ s)
(/\) :: (Ord a, Enum a, Bounded a) => Geo a -> Geo a -> Geo a
(/\) f g =
memo $ \t -> sum
[ sign tau (t โณ tau) * f tau * g (t โณ tau)
| tau <- S.toList (powerSet pseudoscalarBasis)
]
data Basis = X | Y | Z | W deriving (Enum, Bounded, Show, Eq, Ord)
x, y, z, w :: Geo Basis
x b = if b == S.singleton X then 1 else 0
y b = if b == S.singleton Y then 1 else 0
z b = if b == S.singleton Z then 1 else 0
w b = if b == S.singleton W then 1 else 0
instance Num (Geo a) where
( ) = liftA2 ( )
(*) = liftA2 (*)
fromInteger = const . fromInteger
negate = fmap negate
abs = fmap abs
signum = fmap signum
instance Fractional (Geo a) where
(/) = liftA2 (/)
fromRational = const . fromRational
main = do
let system = (3*x y-5*w) /\ (y z-4*w) /\ (x-y 2*z-8*w)
solution = dual system
wValue = solution (S.singleton W)
homogenizedSolution = solution / realToFrac wValue
putStrLn ("x = " show (homogenizedSolution (S.singleton X)))
putStrLn ("y = " show (homogenizedSolution (S.singleton Y)))
putStrLn ("z = " show (homogenizedSolution (S.singleton Z)))