\documentclass[leqno,fleqn,12pt]{article}%
\usepackage{MathSpad,makeidx,ifthen,latexsym}

\usepackage{a4wide, amssymb,latexsym,euler, beton,graphs}
\usepackage{concrete}
\usepackage[T1]{fontenc}


%include lhs2TeX.fmt
%include lhs2TeX.sty
%include polycode.fmt

\begin{document}
\setcounter{page}{13}

 
\appendix{}
\section{Appendix: Haskell Implementation}\label{Haskell Implementation}



This appendix contains an encoding of the enumeration algorithms in Haskell.  The file from which this
printed version was compiled  is a
so-called ``lhs2\TeX{}'' file\footnote{The lhs2\TeX{} system
 has been implemented by Ralf Hinze and Andres L\"oh.} which can be used directly as input to a Haskell 
compiler; this safeguards against
typographical errors in the printed paper.  

%format `multM` = "\times"
%format `multMrow` = "\times"
%format * = "\mskip-4mu\times\mskip-4mu"
%format floor(x) = "\lfloor" x "\rfloor"
%format % = "/"
%format == = "=="
%format $ = "\ "
\def\fraction#1#2{{}^{\smash{#1}}\mskip-4mu/\mskip-3mu_{\smash{#2}}}
\def\mat(#1,#2,#3,#4){\left({#1 \atop #3}\;{#2 \atop #4}\right)}
\def\vector(#1,#2){\left({#1\ms{2}#2}\right)}
\def\vec(#1,#2){\left({#1 \atop #2}\right)}



%if False

> module EnumeratingRationals where

> import Prelude hiding (gcd)
> import Data.List
> import Ratio

%endif

The implementation encodes a matrix as a list of \emph{columns}.   

> type Entry   = Integer
> type Column = [Entry]
> type Matrix  = [Column] 

%if False

> multM :: Matrix -> Matrix -> Matrix
> -- Pre-multiplying a 1x2 vector by a 2x2 matrix
> multM [[x],[y]] [[m,n],[m',n']] = [[ x*m + y*n ] , [ x*m' + y*n' ]]
> -- Post-multiplying a 2x1 vector by a 2x2 matrix
> multM [[m,n],[m',n']] [[x,y]] = [ [ m*x + m'*y , n*x + n'*y ] ]
> -- Multiplying a 2x2 matrix by a 2x2 matrix
> multM [[m,n],[m',n']] [[x,y],[z,w]] = [ [ m*x + m'*y  ,  n*x + n'*y ],
>                                         [ m*z + m'*w  ,  n*z + n'*w ] ]

> matA :: Matrix
> matA = [[1,-1],
>         [0,1]]
> matB :: Matrix
> matB = [[1,0],
>         [-1,1]]
> matId :: Matrix
> matId = [[1,0],
>          [0,1]]

> invA  :: Matrix
> invA  = [[1,1],
>          [0,1]]
> invB  :: Matrix
> invB  = [[1,0],
>          [1,1]]

> matL = invA
> matR = invB

%endif

We define a type of non-empty trees, with associated map, fold and unfold functions.

> data Tree a             = Node (a, Tree a, Tree a) 
> mapt f (Node (a,l,r))   = Node (f a, mapt f l, mapt f r)
> foldt f (Node (a,l,r))  = f (a, foldt f l, foldt f r)
> unfoldt f x             = let (a,y,z) = f x in Node (a, unfoldt f y, unfoldt f z)

%if False

> takeTree                   :: Integer -> Tree a -> Tree a
> takeTree 0 _               = error "There is no empty tree"
> takeTree n (Node (a,l,r))  = Node (a, takeTree (n-1) l, takeTree (n-1) r)

%endif

With matrices $\mathit{matId\/}$,  $\mathit{matL\/}$ and $\mathit{matR\/}$ defined to be the identity matrix, the matrix \setms{0.3em}$\mathbf{L\/}$ and the matrix $\mathbf{R\/}$,
respectively, the tree of matrices is generated as follows.

> mTree  :: Tree Matrix
> mTree  = unfoldt level matId
>  where level m = (m, m `multM` matL, m `multM` matR)


The Calkin-Wilf tree can be obtained by pre-multiplying the matrices by the
vector $\vector(1,1)$.


%if False

> cwTree'  :: Tree Rational
> cwTree'  = foldt rat mTree
>  where rat  (m,l,r)        = Node (mkCWRat ([[1], [1]] `multM` m), l, r)

%endif

> cwTree  :: Tree Rational
> cwTree  = mapt (mkCWRat . ([[1],[1]] `multM`)) mTree

> mkCWRat          :: Matrix -> Rational
> mkCWRat [[m], [n]]  =  n%m



Similarly, the 
Stern-Brocot tree can be obtained by post-multiplying the matrices by
the transpose of the vector $\vector(1,1)$ , i.e.,  $\vec(1,1)$.


%if False

> sbTree'  :: Tree Rational
> sbTree'  = foldt rat mTree
>  where rat  (m,l,r)        = Node (mkSBRat (m `multM` [[1,1]]), l, r)

%endif


> sbTree  :: Tree Rational
> sbTree  = mapt (mkSBRat . (`multM` [[1,1]])) mTree

> mkSBRat            :: Matrix -> Rational
> mkSBRat [[m,n]]  =  m%n



We enumerate the matrices using the |iterate| function, computing each
matrix from the previous one. 

> nextM :: Matrix -> Matrix
> nextM [[1, 0], [n, 1]]  =  [[1,n+1], [0,1]]  
> nextM [c0, c1]            =  let  j  = floor (((sum c1) -1)%(sum c0))
>                                   k  = 2*j+1
>                                   ck = map (k*) c0
>                              in   [zipWith (-) ck c1 , c0]


> mats  :: [Matrix]
> mats  = iterate nextM matId


The (non-optimised)  implementation of the Calkin-Wilf enumeration is then a matter of
premultiplying by $\vector(1,1)$.  

%if False

> cwEnum' :: [Rational]
> cwEnum' = map  mkCW mats
>  where mkCW = mkCWRat . (premul [[1],  [1]])
> premul v m = v `multM` m

%endif

> cwEnum :: [Rational]
> cwEnum = map  mkCW mats
>  where mkCW = mkCWRat . ([[1],  [1]] `multM`)


The Stern-Brocot enumeration can be defined in a similar way, but instead of
premultiplying, we postmultiply by  $\vec(1,1)$:

%if False

> sbEnum' :: [Rational]
> sbEnum' = map  mkSB mats
>  where mkSB = mkSBRat . (posmul [[1, 1]])
> posmul v m = m `multM` v

%endif


> sbEnum :: [Rational]
> sbEnum = map mkSB mats
>  where mkSB = mkSBRat . (`multM` [[1, 1]])


Incorporating the optimisations discussed above, the Calkin-Wilf
enumeration is transformed to the algorithm attributed to Newman.

> cwnEnum :: [Rational]
> cwnEnum = iterate nextCW $ 1%1

> nextCW :: Rational -> Rational
> nextCW r       =  let  (n,m)  = (numerator r, denominator r)
>                        j      = floor (n%m)
>                   in   m%((2*j+1)*m-n)


\end{document}