```{- |
Encoding of ordinals up to epsilon_0 as an iterated multiset:
definition in Basic Proof Theory by Troelstra and Schwichenberg.
Note, this is not the most efficient way to calculate ordinals.
This library is better than having none.
I think CNF representation would be more efficent,
planning to add in the next version of this library.

For further analysis on efficiency of implementations on ordinals see
<http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.91.8089>

FYI, an ordinal calculator that covers wider range beyond epsion_0
can be found at <http://www.mtnmath.com/ord/> which is written C++.
However, I found a serious error in this calculator (ver 0.2):
the property a^b * a^c == a^(b+c) does not hold.  For example,

@
ordCalc> 3 ^(w^w + w + 1) * 3^(w^w) == 3 ^ (w^w + w + 1 + w^w)
FALSE
@

This calculator didn't seem to run QuickCeck style automatic
testing, although it does have hundreds (or maybe more than a thousand)
tests cases but some of them even causes segmentation faults
depending on the machine I built this ordCalc.
-}
module Math.Ordinals.MultiSet (
-- * Types
Ordinal(..)

-- * Operators
, (^:)
-- | There are more operators such as
-- defined in the 'Num' class instance declaration for 'Ordinal'.
-- Although ordinals are not really 'Num' we decided to make it
-- as a 'Num' instance for convenience;  We can use functions such
-- as 'sum' and 'product' over list of ordinals and they behave well.
-- We also plan to implement division and remainder operations
-- in the 'Integral' class instance for similar reason.  For further
-- information on class instances see the source code.

-- * Auxiliary functions
, w, wf

-- * Auxiliary operators
-- | These operators can manipulate the 'Ordinal' newtype data structure
-- internally.  Use with care since it can break the well-formedness ('wf')
-- of ordinal representation.
, (.:), (++.), (.++.)
) where

import Data.List (groupBy, intersperse)

newtype Ordinal = O [Ordinal] deriving Eq

-- | convenience function that takes an argument as the power of omega
-- (the first limit ordinal).
w :: Ordinal -> Ordinal
w o = O [o]

instance Show Ordinal where
show o@(O os)
| o < O[1]  = show (length os)
| otherwise = "("
++ (concat . intersperse " + ")
["w "++ show o' ++ showFactor(length os') | os'@(o':_)<-oss ]
++ ")"
where showFactor k = if k==1 then "" else (" *"++show k)
oss = groupBy (==) os

instance Ord Ordinal where
compare 0         0         = EQ
compare 0         _         = LT
compare _         0         = GT
compare (O(a:as)) (O(b:bs)) = case compare a b of
EQ -> compare as bs
r  -> r

-- | well formedness of ordinals
wf :: Ordinal -> Bool
wf (O [])            = True
wf (O [o])           = wf o
wf (O (o:os@(o':_))) = o >= o' && wf (O os)

instance Enum Ordinal where -- TODO
toEnum = fromInt

instance Real Ordinal where -- TODO

instance Integral Ordinal where
-- works only within Int limit
toInteger o@(O as) | all (0 ==) as = fromIntegral (toInt o)
toInteger _                        = error "ordinal not less than omega"

-- | Left division. not yet defined
div a b = fst (divMod a b)

-- | not yet defined
mod a b = snd (divMod a b)

-- | not yet defined
divMod = quotRem

-- | not yet defined -- TODO -- somebody help figure this out
quotRem _ _ = error "Ordinal quotRem not yet defined"

instance Num Ordinal where
-- to parse in haskell integer literals within Int limit
fromInteger k = fromInt (fromIntegral k)

abs = id

signum 0         = 0
signum (O(o:os)) = O[o]

o    + 0          = o
O as + O bs@(b:_) = O (takeWhile (>=b) as ++ bs)

-- | Left subtraction
-- for a <= b exists r = b - a such that a + r = b
-- i.e., a + (b - a) = b for a <= b
o            - 0            = o
0            - _            = 0
o1@(O(a:as)) - o2@(O(b:bs)) = case compare a b of
LT -> 0
EQ -> O as - O bs
GT -> o1

-- | Multiplication.
-- Implemented as something similar to
-- <http://www.volny.cz/behounek/logic/papers/ordcalc/index.html>
0    * _                   = 0
o1@(O(a:as)) * O bs
| null bs'  = sum [o1 | _ <- zs]     -- finite
| null zs   = O [a+b | b<-bs']       -- limit ordinal
| otherwise = o1 * O bs' + o1 * O zs -- non-limit ordinal
where (bs',zs) = span (>0) bs
(as1,as2) = span (a==) (a:as)

-- | Exponentiation.
-- Defined a new operator since neither '^' nor '^^' will work.
-- Note, @(w o)@ is same as @(w 1 :^ o)@ for any oridnal @o@.
infixr 8 ^:
(^:) :: Ordinal -> Ordinal -> Ordinal
_                ^:  0                = 1
0                ^:  _                = 0
1                ^:  _                = 1
o1@(O osA@(a:_)) ^:  o2@(O osB@(b:_))
| finiteB   = o1^: (o2 - 1) * o1       -- finite power
| mB > 0    = o1^: O osB'  *  o1^: mB  -- non-limit ordinals
| finiteA   = O [o2]
| otherwise = O [a * o2]
where
(osA',zsA) = span (>0) osA
(osB',zsB) = span (>0) osB
mA = O zsA
mB = O zsB
finiteB = null osB'
finiteA = null osA'

toInt (O as) | all (0 ==) as = length as
toInt _                      = error "ordinal not less than omega"

fromInt k = O (replicate k (O []))

infixr 5 .:
o .: O os = O (o:os)
infixr 5 ++.
as ++. O bs = O (as++bs)
infixr 5 .++.
O as .++. O bs = O (as++bs)

```