-- Haskell98! (modulo so convenient and accepted in H98 mode pattern guards)

-- Total stream processors and the quantification over all possible
-- infinite streams.

-- This is the stream formulation of the article by
-- Marti'n Escardo' ``seemingly impossible functional program''.
-- The stream formulation reveals the mystery of deciding the satisfiability
-- of a total computable predicate over the set of all Cantor numerals,
-- (The cardinality of that set is greater than that of natural numbers!)
--
-- As an additional contribution, we show how to write functions over Cantor
-- numerals in a `natural' monadic style so that those functions become
-- self-partially evaluating. The instantiation of the functions in an
-- appropriate (pure, neither ST nor IO) monad gives us transparent 
-- memoization, without any changes to the functions themselves. The monad
-- in question is pure and involves no reference cells.

-- On `dense' functions on numerals (i.e., those that need to examine
-- most of the bits of its argument, up to a limit), our technique 
-- performs about 9 times better than the most sophisticated
-- one by Marti'n Escardo'.

-- See the accompanying article for explanation.

module StreamPEval where


import qualified Data.IntMap as M
import Control.Monad
import Debug.Trace


data Bit = Zero | One deriving (Eq,Show)

type Natural = Int

-- from Marti'n Escardo's article
coerce :: Bit -> Natural
coerce Zero = 0
coerce One = 1

-- Representation for Cantor numerals, infinite streams of bits
-- A Cantor numeral consists of a finite prefix and a function to
-- extend it
type Prefix = M.IntMap Bit
data Cantor = Cantor Prefix (Prefix -> Cantor)

-- sample Cantor numbers
-- One may regard a Cantor number as a _real_ number within [0,1]

zeros, ones :: Natural -> (Prefix -> Cantor)
zeros i pr = Cantor (M.insert i Zero pr) (zeros (succ i))
ones i pr  = Cantor (M.insert i One pr) (ones (succ i))

-- All zeros and all ones
zero, one:: Cantor
zero = Cantor (M.singleton 0 Zero) (zeros 1)
one  = Cantor (M.singleton 0 One) (ones 1)

-- Alternating ones and zeros
twothirds = twothirds1 0 M.empty
 where
 twothirds0 i pr = Cantor (M.insert i Zero pr) (twothirds1 (succ i))
 twothirds1 i pr = Cantor (M.insert i One pr)  (twothirds0 (succ i))


-- Give the finite prefix of a Cantor number (useful for printing)
ctake :: Natural -> Cantor -> [Natural]
ctake n c = ctake' 0 c
 where
 ctake' i _ | i >= n = []
 ctake' i c@(Cantor prefix _) | Just v <- M.lookup i prefix
    = coerce v : ctake' (succ i) c
 ctake' i c@(Cantor prefix step) = ctake' i (step prefix)



-- Functions on Cantor numbers
-- Given a finite prefix, the function returns the value or reports `failure' 
-- (hence the monadic type). The `failure' means the given data are not 
-- sufficient to compute the result and so the function ought to be
-- called with a longer prefix.
type CantorFN m a = Prefix -> m a


-- A computable function guarantees that it successfully computes
-- the value once the prefix becomes long enough.
-- Note, some useful functions aren't in this class, for example, equality.

-- To make monadic expressions easier to write, we introduce a bit of
-- sugar. First, we assume (m a) is in Num provided that a is.

-- Attempt to extract the i-th bit from the prefix
class Monad m => GetBit m where
    getbit :: Prefix -> Int -> m Natural


-- Sample functions on Cantor numbers (from Marti'n Escardo's article)
-- Our notation is essentially the same as his (given in the comments)

-- proj :: Natural -> (Cantor -> Integer)
-- proj i = \a -> coerce(a i)
proj i = \pr -> getbit pr i

-- f, g, h :: Cantor -> Integer
f,g,h :: (GetBit m, Num (m Natural)) => CantorFN m Natural

-- f a = coerce(a(7 * coerce(a 4) +  4 * (coerce(a 7)) + 4))
f a = getbit a =<< (7 * (getbit a 4) +  4 * (getbit a 7) + 4)

-- g a = coerce(a(coerce(a 4) + 11 * (coerce(a 7))))
g a = getbit a =<< (getbit a 4 + 11 * (getbit a 7))

-- h a = if a 7 == Zero
--       then if a 4 == Zero then coerce(a  4) else coerce(a 11)
--       else if a 4 == One  then coerce(a 15) else coerce(a  8)


h a = ifm ((== 0) `liftM` getbit a 7)
          (ifm ((== 0) `liftM` getbit a 4) (getbit a 4) (getbit a 11))
          (ifm ((== 1) `liftM` getbit a 4) (getbit a 15) (getbit a 8))

ifm test tb eb = do b <- test; if b then tb else eb


-- The following function is not present in Marti'n Escardo's article
-- It converts the prefix of width w of a Cantor number to a regular number
-- It is, too, total, and shows an important corner case, when a function
-- uses all of the bits from its argument, up to a certain limit.

-- ton w a = foldl (\ac i -> 2*ac + coerce (a i)) 0 [0..w-1]
ton w a = foldl (\ac i -> 2*ac + getbit a i) 0 [0..w-1]


-- The first approach: the monad in question is Maybe
-- Just x indices success and Nothing indicates failure.

instance Num a => Num (Maybe a) where
    (+) = liftM2 (+)
    (-) = liftM2 (-)
    (*) = liftM2 (*)
    fromInteger = return . fromInteger


instance GetBit Maybe where
    getbit pr i = liftM coerce (M.lookup i pr)
{-
    getbit pr i = 
	maybe (trace ("no bit " ++ show i) Nothing) 
	      (\b -> let b' = coerce b
		     in trace ("bit " ++ show i ++ ": " ++ show b') (Just b'))
	      (M.lookup i pr)
-}

-- Applying a CantorFN to a Cantor number: feeding a CantorFN
-- a progressively longer prefix until the function succeeds
capply1 :: CantorFN Maybe a -> Cantor -> a
capply1 f (Cantor pr next) = maybe (capply1 f (next pr)) id (f pr)
	       

-- Quantification over Cantor numerals
-- Given a total computable predicate, CantorFN Maybe Bool,
-- check to see if it is satisfiable. If so, return the witness,
-- a Cantor numeral that makes the predicate succeed with a True value.
-- By definition, the predicate succeeds on a finite prefix of a Cantor
-- number; the rest of the number is irrelevant. We just set it to be all
-- zeros.

-- The following predicate uses Depth-first search, with no memoization.

search1 :: CantorFN Maybe Bool -> Maybe Cantor
search1 p = s0 M.empty 0
 where s0 pr i = let pr' = M.insert i Zero pr
		     i'  = succ i
		 in case p pr' of
		     Just True  -> Just $ Cantor pr' (zeros i')
		     Nothing | c@(Just _) <- s0 pr' i' -> c
		     _ -> s1 pr i
       s1 pr i = let pr' = M.insert i One pr
		     i'  = succ i
		 in case p pr' of
		     Just True  -> Just $ Cantor pr' (zeros i')
		     Nothing    -> s0 pr' i'
		     Just False -> Nothing

-- conventional quantifiers
forsome1,forevery1 :: CantorFN Maybe Bool -> Bool
forsome1  p  = maybe False (const True) (search1 p)
forevery1 p = not (forsome1 (liftM not . p))

-- CantorFN m a admit decidable extensional equality (provided 'a' admits it)

-- equal :: Eq y => (Cantor -> y) -> (Cantor -> y) -> Bool
-- equal f g = forevery(\a -> f a == g a)

-- equal :: Eq y => (Cantor -> y) -> (Cantor -> y) -> Bool
equal1 f g = forevery1(\a -> liftM2 (==) (f a) (g a))


-- tests

-- The following tests look better with the tracing enabled
test11 = capply1  (ton 3) twothirds
test12 = forsome1 (liftM (== 21) . (ton 3))
test13 = equal1   (ton 3) (ton 3)
test14 = forsome1 (liftM (== 1) . f)

-- The following are stress tests. It's better to disable the tracing
testr11 = equal1 f g
-- (0.01 secs, 1050148 bytes)
testr12 = equal1 f h
-- (0.86 secs, 47096912 bytes)
testr13 = equal1 f f
-- (0.80 secs, 50578068 bytes)
testr14 = equal1 g g
-- (0.13 secs, 6407772 bytes)
testr15 = equal1 h h
-- (0.95 secs, 41502680 bytes)

-- Real stress tests
-- The tests13-15 show the exponential behavior

tests11 = equal1 (proj 14) (proj 14)
-- Marti'n Escardo's find_i result: (1.82 secs, 39783148 bytes)
-- ours: (0.42 secs, 22456632 bytes)
tests12 = equal1 (proj 14) (proj 15)
-- Marti'n Escardo's find_i result: (3.64 secs, 79934000 bytes)
-- ours: (0.00 secs, 0 bytes)
tests13 = equal1 (ton 14) (ton 14)
-- Marti'n Escardo's find_i result: (5.98 secs, 153746524 bytes)
-- ours: (2.96 secs, 245394268 bytes)
tests14 = equal1 (ton 14) (ton 15)
-- Marti'n Escardo's find_i result: (12.85 secs, 325268744 bytes)
-- ours: (0.00 secs, 0 bytes)
tests15 = equal1 (ton 15) (ton 15)
-- Marti'n Escardo's find_i result: (13.61 secs, 342230580 bytes)
-- ours: (6.03 secs, 524111536 bytes)

-- The second attempt: a self-partially evaluating monad
-- In the case of `failure', it returns potentially a new function
-- to apply to a bigger prefix.
data Inc a = Done a | More (Prefix -> Inc a)

instance Monad Inc where
    return = Done
    (Done x) >>= k = k x
    (More f) >>= k = More (\pr -> f pr >>= k)

instance Eq (Inc a)
instance Show (Inc a)

instance Num a => Num (Inc a) where
    (+) = liftM2 (+)
    (-) = liftM2 (-)
    (*) = liftM2 (*)
    fromInteger = return . fromInteger

instance GetBit Inc where
-- {-
    getbit pr i = 
	maybe (More (\pr -> getbit pr i)) (return . coerce) (M.lookup i pr)
-- -}
{-
    getbit pr i = 
	maybe (trace ("no bit " ++ show i) (More (\pr -> getbit pr i)))
	      (\b -> let b' = coerce b
		     in trace ("bit "++ show i ++ ": " ++ show b') (return b'))
	      (M.lookup i pr)
-}

-- Applying a CantorFN to a Cantor number: feeding a CantorFN
-- a progressively longer prefix until the function succeeds
-- Please compare with capply1
capply2 :: CantorFN Inc a -> Cantor -> a
capply2 f (Cantor pr next) = case f pr of
			       Done x -> x
			       More f -> capply2 f (next pr)

-- The version of search1 with memoization...
search2 :: CantorFN Inc Bool -> Maybe Cantor
search2 p = s0 p M.empty 0
 where s0 p pr i = let pr' = M.insert i Zero pr
		       i'  = succ i
		   in case p pr' of
		       Done True  -> Just $ Cantor pr' (zeros i')
		       More p | c@(Just _) <- s0 p pr' i' -> c
		       _ -> s1 p pr i
       s1 p pr i = let pr' = M.insert i One pr
		       i'  = succ i
		   in case p pr' of
		       Done True  -> Just $ Cantor pr' (zeros i')
		       More p     -> s0 p pr' i'
		       _          -> Nothing


forsome2 p  = maybe False (const True) (search2 p)
forevery2 p = not (forsome2 (liftM not . p))

-- One may think we can use the same formula for equality:
-- equal2 f g = forevery2(\a -> liftM2 (==) (f a) (g a))
-- Alas, it does not take benefit of memoization.
-- We need a special version of liftM2 for our monad, below:

combine op f g = \pr -> case (f pr,g pr) of
                         (Done fv, Done gv) -> Done (op fv gv)
			 (More f, Done gv)  -> More (liftM (flip op gv) . f)
			 (Done fv, More g)  -> More (liftM (op fv) . g)
			 (More f, More g)   -> More (combine op f g)

equal2 f g = forevery2 (combine (==) f g)


-- tests

-- The following tests look better with the tracing enabled
test21 = capply2  (ton 3) twothirds
test22 = forsome2 (liftM (== 21) . (ton 3))
test23 = equal2   (ton 3) (ton 3)
test24 = forsome2 (liftM (== 1) . f)

-- The following are stress tests. It's better to disable the tracing
testr21 = equal2 f g
-- (0.00 secs, 0 bytes)
testr22 = equal2 f h
-- (0.41 secs, 15737072 bytes)
testr23 = equal2 f f
-- (0.40 secs, 15683040 bytes)
testr24 = equal2 g g
-- (0.08 secs, 2488052 bytes)
testr25 = equal2 h h
-- (0.40 secs, 15709976 bytes)

-- Real stress tests
-- The tests23-25 show the exponential behavior

tests21 = equal2 (proj 14) (proj 14)
-- Marti'n Escardo's find_i result: (1.82 secs, 39783148 bytes)
-- Ours: (0.72 secs, 27340500 bytes)
tests22 = equal2 (proj 14) (proj 15)
-- Marti'n Escardo's find_i result: (3.64 secs, 79934000 bytes)
-- Ours: (0.00 secs, 0 bytes)
tests23 = equal2 (ton 14) (ton 14)
-- Marti'n Escardo's find_i result: (5.98 secs, 153746524 bytes)
-- Ours: (0.70 secs, 27485480 bytes)
tests24 = equal2 (ton 14) (ton 15)
-- Marti'n Escardo's find_i result: (12.85 secs, 325268744 bytes)
-- ours: (0.00 secs, 0 bytes)
tests25 = equal2 (ton 15) (ton 15)
-- Marti'n Escardo's find_i result: (13.61 secs, 342230580 bytes)
-- ours: (1.41 secs, 55052820 bytes)

tests26 = equal2 (ton 17) (ton 17)
-- (5.60 secs, 218557248 bytes)

-- Incidentally, even with the most sophisticated find_vii,
-- Marti'n Escardo's code
-- equal (ton 14) (ton 14) still takes
-- (6.87 secs, 186692916 bytes)
-- Thus, our technique performs about 9 times faster on `dense' functions.