-- Programs discussed in the paper -- "From coinductive proofs to exact real arithmetic" -- Contents -- 1. Natural numbers and the parity example -- 2. Representation of real numbers in [-1,1] -- 3. Memoized representation of u.c. functions -- 4. Wellfounded induction/recursion and digital systems -- 5. Example: linear affine maps f_(u,v)(x) = u*x+v -- 6. Example: quadratic maps f_(u,v,w)(x) = u*x^2+v*x+w -- 7. Example: logistic maps and their iterations -- 8. Auxiliary definitions -- 9. Programs following strictly the scheme of the last part of Sect. 3 (of the paper) -- 1. Natural numbers and the parity example -- ========================================= data Nat = Zero | Succ Nat -- data deriving Show -- "deriving Show" enables a default display parity :: Nat -> Bool parity Zero = True parity (Succ x) = case parity x of {True -> False ; False -> True} parity1 :: Nat -> (Nat,Bool) parity1 Zero = (Zero,True) parity1 (Succ x) = case parity1 x of {(y,True) -> (y,False) ; (y,False) -> (Succ y,True)} -- 2. Representation of real numbers in [-1,1] -- ========================================== -- signed digits data SD = N | Z | P deriving Show -- N = -1 , Z = 0, P = 1 type J0 alpha = (SD,alpha) -- functoriality of J0 mapJ0 :: (alpha -> beta) -> J0 alpha -> J0 beta mapJ0 f (d,x) = (d,f x) -- signed digit streams data C0 = ConsC0 (J0 C0) deriving Show -- codata -- coiteration coitC0 :: (alpha -> J0 alpha) -> alpha -> C0 coitC0 s x = ConsC0 (mapJ0 (coitC0 s) (s x)) -- this is equivalent to -- coitC0 s x = ConsC0 (d,coitC0 s y) where (d,y) = s x -- extracted programs sd2Rational :: SD -> Rational sd2Rational d = case d of {N -> -1 ; Z -> 0 ; P -> 1} cauchy2sd :: (Nat -> Rational) -> C0 cauchy2sd = coitC0 step where step :: (Nat -> Rational) -> J0(Nat -> Rational) step f = (d,f') where q = f (Succ (Succ Zero)) d = if q > 1/4 then P else if abs q <= 1/4 then Z else N f' n = 2 * f (Succ n) - sd2Rational d -- 3. Memoized representation of u.c. functions -- ============================================ -- (see the end of this script for auxiliary definitions) data K1 alpha beta = W1 SD alpha | R1 (Triple beta) deriving Show -- functoriality of K1 mapK1 :: (alpha -> alpha') -> (beta -> beta') -> K1 alpha beta -> K1 alpha' beta' mapK1 f g (W1 d a) = W1 d (f a) mapK1 f g (R1 (bN,bZ,bP)) = R1 (g bN,g bZ,g bP) data J1 alpha = ConsJ1 (K1 alpha (J1 alpha)) deriving Show -- data -- iteration for J1 alpha itJ1 :: (K1 alpha beta -> beta) -> J1 alpha -> beta itJ1 s (ConsJ1 z) = s (mapK1 id (itJ1 s) z) -- functoriality of J1 mapJ1 :: (alpha -> alpha') -> J1 alpha -> J1 alpha' mapJ1 f (ConsJ1 x) = ConsJ1 (mapK1 f (mapJ1 f) x) data C1 = ConsC1 (J1 C1) deriving Show -- codata -- coiteration for C1 coitC1 :: (alpha -> J1 alpha) -> alpha -> C1 coitC1 s x = ConsC1 (mapJ1 (coitC1 s) (s x)) -- extracted programs (Proposition 1) -- application of C1 to C0 (n=1,m=0) appC :: C1 -> C0 -> C0 appC c ds = coitC0 costep (c,ds) where costep :: (C1,C0) -> J0 (C1,C0) costep (ConsC1 x,ds) = aux x ds aux :: J1 C1 -> C0 -> J0 (C1,C0) aux = itJ1 step step :: K1 C1 (C0 -> J0 (C1,C0)) -> C0 -> J0 (C1,C0) step (W1 d c') ds = (d,(c',ds)) step (R1 es) (ConsC0 (d0,ds')) = appTriple es d0 ds' -- composition on C1 (n=m=1) compC1 :: C1 -> C1 -> C1 compC1 c1 c2 = coitC1 costep (c1,c2) where costep :: (C1,C1) -> J1 (C1,C1) costep (ConsC1 x1,c2) = aux x1 c2 aux :: J1 C1 -> C1 -> J1(C1,C1) aux = itJ1 step step :: K1 C1 (C1 -> J1 (C1,C1)) -> C1 -> J1 (C1,C1) step (W1 d1 c1') c2 = ConsJ1 (W1 d1 (c1',c2)) step (R1 es) (ConsC1 x2) = subaux x2 where subaux :: J1 C1 -> J1 (C1,C1) subaux = itJ1 substep substep :: K1 C1 (J1 (C1,C1)) -> J1 (C1,C1) substep (W1 d2 c2') = appTriple es d2 c2' substep (R1 fs) = ConsJ1 (R1 fs) -- 4. Wellfounded induction/recursion and digital systems -- ====================================================== -- induction on a wellfounded relation < on A: -- if X is < progressive, -- i.e. all x in A (all y in A (y < x => X(y)) => X(x)), -- then A is a subset of X. -- If the definition of < has no computational content, -- then this is realised by wfrec :: (a -> (a -> b) -> b) -> a -> b wfrec prog = h where h x = prog x h -- digital systems (Proposition 2, n=1) digitsys1 :: (a -> Either (SD,a) (Triple a)) -> a -> C1 digitsys1 s = coitC1 (wfrec prog) where -- prog :: a -> (a -> J1 a) -> J1 a prog x ih = case s x of {Left (d,a') -> ConsJ1 (W1 d a') ; Right as -> ConsJ1 (R1 (mapTriple ih as))} -- 5. Example: linear affine maps f_(u,v)(x) = u*x+v -- ================================================= linC1 :: Rat2 -> C1 linC1 = digitsys1 s where s :: Rat2 -> Either (SD,Rat2) (Triple Rat2) s (u,v) = if u <= 1/4 then let w = if abs v <= 1/4 then 0 else signum v e = signumSD 0 w in Left (e,(2*u,2*v-w)) else Right (abstTriple (\d -> (u/2,u*fromSD d/2+v))) -- runC (linC (1/15,1/3)) (period [P,Z]) -- 6. Example: quadratic maps f_(u,v,w)(x) = u*x^2+v*x+w -- ===================================================== quadC1 :: Rat3 -> C1 quadC1 = digitsys1 s where s :: Rat3 -> Either (SD,Rat3) (Triple Rat3) s uvw = case (filter (quadTest uvw) [N,Z,P]) of (e:_) -> Left (e,quadWrite uvw e) [] -> Right (abstTriple (quadRead uvw)) quadWrite :: Rat3 -> SD -> Rat3 quadWrite (u,v,w) e = (2*u , 2*v , 2*w - e') where e' = fromSD e quadRead :: Rat3 -> SD -> Rat3 quadRead (u,v,w) d = (u/4 , (u*d'+v)/2 , u*d'^2/4 + v*d'/2 + w) where d' = fromSD d quadTest :: Rat3 -> SD -> Bool quadTest (u,v,w) e = (e'-1)/2 <= low && high <= (e'+1)/2 where e' = fromSD e low = minimum criticals -- max (f_(u,v,w) I) high = maximum criticals -- min (f_(u,v,w) I) criticals = [ u+v+w -- f_(u,v,w) 1 , u-v+w ] -- f_(u,v,w) (-1) ++ (if u == 0 then [] else let x = -v/(2*u) -- extremal point in if -1 <= x && x <= 1 then [ u*x^2 + v*x + w ] -- f_(u,v,w) x else []) -- 7. Example: logistic maps and their iterations -- =========================================== lma :: (Num a) => a -> a -> a lma a x = a*(1-x^2)-1 -- 0 <= a <= 2 -- = -a*x^2 + a - 1 lmaC1 :: Rational -> C1 lmaC1 a = quadC1 (-a,0,a-1) lm12C1,lm23C1,lm34C1,lm2C1 :: C1 lm12C1 = lmaC1 (1/2) lm23C1 = lmaC1 (2/3) lm34C1 = lmaC1 (3/4) lm2C1 = lmaC1 2 -- Iterating composition iterC1 :: C1 -> Int -> C1 iterC1 c n = iterate (compC1 c) c !! (n-1) ---------------------------------------------------------- -- 8. Auxiliary definitions -- ======================== type Triple alpha = (alpha,alpha,alpha) -- Triple a is isomorphic to SD -> a appTriple :: Triple alpha -> SD -> alpha appTriple (xN,xZ,xP) d = case d of {N -> xN ; Z -> xZ ; P -> xP} abstTriple :: (SD -> alpha) -> Triple alpha abstTriple f = (f N, f Z, f P) -- functoriality of Triple mapTriple :: (alpha -> beta) -> Triple alpha -> Triple beta mapTriple f (xN,xZ,xP) = (f xN,f xZ,f xP) -- every rational number in [-1,1] is in C0: ratC0 :: Rational -> C0 ratC0 = coitC0 step where step :: Rational -> J0 Rational step q = let e = signumSD (1/4) q in (e , 2 * fromSD e - q) signumSD :: (Num a,Ord a) => a -> a -> SD signumSD eps x | x < -eps = N | x > eps = P | otherwise = Z type Rat2 = (Rational,Rational) type Rat3 = (Rational,Rational,Rational) fromSD :: Num a => SD -> a fromSD d = case d of {N -> -1 ; Z -> 0 ; P -> 1} toSD :: (Num a,Ord a) => a -> SD toSD = signumSD 0 -- 9. Programs following strictly the scheme of the last part of Sect. 3 -- ===================================================================== type N2 = Int -- {1,2} data K2 alpha beta = W2 SD alpha | R2 N2 (Triple beta) deriving Show data J2 alpha = ConsJ2 (K2 alpha (J2 alpha)) -- data deriving Show data C2 = ConsC2 (J2 C2) -- codata deriving Show -------------------------------------------------------- -- In general: if -- Phi :: * -> * -- is a functor with action on morphism -- mapPhi :: (alpha -> beta) -> Phi alpha -> Phi beta -- both given here by dummy definitions data Phi a = PhiDef mapPhi :: (alpha -> beta) -> Phi alpha -> Phi beta mapPhi = undefined -- then data Fix = ConsFix (Phi Fix) coitFix :: (alpha -> Phi alpha) -> alpha -> Fix -- Fix viewed as codata coitFix s x = ConsFix (mapPhi (coitFix s) (s x)) itFix :: (Phi alpha -> alpha) -> Fix -> alpha -- Fix viewed as data itFix s (ConsFix z) = s (mapPhi (itFix s) z) -------------------------------------------------------- -- Programming the type Nat and the functions parity and parity1, -- following pedantically the general pattern above, gives: data PhiNat alpha = PhiZero | PhiSucc alpha deriving Show -- functoriality of PhiNat mapPhiNat :: (alpha -> beta) -> PhiNat alpha -> PhiNat beta mapPhiNat f PhiZero = PhiZero mapPhiNat f (PhiSucc x) = PhiSucc (f x) data FixNat = ConsFixNat (PhiNat FixNat) deriving Show -- data -- zero and successor for Fixnat zeroFixNat :: FixNat zeroFixNat = ConsFixNat PhiZero succFixNat :: FixNat -> FixNat succFixNat = ConsFixNat . PhiSucc -- iterator itFixNat :: (PhiNat alpha -> alpha) -> FixNat -> alpha itFixNat s (ConsFixNat z) = s (mapPhiNat (itFixNat s) z) -- extracted programs parity' :: FixNat -> Bool parity' = itFixNat step where step :: PhiNat Bool -> Bool step PhiZero = True step (PhiSucc b) = case b of {True -> False; False -> True} parity1' :: FixNat -> (FixNat,Bool) parity1' = itFixNat step1 where step1 :: PhiNat (FixNat,Bool) -> (FixNat,Bool) step1 PhiZero = (zeroFixNat,True) step1 (PhiSucc (y,b)) = case b of {True -> (y,False); False -> (succFixNat y,True)} -- Remark: It would be better programming style to define -- iterator and coiterator once and for all in a type class. -- We don't do this for the following reasons: -- 1. Bringing type classes into play would distract from the -- real issue of these demo programs. -- 2. The programs would become harder to understand for people less -- familiar with Haskell. -- 3. Type classes are good for hand crafted programs, but we want -- to extract programs mechanically, eventually. Hence high level -- concepts such as type classes are not needed (moreover, they -- make programs slower). -- mapJ1 f = itJ1 (ConsJ1 . mapK1 f id) -- beta = J1 alpha' -- mapK1 f id :: K1 alpha beta -> K1 alpha' beta -- ConsJ1 . mapK1 f id :: K1 alpha beta -> beta -- itJ1 (ConsJ1 . mapK1 f id) :: J1 alpha -> beta -------------------------------- -- latexing trees data Tree a = Node a [Tree a] -- codata -- Add positions to the nodes of a tree given with of tree and position of root posTree :: Tree a -> Float -> Point -> Tree (a,Point) posTree (Node x ts) width p = -- p = position of root if width <= thresholdWidth then Node (x,p) [] else let {l = length ts ; width' = width / fromIntegral l ; is = [0..l-1] ; f i = smooth2 (addPoint p (-width/2 + width'/2 + fromIntegral i*width', verticalOffset)) } in Node (x,p) [posTree (ts !! i) width' (f i) | i <- is] {- The goal is to create from a tree a tikz picture using the syntax of the example below: \begin{tikzpicture}[fill=blue!20] \path (0,0) node(a) [circle,draw,fill] {A} (3,-1) node(b) [circle,draw,fill] {B} (5,2) node(c) [circle,draw,fill] {C}; \draw[thick] (a) -- (b); \draw[thick] (c) -- (b); \end{tikzpicture} To use tikz add \usepackage{tikz} to the preamble of your latex document. -} -- Computing the vertices (paths) of a tree given -- a vertex for the root. vertices :: Tree a -> [Int] -> [(a,[Int])] vertices (Node x ts) path = -- path is the vertex attached to the root (x,path) : concat [vertices (ts!!i) (i:path) | i <- [0..length ts-1]] -- Computing the edges of a tree (as pairs of vertices) -- given a vertex for the root. edges :: Tree a -> [Int] -> [([Int],[Int])] edges (Node _ ts) path = [(path,i:path) | i <- [0..length ts-1]] ++ concat [edges (ts!!i) (i:path) | i <- [0..length ts-1]] -- Some global constants and utility types and functions. verticalOffset :: Float verticalOffset = -1 -- verticalOffset = -1.5 thresholdWidth :: Float thresholdWidth = 1.2 unitlength :: String unitlength = "1.2em" type Point = (Float,Float) addPoint :: Point -> Point -> Point addPoint (x,y) (x',y') = (x+x',y+y') round1dec :: Float -> Float round1dec x = fromIntegral (round (10 * x)) / 10 round2dec :: Float -> Float round2dec x = fromIntegral (round (100 * x)) / 100 smooth :: Point -> Point smooth (x,y) = (round1dec x, round1dec y) smooth2 :: Point -> Point smooth2 (x,y) = (round2dec x, round2dec y) -- Computing a tikz picture from a tree. treetikz :: Tree String -> Float -> Point -> String treetikz tree width p = let { pt = posTree tree width p ; vs = vertices pt [0] ; es = edges pt [0] ; beg = "\\begin{tikzpicture}[inner sep = 0.4mm,fill=blue!20]\n\\path" ; end = "\\end{tikzpicture}" ; showpath = map (toEnum . (fromEnum 'a' +)) ; nod ((label,q),path) = show q ++ " node(" ++ showpath path ++ ") " -- ++ "[circle,draw,fill] {" ++ label ++ "}" ; ++ "[circle,draw] {" ++ label ++ "}" ; edg (path1,path2) = "\\draw[thick] (" ++ showpath path1 ++ ") -- (" ++ showpath path2 ++ ");" } in beg ++ unlines (map nod vs) ++ ";\n" ++ unlines (map edg es) ++ end -- Creating a complete latex document for latex code using tikz mktikzdoc :: String -> String mktikzdoc str = "\\documentclass{article}\n" ++ "\\setlength{\\unitlength}{" ++ unitlength ++ "}\n" ++ "\\usepackage{tikz}\n" ++ "\\begin{document}\n" ++ str ++ "\n\\end{document}" treetest :: FilePath -> Float -> Point -> Rational -> IO () treetest fn width p q = writeFile fn (treetikz (treeC1 (lmaC1 q)) width p) tikztest1 :: FilePath -> Float -> Point -> Rational -> IO () tikztest1 fn width p q = writeFile fn (mktikzdoc (treetikz (treeC1 (lmaC1 q)) width p)) -- C1 into Tree treeC1 :: C1 -> Tree String treeC1 (ConsC1 x) = aux x where aux :: J1 C1 -> Tree String aux = itJ1 step step :: K1 C1 (Tree String) -> Tree String step (W1 d c) = Node (show d) [treeC1 c] step (R1 (tN,tZ,tP)) = Node "" [tN,tZ,tP]