|| ***************************************************************************
|| Author       : Peter John Potts
|| File         : thesis.m
|| Description  : Exact Real Arithmetic using Mobius Transformations
|| Example      : "eshow ee 40" gives the natural number to
||                40 exact floating point digits of precision
||                converted to decimal, namely "0.27182818284e1".
|| ***************************************************************************

|| ***************************************************************************
|| Type Definitions
|| ***************************************************************************

|| Linear Fractional Transformations

vector == (num,num)
matrix == (vector,vector)
tensor == (matrix,matrix)
lft ::= LftV vector | LftM matrix | LftT tensor num

|| Expression Tree

expression ::=
    ExpV vector |
    ExpM matrix expression |
    ExpT tensor num expression expression

|| Exact Floating Point

sefp ::= Spos uefp | Sinf uefp | Sneg uefp | Szer uefp
digits == (num,num)
uefp == (digits,expression)

|| ***************************************************************************
|| Term Definitions
|| ***************************************************************************

|| Basic Functions

one :: * -> num -> *
one f 1 = f

identity :: matrix
identity = ((1,0),(0,1))

trans :: ((*,*),(*,*)) -> ((*,*),(*,*))
trans ((a,b),(c,d)) = ((a,c),(b,d))

determinant :: matrix -> num
determinant ((a,b),(c,d)) = a * d - b * c

inverse :: matrix -> matrix
inverse ((a,b),(c,d)) = ((d,-b),(-c,a))

|| Binary Scaling Functions

vscale :: vector -> vector
vscale (a,b)
  = vscale (a div 2,b div 2), a mod 2 = 0 & b mod 2 = 0
  = (a,b), otherwise

mscale :: matrix -> matrix
mscale ((a,b),(c,d))
  = mscale ((a div 2,b div 2),(c div 2,d div 2)),
        a mod 2 = 0 & b mod 2 = 0 & c mod 2 = 0 & d mod 2 = 0
  = ((a,b),(c,d)), otherwise

tscale :: tensor -> tensor
tscale (((a,b),(c,d)),((e,f),(g,h)))
  = tscale (((a div 2,b div 2),(c div 2,d div 2)),
        ((e div 2,f div 2),(g div 2,h div 2))),
        a mod 2 = 0 & b mod 2 = 0 & c mod 2 = 0 & d mod 2 = 0 &
        e mod 2 = 0 & f mod 2 = 0 & g mod 2 = 0 & h mod 2 = 0
  = (((a,b),(c,d)),((e,f),(g,h))), otherwise

|| Exact Floating Point

spos = ((1,0),(0,1))
sinf = ((1,-1),(1,1))
sneg = ((0,1),(-1,0))
szer = ((1,1),(-1,1))

ispos = LftM (inverse spos)
isinf = LftM (inverse sinf)
isneg = LftM (inverse sneg)
iszer = LftM (inverse szer)

dneg = ((1,1),(0,2))
dzer = ((3,1),(1,3))
dpos = ((2,0),(1,1))

idneg = LftM (inverse dneg)
idzer = LftM (inverse dzer)
idpos = LftM (inverse dpos)

|| Type Cast Functions

dtom :: digits -> matrix
dtom (n,c) = mscale ((2^n+c+1,2^n-c-1),(2^n+c-1,2^n-c+1))

utoe :: uefp -> expression
utoe (d,e) = ExpM (dtom d) e

stom :: sefp -> matrix
stom (Spos u) = mdotm spos (dtom (fst u))
stom (Sinf u) = mdotm sinf (dtom (fst u))
stom (Sneg u) = mdotm sneg (dtom (fst u))
stom (Szer u) = mdotm szer (dtom (fst u))

|| Basic Arithmetic Operations

tadd = (((0,0),(1,0)),((1,0),(0,1)))
tsub = (((0,0),(1,0)),((-1,0),(0,1)))
tmul = (((1,0),(0,0)),((0,0),(0,1)))
tdiv = (((0,0),(1,0)),((0,1),(0,0)))

srec :: sefp -> sefp
srec (Spos u) = Spos (urec u)
srec (Sneg u) = Sneg (urec u)
srec (Szer u) = Sinf (urec u)
srec (Sinf u) = Szer (urec u)

urec :: uefp -> uefp
urec ((n,c),e) = ((n,-c),erec e)

erec :: expression -> expression
erec e = ExpM ((0,1),(1,0)) e

|| Linear Fractional Transformation Products

mdotv :: matrix -> vector -> vector
mdotv ((a,b),(c,d)) (e,f) = (a * e + c * f,b * e + d * f)

mdotm :: matrix -> matrix -> matrix
mdotm m (v,w) = (mdotv m v,mdotv m w)

mdott :: matrix -> tensor -> tensor
mdott m (n,o) = (mdotm m n,mdotm m o)

tleftv :: tensor -> vector -> matrix
tleftv t v = trightv (trans t) v

tleftm :: tensor -> matrix -> tensor
tleftm t m = trans (trightm (trans t) m)
 
trightv :: tensor -> vector -> matrix
trightv (m,n) v = (mdotv m v,mdotv n v)
 
trightm :: tensor -> matrix -> tensor
trightm (m,n) o = (mdotm m o,mdotm n o)

dot :: num -> lft -> lft -> lft
dot 1 (LftM m) (LftV v) = LftV (vscale (mdotv m v))
dot 1 (LftM m) (LftM n) = LftM (mscale (mdotm m n))
dot 1 (LftM m) (LftT t i) = LftT (tscale (mdott m t)) i
dot 1 (LftT t i) (LftV v) = LftM (mscale (tleftv t v))
dot 1 (LftT t i) (LftM m)
  = LftT t i, m = identity
  = LftT (tscale (tleftm t m)) (i+1), otherwise
dot 2 (LftT t i) (LftV v) = LftM (mscale (trightv t v))
dot 2 (LftT t i) (LftM m)
  = LftT t i, m = identity
  = LftT (tscale (trightm t m)) (i+1), otherwise

|| The Refinement Property

sign :: vector -> num
sign (a,b)
  = -1, a<0 & b<=0
  =  0, a<0 & b> 0
  = -1, a=0 & b< 0
  =  0, a=0 & b= 0
  =  1, a=0 & b> 0
  =  0, a>0 & b< 0
  =  1, a>0 & b>=0

vrefine :: vector -> bool
vrefine v = sign v ~= 0

mrefine :: matrix -> bool
mrefine (v,w) = a = b & b ~= 0
                 where a = sign v
                       b = sign w

trefine :: tensor -> bool
trefine ((v,w),(x,y)) = a = b & b = c & c = d & d ~= 0
                           where a = sign v
                                 b = sign w
                                 c = sign x
                                 d = sign y

refine :: lft -> bool
refine (LftV v) = vrefine v
refine (LftM m) = mrefine m
refine (LftT t i) = trefine t

|| Square Bracket Application

app :: lft -> (num -> expression) -> expression
app (LftM m) g
  = cons (dot 1 (LftM m) (head (g 1))) (tail (g 1))
app (LftT t i) g
  = cons (dot 1 (dot 2 (LftT t i) (head (g 2))) (head (g 1))) h
    where c = branch (head (g 1))
          h i = tail (g 1) i, i <= c
              = tail (g 2) (i-c), otherwise

|| Tensor Absorption Strategy

vlessv :: vector -> vector -> bool
vlessv v w = determinant (v,w) < 0

mlessv :: matrix -> vector -> bool
mlessv (v,w) x = vlessv v x & vlessv w x

mlessm :: matrix -> matrix -> bool
mlessm m (v,w) = mlessv m v & mlessv m w

mdisjointm :: matrix -> matrix -> bool
mdisjointm m n = mlessm m n \/ mlessm n m

strategyf :: tensor -> num -> num
strategyf t i = i mod 2 + 1

strategyo :: tensor -> num -> num
strategyo t = strategyr t, trefine t
            = strategyf t, otherwise

strategyr :: tensor -> num -> num
strategyr t i = 2, mdisjointm (fst (trans t)) (snd (trans t))
              = 1, otherwise

decision :: num -> lft -> bool
decision 1 (LftM m) = True
decision 1 (LftT t i) = strategyo t i = 1
decision 2 (LftT t i) = strategyo t i = 2

|| Basic Expression Tree Functions

branch :: lft -> num
branch (LftV v) = 0
branch (LftM m) = 1
branch (LftT t i) = 2

vis :: lft -> bool
vis (LftV v) = True
vis e = False

mis :: lft -> bool
mis (LftM m) = True
mis e = False

tis :: lft -> bool
tis (LftT t i) = True
tis e = False

cons :: lft -> (num -> expression) -> expression
cons (LftV v) f = ExpV v
cons (LftM m) f = ExpM m (f 1)
cons (LftT t i) f = ExpT t i (f 1) (f 2)

head :: expression -> lft
head (ExpV v) = LftV v
head (ExpM m e) = LftM m
head (ExpT t i e f) = LftT t i

tail :: expression -> num -> expression
tail (ExpM m e) 1 = e
tail (ExpT t i e f) 1 = e
tail (ExpT t i e f) 2 = f

|| Normalization Functions

sem :: expression -> num -> sefp
sem e i
  = Spos (dem (0,0) (app ispos (one e)) i), refine (dot 1 ispos l)
  = Sneg (dem (0,0) (app isneg (one e)) i), refine (dot 1 isneg l)
  = Szer (dem (0,0) (app iszer (one e)) i), refine (dot 1 iszer l)
  = Sinf (dem (0,0) (app isinf (one e)) i), refine (dot 1 isinf l)
  = sem (app l f) i, otherwise
    where l = head e
          f d = ab l (tail e d) (decision d l)

dem :: digits -> expression -> num -> uefp
dem (i,c) e j
  = ((i,c),e), j = 0 \/ vis l
  = dem (i+1,2*c-1) (app idneg (one e)) (j-1), refine (dot 1 idneg l)
  = dem (i+1,2*c+1) (app idpos (one e)) (j-1), refine (dot 1 idpos l)
  = dem (i+1,2*c) (app idzer (one e)) (j-1), refine (dot 1 idzer l)
  = dem (i,c) (app l f) j, otherwise
    where l = head e
          f d = ab l (tail e d) (decision d l)

ab :: lft -> expression -> bool -> expression
ab k e b
  = utoe ((0,0),e), b = False
  = utoe (dem (0,0) e 1), tis k & tis (head e)
  = e, otherwise

|| Decimal Output Function

eshow :: expression -> num -> [char]
eshow e i = mshow (stom (sem e i))

mshow :: matrix -> [char]
mshow m
  = show p, d=0 & q=1
  = (show p) ++ "/" ++ (show q), d=0 & q~=1
  = sshow (scientific m 0), d~=0
    where d = determinant m
          (p,q) = vscale (fst m)

sshow :: [num] -> [char] 
sshow [] = "unbounded"
sshow (e : m)
  = (shows v) ++ (showm v) ++ (showe h)
    where f = (foldr g 0) . reverse
          g d c = d+10*c
          (h,l,v) = normalize (e,#m,f m)

normalize :: (num,num,num) -> (num,num,num) 
normalize (e,l,v)
  = normalize (e-1,l-1,v), l>0 & (abs v)<10^(l-1)
  = (e,l,v), otherwise

shows :: num -> [char] 
shows v = "-", v < 0
        = "", v >= 0

showm :: num -> [char] 
showm v = "0", v = 0
        = "0." ++ show (abs v), v ~= 0

showe :: num -> [char]
showe e = "e" ++ (show e)

scientific :: matrix -> num -> [num] 
scientific m n
  = [], vrefine (mdotv (inverse m) (1,0))
  = n : (mantissa (-9) 9 m), mrefine (mdotm (inverse szer) m)
  = scientific (mdotm ((1,0),(0,10)) m) (n+1), otherwise

mantissa :: num -> num -> matrix -> [num]
mantissa i n m
  = i : mantissa (-9) 9 (e i), c i
  = mantissa (i+1) n m, i < n
  = [], otherwise
    where c n = mrefine (mdotm (inverse (d n)) m)
          d n = ((n+1,10),(n-1,10))
          e n = mdotm ((10,0),(-n,1)) m

|| Elementary functions

eiterate :: (num -> matrix) -> num -> expression
eiterate i n = ExpM (i n) (eiterate i (n+1))

eiteratex :: (num -> tensor) ->num -> expression -> expression
eiteratex i n x = ExpT (i n) 0 x (eiteratex i (n+1) x)

esqrtrat :: num -> num -> expression
esqrtrat p q = rollover p q (p-q)

rollover :: num -> num -> num -> expression
rollover a b c = ExpM dneg (rollover (4*a) d c), d >= 0
               = ExpM dpos (rollover (-d) (4*b) c), otherwise
                 where d = 2*(b-a)+c

esqrtspos :: expression -> expression
esqrtspos = eiteratex itersqrtspos 0
itersqrtspos :: num -> tensor
itersqrtspos n = (((1,0),(2,1)),((1,2),(0,1)))

elogspos :: expression -> expression
elogspos = eiteratex iterlogspos 0
iterlogspos :: num -> tensor
iterlogspos 0 = (((1,0),(1,1)),((-1,1),(-1,0)))
iterlogspos n = (((n,0),(2*n+1,n+1)),((n+1,2*n+1),(0,n)))

ee :: expression
ee = eiterate itere 0
itere :: num -> matrix
itere n = ((2*n+2,2*n+1),(2*n+1,2*n))

eexpszer :: expression -> expression
eexpszer = eiteratex iterexpszer 0
iterexpszer :: num -> tensor
iterexpszer n = (((2*n+2,2*n+1),(2*n+1,2*n)),
                ((2*n,2*n+1),(2*n+1,2*n+2)))

epi :: expression
epi = ExpT tdiv 0 (esqrtrat 10005 1) eomega
eomega :: expression
eomega = eiterate iteromega 0
iteromega :: num -> matrix
iteromega 0 = ((6795705,213440),(6795704,213440))
iteromega n = ((e-d-c,e+d+c),(e+d-c,e-d+c))
              where b = (2*n-1)*(6*n-5)*(6*n-1)
                    c = b*(545140134*n + 13591409)
                    d = b*(n+1)
                    e = 10939058860032000*n^4

etanszer :: expression -> expression
etanszer = eiteratex itertanszer 0
itertanszer :: num -> tensor
itertanszer 0 = (((1,2),(1,0)),((-1,0),(-1,2)))
itertanszer n = (((2*n+1,2*n+3),(2*n-1,2*n+1)),
                ((2*n+1,2*n-1),(2*n+3,2*n+1)))

earctanszer :: expression -> expression
earctanszer = eiteratex iterarctanszer 0
iterarctanszer :: num -> tensor
iterarctanszer 0 = (((1,2),(1,0)),((-1,0),(-1,2)))
iterarctanszer n = (((2*n+1,n+1),(n,0)),((0,n),(n+1,2*n+1)))