-- 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 errors in this calculator:
-- the property a^b * a^c == a^(b+c) doesn't hold.
-- (e.g., 
-- 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 tests
-- although it does have large tests cases some of them seg faults sometimes
-- depending on the machine I built this ord calcualtor.



module Math.Ordinals.MultiSet where

import Data.List (groupBy, intersperse)

newtype Ordinal = O [Ordinal] deriving Eq

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 (O [])            = True
wf (O [o])           = wf o
wf (O (o:os@(o':_))) = o >= o' && wf (O os)


instance Enum Ordinal where -- TODO
  toEnum = fromInt
  fromEnum = toInt

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"

  div a b = fst (divMod a b)

  mod a b = snd (divMod a b)

  divMod = quotRem -- 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)

  -- 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

  -- 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 : define new op since neither ^ nor ^^ will work
infixr 8 ^:
_                ^:  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)