{- |
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
 -- '+' (addition), '*' (multiplication), '-', (left addtion)
 -- 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
  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"

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

  -- | Addition.
  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)