spectral: remove compreals
authorMichal Terepeta <michal.terepeta@gmail.com>
Mon, 13 Mar 2017 22:37:14 +0000 (18:37 -0400)
committerBen Gamari <ben@smart-cactus.org>
Mon, 13 Mar 2017 23:44:06 +0000 (19:44 -0400)
It's been disabled for a long time and there's no input to actually
run it.

Signed-off-by: Michal Terepeta <michal.terepeta@gmail.com>
Test Plan: build & run

Reviewers: bgamari

Differential Revision: https://phabricator.haskell.org/D3330

spectral/Makefile
spectral/compreals/ContinuedFractions.lhs [deleted file]
spectral/compreals/Doc.lhs [deleted file]
spectral/compreals/Main.lhs [deleted file]
spectral/compreals/QRationals.lhs [deleted file]
spectral/compreals/RealReals.lhs [deleted file]
spectral/compreals/Transcendentals.lhs [deleted file]
spectral/compreals/paper.tex [deleted file]

index c66d3bd..64310af 100644 (file)
@@ -7,8 +7,7 @@ SUBDIRS = ansi atom awards banner boyer boyer2 calendar cichelli circsim \
          mandel mandel2 mate minimax multiplier para power pretty primetest \
          puzzle rewrite scc secretary simple sorting sphere treejoin
 
-# compreals    no suitable test data
-OTHER_SUBDIRS = compreals lambda last-piece triangle
+OTHER_SUBDIRS = lambda last-piece triangle
 
 include $(TOP)/mk/target.mk
 
diff --git a/spectral/compreals/ContinuedFractions.lhs b/spectral/compreals/ContinuedFractions.lhs
deleted file mode 100644 (file)
index 85d85d2..0000000
+++ /dev/null
@@ -1,458 +0,0 @@
-\chapter{Continued Fractions}
-
-\begin{verbatim}
-$Log: ContinuedFractions.lhs,v $
-Revision 1.2  1999/11/02 16:10:42  simonpj
-Haskell 98 changes
-
-Revision 1.1  1996/01/08 20:05:19  partain
-Initial revision
-
-\end{verbatim}
-
-> module ContinuedFractions
->  (module Maybe, module QRationals,
->   ContinuedFraction(..), Homography(..), Interval(..), rat2cf,
->   algebraicAlgorithm, algebraicOutput,
->   quadraticAlgorithm, quadraticOutput,
->   cfRat2CFSqrt, decimals, accuracy, cf2Rat, integerFraction)
-> where
-> import Data.Maybe
-> import Data.List( sort )
-> import QRationals
-
-\section{Representation of Continued Fractions}
-
-The {\em Continued Fraction} representation of a number $r$ is a list of numbers
-\[r = [x_0,~x_1,~x_2,~\ldots~x_i,~]\]
-with the following interpretation:
-\[r = x_0 + \frac{1}{x_1 + \frac{1}{x_2 +
-            \frac{1}{\ddots +\frac{1}{x_i + \frac{1}{ddots}}}}}\]
-
-In this section we will describe all of the representational issues
-associated with continued fractions.
-
-We represent continued fractions as a list of @Integer@. The
-@Homography@ is described seperately under the Algebraic and Quadratic
-Algorithms. An @Interval@ is simply an interval with rational
-end-points.
-
-> type ContinuedFraction = [Integer]
-> type Homography        = ([Integer],[Integer])
-> type Interval          = (QRational,QRational)
-
-To make it easy to negate a continued fraction, we are assuming that
-the continued fractions are symmetric in the following sense:
-\[-[x_0,\,x1,\,\ldots] = [ -x_0,\,-x_1,\,\ldots ] \]
-
-\subsection{Creating a (Finite) Continued Fraction from a QRational}
-
-The function @rat2cf@ turns a @QRational@ into the equivalent {\bf
-Z}-fraction. This is acceptable provided that the invariant that we
-choose for the relationship between the leading term of a continued
-fraction and its interval is such that:
-
-\[[\lfloor r\rceil - \frac{1}{2},\,\lfloor r\rceil +\frac{1}{2}) \subseteq
-@cfBound@\,r\, \]
-
-> rat2cf :: QRational -> ContinuedFraction
-> rat2cf q
->  = if qUndefined q then error "rat2cf: undefined" else
->    if qInfinite  q then []                        else
->    rat2cf' q
->    where rat2cf' q = x : if q' == 0 then [] else rat2cf' (1 / q')
->                      where x  = qRound q
->                            q' = q - (x%%1)
-
-\subsection{Enclosing Bounds for Continued Fractions}
-
-The function @cfBound@ gives @QRational@ bounds that enclose the
-Euclidean Continued Fraction with leading term @x@. Because @x@ is
-known to be @QRational@ -- and hence the denominator of the solution
-to the quadratic equation is non-zero -- both of the roots are finite.
-
-We observe the following conventions:
-\begin{itemize}
-\item The interval $(i,\,s)$ where $s<i$ is the interval
-$(-\infty,\,s)\cup(i,\,\infty)$.
-\end{itemize}
-
-> cfBound :: QRational -> Interval
-> cfBound x
->  = if qInfinite x then (1%%0,   1%%0)    else
->    if abs x < 2   then (-2%%1,  2%%1)    else
->    if x < 0       then (1%%2-x, x+1%%2)  else
->                        (x-1%%2, -x-1%%2)
-
-
- > cfBound x
- >  = if qInfinite x then (1%%0, 1%%0) else 
- >    if n > 0       then (c-pm,c+pm)  else
- >                        intersect (c-pm,c+pm)
- >    where x2 = x*x
- >          n  = 3 - x2 -- never 0 because x is QRational!
- >          c  = 4*x / n
- >          pm = r*(x2+1)/n
- >          r  = if x2 < 3 then rootThreePlus else rootThreeMinus
- >          ax = abs x - (1%%2)
- >          intersect (i,s)
- >            = if i >= s  then (max i ax, min s (-ax)) else
- >              if i < ax  then (ax,s) else
- >              if i > -ax then (i,-ax) else
- >                              error "cfBound: multiple intervals"
-
-Two guesses for $\sqrt 3$ are:
-
-> rootThreeMinus = 1732050807 %% 1000000000 :: QRational
-> rootThreePlus  = 1732050808 %% 1000000000 :: QRational
-
-These satisfy $@rootThreeMinus@ < \sqrt 3$ and $@rootThreePlus@ > \sqrt 3$.
-
-\subsection{Generating Output Terms of Continued Fractions}
-
-We need a way to generate the next term to be emmitted by the various
-algorithms that process the continued fractions; this is accomplished
-by the @generateOutput@ function.
-
-If the answer is now known with sufficient accuracy, the function
-returns @Just o@ where @o@ is the leading term of the output continued
-fraction.  If it is not possible to determine the leading term of the
-output continued fraction, @generateOutput@ returns @Nothing@.
-
-> generateOutput
->   :: (Homography -> [QRational] -> Interval)  -> -- Bound
->      (Homography -> [QRational] -> QRational) -> -- Evaluator
->      Homography                               -> -- Homography
->      [QRational]                              -> -- Leading Input Terms
->      Maybe Integer                               -- Possible Output
-
-> generateOutput bound eval homography qs
->  = if not (qInfinite q0 && qInfinite q1) --  && signum (o%%1) == signum o'
->       && acceptable int (o%%1) then Just o else Nothing
->    where int@(q0,q1) = bound homography qs
->          o' = eval homography qs
->          o  = qRound (if abs o' < 2 then o' else
->                       if abs q0 < abs q1 then q0 else q1)
-
-> acceptable :: Interval -> QRational -> Bool
-> acceptable int o = cfBound o `encloses` int
-
-\subsection{Interval Operations}
-
-The function @encloses@ tests whether its first interval argument
-encloses its second.
-
-> encloses :: Interval -> Interval -> Bool
-> (i0,s0) `encloses` (i1,s1)
->  = if i0 > s0 then 
->       if i1 < s1 then i0 <= i1 || s1 <= s0 
->                  else i0 <= i1 && s1 <= s0
->    else i1 < s1 && i0 <= i1 && s1 <= s0
-
-Using $A \cap B = (A^c \cup B^c)^c$ ...
-
-> intersect :: Interval -> Interval -> Interval
-> intersect (i,s) (i',s') = (rs,ri) where (ri,rs) = union [(s,i),(s',i')]
-
-> union :: [Interval] -> Interval
-> union = result . foldl join [] . sort . finitize
->         where result [int]         = int
->               result [(i,_),(_,s)] = (i,s)
->               result ints          = (0,0)
-
-Assuming that the intervals are finite and ordered....
-
-Then $@i@ \leq @i'@$ and if they're equal $@s@ \leq @s'@$
-
-> join :: [Interval] -> Interval -> [Interval]
-> join []             y         = [y]
-> join xs@((i,s):iss) y@(i',s') = if s >= s' then xs else
->                                 if s >= i' then (i,s'):iss else y:xs
-
-> finitize
->  = concat . map split
->    where split x@(i,s) = if i < s then [x] else [(-1%%0,s),(i,1%%0)]
-
-> finiteInterval   (i,s) = i < s
-> infiniteInterval (i,s) = i > s
-
-\section{The Algebraic Algorithm}
-
-The algebraic algorithm takes a homography and calculates the following:
-\[@algebraicAlgorithm@\,\left(\begin{array}[c]{cc} 
-n_0 & n_1 \\ d_0 & d_1\end{array}\right)\,x \, = \,
- \frac{n_0x+n_1}{d_0x+d1}\]
-
-Given that $x$ isn't $\bot$, this is well-defined provided that the
-determinant of the homography is non-zero. If the determinant {\em is}
-zero then the result is independent of the value of $x$.
-
-Because we notice that absorbing a (well-defined and finite) term from
-the input continued fraction or emitting a (well-defined and finite)
-term of the output continued fraction doesn't change the determinant
-of the homography we need only check this once.
-
-> algebraicAlgorithm :: Homography        -> -- Homography
->                       ContinuedFraction -> -- Input continued fraction
->                       ContinuedFraction    -- Output continued fraction
-
-> algebraicAlgorithm homography@([n0,n1],[d0,d1]) xs
->  = if n0*d1 /= n1*d0 then algebraicAlgorithmNoCheck homography xs
->    else if d1 /= 0   then rat2cf (n1%%d1)
->    else if d0 /= 0   then rat2cf (n0%%d0) 
->    else                   error "algebraicAlgorithm: undefined homography"
-
-> algebraicAlgorithmNoCheck :: Homography        -> -- Homography
->                              ContinuedFraction -> -- Input  CF
->                              ContinuedFraction    -- Output CF
-
-> algebraicAlgorithmNoCheck homography@([n0,n1],[d0,d1]) []
->  = rat2cf (n0%%d0)
-> algebraicAlgorithmNoCheck homography@([n0,n1],[d0,d1]) xs'@(x:xs)
->  = case algebraicOutput homography (x%%1) of
->      Just o  -> o : algebraicAlgorithmNoCheck ([d0,d1],[n0-d0*o,n1-d1*o]) xs'
->      Nothing ->     algebraicAlgorithmNoCheck ([n0*x+n1,n0],[d0*x+d1,d0]) xs
-
-> algebraicOutput :: Homography ->    -- Homography
->                    QRational  ->    -- Leading Term of Input CF
->                    Maybe Integer    -- Leading Term of Output CF 
-
-> algebraicOutput homography q
->  = if q == 0 then Nothing else a
->    where a = generateOutput algebraicBound algebraicEval homography [q]
-
-\subsection{Bounds on the Output Continued Fraction}
-
-The function @algebraicBound@, given a homography and the first
-(@QRational@) term of a continued fraction determines the interval
-within which the transformed continued fraction {\em must} lie.
-
-> algebraicBound :: Homography  -> -- Homography
->                   [QRational] -> -- Leading term of input CF
->                   Interval       -- Interval
-
-> algebraicBound homography@([n0,n1],[d0,d1]) [q]
->  = if n0*d1 == n1*d0
->    then (indp, indp) -- bound independent of x
->    else if xi == dcross || xs == dcross -- infinity at endpoint
->         then if xi > xs && (ncross <= xs || xi <= ncross) ||
->                 xi < xs &&  ncross <= xs && xi <= ncross
->              then zint  -- zero in interval
->              else nzint -- zero not in interval
->    else if xi > xs && (d0 == 0 || xi < dcross || dcross < xs) ||
->            xi < xs &&  d0 /= 0 && xi < dcross && dcross < xs
->         then (ys,yi) -- infinity in interval
->         else (yi,ys) -- infinity not in interval
->    where dcross  = (-d1) %% d0
->          ncross  = (-n1) %% n0
->          (xi,xs) = cfBound q
->          yi'     = algebraicEval homography [xi]
->          ys'     = algebraicEval homography [xs]
->          (yi,ys) = (min yi' ys', max yi' ys')
->          indp    = if d0 /= 0 then n0%%d0 else
->                    if d1 /= 0 then n1%%d1 else
->                       error "algebraicBound: denominator 0"
->          zint    = if qInfinite yi then if ys < 0%%1 then (ys,  1%%0)
->                                                      else (-1%%0, ys)
->                                    else if yi < 0%%1 then (yi,  1%%0)
->                                                      else (-1%%0, yi)
->          nzint   = if qInfinite yi then if ys < 0%%1 then (-1%%0, ys)
->                                                      else (ys,  1%%0)
->                                    else if yi < 0%%1 then (-1%%0, yi)
->                                                      else (yi,  1%%0)
-
-\subsection{Evaluating the homography}
-
-The function @algebraicEval@ evaluates the homography applied to a
-QRational number.
-
-> algebraicEval :: Homography  -> -- Homography
->                  [QRational] -> -- Input QRational
->                  QRational      -- Output QRational
-
-> algebraicEval ([n0,n1],[d0,d1]) [x] = (n0*xn + n1*xd) %% (d0*xn + d1*xd)
->                                        where xn = qNumerator   x
->                                              xd = qDenominator x
-
-\section{Quadratic Functions}
-
-In a similar fashion we can define the Quadratic Algorithm. This
-calculates the following sum:
-
-\[@quadraticAlgorithm@~
-\left(\begin{array}[c]{cccc}n_0 & n_1 & n_2 & n_3\\d_0 & d_1 & d_2 & d_3
-\end{array}\right)~x~y
-= \frac{n_0xy + n_1y + n_2x + n_3}{d_0xy +d_1y +d_2x + d_3}
-\]
-
-I could believe that the @quadraticAlgorithm@ might not work correctly
-for finite continued fractions, as I haven't checked the conditions on
-the homography for definedness. On the positive side, there shouldn't
-be any finite continued fractions anyway.
-
-> quadraticAlgorithm :: Homography         -> -- Homography
->                       ContinuedFraction  -> -- First Input CF
->                       ContinuedFraction  -> -- Second Input CF
->                       ContinuedFraction     -- Output CF
-
-> quadraticAlgorithm homography@([n0,n1,n2,n3],[d0,d1,d2,d3]) []  []
->  = rat2cf (n0 %% d0)
-> quadraticAlgorithm homography@([n0,n1,n2,n3],[d0,d1,d2,d3]) []  ys'@(y:ys)
->  = if aaDet homography' == 0
->    then quadraticAlgorithm (qaInRight homography y) [] ys
->    else algebraicAlgorithmNoCheck homography' ys'
->    where homography' = ([n0,n2],[d0,d2])
-> quadraticAlgorithm homography@([n0,n1,n2,n3],[d0,d1,d2,d3]) xs'@(x:xs) []
->  = if aaDet homography' == 0
->    then quadraticAlgorithm (qaInLeft homography x) xs []
->    else algebraicAlgorithmNoCheck homography' xs'
->    where homography' = ([n0,n1],[d0,d1])
-> quadraticAlgorithm homography xs'@(x:xs) ys'@(y:ys)
->  = case quadraticOutput homography [x%%1,y%%1] of
->      Just o  -> o : quadraticAlgorithm (qaOut homography o)   xs' ys'
->      Nothing ->     quadraticAlgorithm (qaIn  homography x y) xs  ys
-
-> quadraticOutput :: Homography   -> -- Homography
->                    [QRational]  -> -- Leading Term of Input CF
->                    Maybe Integer   -- Leading Term of Output CF 
-> quadraticOutput = generateOutput quadraticBound quadraticEval
-
-> qaIn :: Homography -> Integer -> Integer -> Homography
-> qaIn ([n0,n1,n2,n3],[d0,d1,d2,d3]) x y
->  = ([(n0*x+n1)*y+n2*x+n3, n0*y+n2, n0*x+n1, n0],
->     [(d0*x+d1)*y+d2*x+d3, d0*y+d2, d0*x+d1, d0])
-
-> qaOut :: Homography -> Integer -> Homography
-> qaOut ([n0,n1,n2,n3],[d0,d1,d2,d3]) o
->  = ([d0,d1,d2,d3], [n0-d0*o, n1-d1*o, n2-d2*o, n3-d3*o])
-
-> qaInLeft :: Homography -> Integer -> Homography
-> qaInLeft ([n0,n1,n2,n3],[d0,d1,d2,d3]) x
->  = ([n0*x+n1, n0, n2*x+n3, n2], [d0*x+d1, d0, d2*x+d3, d2])
-
-> qaInRight :: Homography -> Integer -> Homography
-> qaInRight ([n0,n1,n2,n3],[d0,d1,d2,d3]) y
->  = ([n0*y+n2, n1*y+n3, n0, n1], [d0*y+d2, d1*y+d3, d0, d1])
-
-> aaDet :: Homography -> Integer
-> aaDet ([n0,n1],[d0,d1]) = n0*d1 - n1*d0
-
-[Worry that it would be quicker to do Integer arithmetic in qaEval.]
-
-> quadraticEval :: Homography       -> -- homography
->                  [QRational]      -> -- QRational
->                  QRational                -- result
-
-> quadraticEval (ns,ds) [x, y]
->  = ((n0*x+n1)*y + (n2*x+n3)) / ((d0*x+d1)*y + (d2*x+d3))
->    where [n0,n1,n2,n3] = map (%% 1) ns
->          [d0,d1,d2,d3] = map (%% 1) ds
-
-> quadraticBound :: Homography  -> -- homography
->                   [QRational] -> -- terms of CFs
->                   Interval       -- interval
-
-> quadraticBound ([n0,n1,n2,n3],[d0,d1,d2,d3]) [x, y]
->  = union ints
->    where (xi,xs) = cfBound x
->          (yi,ys) = cfBound y
->          ints    = [algebraicBound (aax x) [y] | x <- [xi,xs]] ++
->                    [algebraicBound (aay y) [x] | y <- [yi,ys]]
->          aax x   = ([n0*xn+n1*xd,n2*xn+n3*xd],[d0*xn+d1*xd,d2*xn+d3*xd])
->                    where xn = qNumerator x
->                          xd = qDenominator x
->          aay y   = ([n0*yn+n2*yd,n1*yn+n3*yd],[d0*yn+d2*yd,d1*yn+d3*yd])
->                    where yn = qNumerator y
->                          yd = qDenominator y
-
-\section{Square Roots}
-
-\subsection{Rational Square Roots}
-
-Notice that there is no way we can test for rationalness of these
-continued fractions.
-
-> cfRat2CFSqrt :: QRational -> ContinuedFraction
-> cfRat2CFSqrt q
->  = map snd (next ([2*d*z,d,n-d*z*z,0], z)) 
->    where n = qNumerator q
->          d = qDenominator q
->          s = intSquareRoot (n*d)
->          z = qRound (s%%d)
->          t = s - z*d
->          next :: ([Integer],Integer) -> [([Integer],Integer)]
->          next (s,z) = if d0 == 0
->                       then [(s,z)]
->                       else (s,z): next ([n0',n1',d0',d1'],z')
->                       where [n0,n1,d0,d1] = s
->                             n0' = d0*z'+d1
->                             n1' = d0
->                             d0' = n0*z'+n1-z'*n0'
->                             d1' = n0-z'*d0
->                             z'  = qRound ((n0+t) %% d0)
-
-\subsection{Integer Square Roots}
-
-The function @intSquareRoot@ returns a non-negative integer $n$, such
-that $n = \lfloor\sqrt{@x@}\, \rceil$.
-
-> intSquareRoot :: Integer -> Integer
-> intSquareRoot x = until satisfy (improve x) x
->                   where satisfy y    = y2+y >= x && y2-y < x
->                                        where y2 = y*y
->                         improve x y  = (y*y+x) `ddiv` (2*y)
->                         ddiv x y     = if (r <= y `div` 2) then q else q+1
->                                        where (q,r) = divMod x y
-
-> root2 = cfRat2CFSqrt (2%%1)
-> root3 = cfRat2CFSqrt (3%%1)
-> test0 :: ContinuedFraction
-> test0 = f 1 where f n = 0: 10^n: f (n+1)
-
-For debugging.
-
-\section{Conversion of Continued Fractions to QRationals}
-
-I've chosen to print out 10 decimal places; you can add more (or less)
-by changing @decimals@.
-
-> decimals = 10 :: Int
-
-> accuracy = 1 %% (10^(toInteger decimals)) :: QRational
-
-To convert a continued fraction to a rational within the given
-accuracy, we use a variant of the Algebraic Algorithm.
-
-> cf2Rat :: ContinuedFraction -> QRational
-> cf2Rat = outputAlgorithm ([1,0],[0,1])
-
-> outputAlgorithm :: Homography        -> -- Homography
->                    ContinuedFraction -> -- Input CF
->                    QRational            -- Final Approximation
-
-> outputAlgorithm homography@([n0,n1],[d0,d1]) [] = n0 %% d0
-> outputAlgorithm homography@([n0,n1],[d0,d1]) xs'@(x:xs)
->  = if i <= s && s-i < accuracy
->    then algebraicEval homography [x%%1]
->    else outputAlgorithm ([n0*x+n1,n0],[d0*x+d1,d0]) xs
->    where (i,s) = algebraicBound homography [x%%1]
-
-
-> integerFraction :: ContinuedFraction           -> -- Input CF
->                    (Integer,ContinuedFraction)    -- Answer
-> integerFraction x = (o, algebraicAlgorithm h' f)
->                     where (h,f) = intFrac (([1,0],[0,1]),x)
->                           o     = qRound (algebraicEval h [head f%%1])
->                           h'    = ([n0-d0*o,n1-d1*o],[d0,d1])
->                                   where ([n0,n1],[d0,d1]) = h
-
-> intFrac :: (Homography,ContinuedFraction) -> (Homography,ContinuedFraction)
-> intFrac st@(_,[_]) = st
-> intFrac st@(homography@([n0,n1],[d0,d1]),(x:xs))
->  = if i <= s && s-i < 1 %% 10
->    then st
->    else intFrac (([n0*x+n1,n0],[d0*x+d1,d0]), xs)
->    where (i,s) = algebraicBound homography [x%%1]
-
-
-
diff --git a/spectral/compreals/Doc.lhs b/spectral/compreals/Doc.lhs
deleted file mode 100644 (file)
index cbb8a53..0000000
+++ /dev/null
@@ -1,13 +0,0 @@
-\documentstyle[11pt,a4]{report}
-%\includeonly{ContinuedFractions}
-
-\begin{document}
-\title{Perfect Arithmetic for the Irrational}
-\author{David Lester\\Department of Computer Science\\Manchester University, UK}
-\maketitle
-\include{Main}
-\include{RealReals}
-\include{Transcendentals}
-\include{ContinuedFractions}
-\include{QRationals}
-\end{document}
diff --git a/spectral/compreals/Main.lhs b/spectral/compreals/Main.lhs
deleted file mode 100644 (file)
index 0c13271..0000000
+++ /dev/null
@@ -1,100 +0,0 @@
-\chapter{Command Interpreter}
-
-\begin{verbatim}
-$Log: Main.lhs,v $
-Revision 1.3  1999/11/02 16:10:42  simonpj
-Haskell 98 changes
-
-Revision 1.2  1996/07/25 21:30:47  partain
-Bulk of final changes for 2.01
-
-Revision 1.1  1996/01/08 20:05:20  partain
-Initial revision
-
-\end{verbatim}
-
-Eventually we will build a command interpreter here, for a desk
-calculator type language.
-
-> module Main where
-> import RealReals
-> import Data.Char
-
-For the moment on the other hand it just prints a given number.
-
-> main = getContents >>= (foldr output (return ()) . map doLine . lines)
-
-Printing out @String@'s is easy:
-
-> output :: String -> IO () -> IO ()
-> output string dialogue = putStr (string++"\n") >> dialogue
-
-The @doLine@ function parses the line, if it is syntactically correct
-it then evaluates the expression returning the answer string.
-
-> doLine :: String -> String
-> doLine = eval [] . tokenize
-
-> tokenize ""                 = []
-> tokenize (c:cs) | isSpace c = tokenize (dropWhile isSpace cs)
-> tokenize (c:cs) | isSymb  c = [c]: tokenize cs
-> tokenize (c:cs) | isAlpha c = case (span isAlphaNum cs) of
->                                (nam,t) -> (c:nam): tokenize t
-> tokenize (c:cs) | isDigit c = case (span isDigit cs)  of
->                                (num,t) -> (c:num): tokenize t
-> tokenize _                  = ["Error"]
-
-> isSymb c = c `elem` "*+-/"
-
-> eval :: [RealReal] -> [String] -> String
-> eval [n] []     = show n
-> eval ns (t:ts) | isSymb  (head t)
->  = case head t of
->     '+' -> check2 (+) ns ts
->     '-' -> check2 (-) ns ts
->     '*' -> check2 (*) ns ts
->     '/' -> check2 (/) ns ts
->     _   -> "Error"
-> eval ns (t:ts) | isDigit (head t)
->  = eval (fromInteger (parseInteger t): ns) ts
-> eval ns (t:ts) | isAlpha (head t)
->  = case t of
->      "abs"    -> check1 abs    ns ts
->      "signum" -> check1 signum ns ts
->      "pi"     -> eval  (pi:ns)    ts
->      "exp"    -> check1 exp    ns ts
->      "log"    -> check1 log    ns ts
->      "sqrt"   -> check1 sqrt   ns ts
->      "sin"    -> check1 sin    ns ts
->      "cos"    -> check1 cos    ns ts
->      "tan"    -> check1 tan    ns ts
->      "asin"   -> check1 asin   ns ts
->      "acos"   -> check1 acos   ns ts
->      "atan"   -> check1 atan   ns ts
->      "sinh"   -> check1 sinh   ns ts
->      "cosh"   -> check1 cosh   ns ts
->      "tanh"   -> check1 tanh   ns ts
->      "asinh"  -> check1 asinh  ns ts
->      "acosh"  -> check1 acosh  ns ts
->      "atanh"  -> check1 atanh  ns ts
->      _        -> "Error"
-> eval _  _ = "Error"
-
-> check1 :: (RealReal -> RealReal) -> [RealReal] -> [String] -> String
-> check1 f (n:ns) ts = eval (f n: ns) ts
-> check1 f _      ts = "Error"
-
-> check2 :: (RealReal -> RealReal -> RealReal) ->
->           [RealReal] -> [String] -> String
-> check2 f (n0:n1:ns) ts = eval (f n1 n0 : ns) ts
-> check2 f _          ts = "Error"
-
-> parseInteger :: String -> Integer
-> parseInteger = makeNumber 10 . map number
->                where number :: Char -> Integer
->                      number c = toInteger (fromEnum c - fromEnum '0')
-
-> makeNumber :: Integer -> [Integer] -> Integer
-> makeNumber m = foldl f 0
->                where f a x = a * m + x
-
diff --git a/spectral/compreals/QRationals.lhs b/spectral/compreals/QRationals.lhs
deleted file mode 100644 (file)
index a065ed9..0000000
+++ /dev/null
@@ -1,95 +0,0 @@
-\chapter{Rational Numbers (with Infinity)}
-
-\begin{verbatim}
-$Log: QRationals.lhs,v $
-Revision 1.3  1999/11/02 16:10:42  simonpj
-Haskell 98 changes
-
-Revision 1.2  1996/07/25 21:30:48  partain
-Bulk of final changes for 2.01
-
-Revision 1.1  1996/01/08 20:05:18  partain
-Initial revision
-
-\end{verbatim}
-
-> module QRationals
->  (QRational, (%%), qNumerator, qDenominator,
->   qInfinite, qUndefined, qRound, qFinite) where
-> import Ratio
-> infixl 7  %% , :%%
-> prec = 7::Int
-
-> data  QRational              =  Integer :%% Integer deriving (Eq)
-
-> qReduce, (%%)                        :: Integer -> Integer -> QRational
-> qNumerator, qDenominator     :: QRational -> Integer
-
-Dealing with the following two functions' cases for $0$ is the only
-change to the standard prelude's rational package.
-
-> qReduce x 0                  =  signum x :%% 0
-> qReduce x y                  =  (signum x * (abs x `div` d)) :%% (y `div` d)
->                                 where d = gcd x y
-
-> x %% 0                       =  signum x :%% 0
-> x %% y                       =  qReduce (x * signum y) (abs y)
-
-> qNumerator (x:%%y)           =  x
-
-> qDenominator (x:%%y)         =  y
-
-> instance Ord QRational where
->     (x:%%y) <= (x':%%y')     =  x * y' <= x' * y
->     (x:%%y) <  (x':%%y')     =  x * y' <  x' * y
-
-> instance Num QRational where
->     (x:%%y) + (x':%%y')      =  qReduce (x*y' + x'*y) (y*y')
->     (x:%%y) - (x':%%y')      =  qReduce (x*y' - x'*y) (y*y')
->     (x:%%y) * (x':%%y')      =  qReduce (x * x') (y * y')
->     negate (x:%% y)          =  (-x) :%% y
->     abs (x:%% y)             =  abs x :%% y
->     signum (x:%% y)          =  signum x :%% 1
->     fromInteger x            =  fromInteger x :%% 1
-
-> instance Real QRational where
->     toRational (x:%%y)       =  x%y
-
-> instance Fractional QRational  where
->     (x:%%y) / (x':%%y')      =  (x*y') %% (y*x')
->     recip (x:%%y)            =  if x < 0 then (-y):%% (-x) else y:%%x
->     fromRational x           =  numerator x :%% denominator x
-
-> instance  Enum QRational where
->      enumFrom x              = enumFromBy x 1
->      enumFromThen x y        = enumFromBy x (y - x)
-
-> enumFromBy n k               =  n : enumFromBy (n+k) k
-
-> instance  Read QRational  where
->      readsPrec p             =  readParen (p > prec)
->                                   (\r -> [(x%%y,u) | (x,s)    <- reads r,
->                                                      ("%%",t) <- lex s,
->                                                      (y,u)    <- reads t ])
-> instance  Show QRational  where
->      showsPrec p (x:%%y)     = showParen (p > prec) (shows x .
->                                                       showString " %% " .
->                                                       shows y)
-
-\section{General extras for Rationals}
-
-> qInfinite   :: QRational -> Bool
-> qInfinite x  = qDenominator x == 0 && qNumerator x /= 0
-
-> qUndefined  :: QRational -> Bool
-> qUndefined x = qDenominator x == 0 && qNumerator x == 0
-
-> qRound      :: QRational -> Integer
-> qRound x     = if qUndefined x || qInfinite x
->                then error "qRound: undefined or infinite"
->                else (2*qNumerator x + d) `div` (2*d) where d = qDenominator x
-
-> qFinite     :: QRational -> Bool
-> qFinite x    = qDenominator x /= 0
-
-This concludes the slightly extended rational package.
diff --git a/spectral/compreals/RealReals.lhs b/spectral/compreals/RealReals.lhs
deleted file mode 100644 (file)
index 22d355c..0000000
+++ /dev/null
@@ -1,373 +0,0 @@
-\chapter{Implementation of @RealReal@'s}
-
-\begin{verbatim}
-$Log: RealReals.lhs,v $
-Revision 1.3  1999/11/02 16:10:42  simonpj
-Haskell 98 changes
-
-Revision 1.2  1996/07/25 21:30:49  partain
-Bulk of final changes for 2.01
-
-Revision 1.1  1996/01/08 20:05:19  partain
-Initial revision
-
-\end{verbatim}
-
-> module RealReals (RealReal) where
-> import Transcendentals
-> import Data.Ratio
-> import Data.List( genericLength )
-> import Numeric( readSigned )
-
-> data RealReal = RealInt Integer           |
->                 RealRat QRational         |
->                 RealCF  ContinuedFraction
-
-\section{@RealReal@ is an instance of @Eq@}
-
-It should be remembered that computable real numbers do not permit 
-computable comparisons, but never mind!
-
-> instance Eq RealReal where
->   (==) = realRealComp (==)
->   (/=) = realRealComp (/=)
-
-The actual comparison is carried out by subtracting @y@ from @x@ and
-testing the result at our current accuracy.
-
-> realRealComp :: (QRational -> QRational -> Bool) -> -- Operator
->                 RealReal                         -> -- First Argument
->                 RealReal                         -> -- Second Argument
->                 Bool                                -- Result
-> realRealComp op x y = op (makeRational (x-y)) 0
-
-\section{@RealReal@ is an instance of @Ord@}
-
-Once more, these operations are not computable.
-
-> instance Ord RealReal where
->   (<=) = realRealComp (<=)
->   (<)  = realRealComp (<)
->   (>=) = realRealComp (>=)
->   (>)  = realRealComp (>)
-
-\section{@RealReal@ is an instance of @Enum@}
-
-We make @RealReal@ an instance of @Enum@.
-
-> instance Enum RealReal where
->   enumFrom         = iterate (+1)
->   enumFromThen n m = iterate (+(m-n)) n
-
-\section{@RealReal@ is an instance of @Num@}
-
-In this section we come to operations that {\em are} computable.
-
-> instance Num RealReal where
->   x + y         = realRealAdd x y
->   negate x      = realRealNegate x
->   x * y         = realRealMultiply x y
->   abs x         = realRealAbs x
->   signum x      = realRealSignum x
->   fromInteger x = RealInt x
-
-\subsection{Addition}
-
-> realRealAdd :: RealReal -> RealReal -> RealReal
-> realRealAdd (RealInt x) (RealInt y) = RealInt (x+y)
-> realRealAdd (RealInt x) (RealRat y) = wrapRat ((x%%1)+y)
-> realRealAdd (RealInt x) (RealCF  y) = RealCF (algebraicAlgorithm h y)
->                                       where h  = ([1,x],[0,1])
-> realRealAdd (RealRat x) (RealInt y) = wrapRat (x+(y%%1))
-> realRealAdd (RealRat x) (RealRat y) = wrapRat (x+y)
-> realRealAdd (RealRat x) (RealCF  y) = RealCF (algebraicAlgorithm h y)
->                                       where nx = qNumerator   x
->                                             dx = qDenominator x
->                                             h  = ([dx,nx],[0,dx])
-> realRealAdd (RealCF  x) (RealInt y) = RealCF (algebraicAlgorithm h x)
->                                       where h  = ([1,y],[0,1])
-> realRealAdd (RealCF  x) (RealRat y) = RealCF (algebraicAlgorithm h x)
->                                       where ny = qNumerator   y
->                                             dy = qDenominator y
->                                             h  = ([dy,ny],[0,dy])
-> realRealAdd (RealCF  x) (RealCF  y) = RealCF (quadraticAlgorithm h x y)
->                                       where h  = ([0,1,1,0],[0,0,0,1])
-
-The @wrapRat@ function ensures that a rational number that is in fact
-an integer, is so treated.
-
-> wrapRat :: QRational -> RealReal
-> wrapRat x = if qDenominator x == 1 then RealInt (qNumerator x) else
->                                         RealRat x
-
-\subsection{Negation}
-
-Here we are using the symmetry of the continued fraction
-representation, so that we can simply map @negate@ over the continued
-fraction.
-
-> realRealNegate :: RealReal -> RealReal
-> realRealNegate (RealInt x) = RealInt (negate x)
-> realRealNegate (RealRat x) = RealRat (negate x)
-> realRealNegate (RealCF  x) = RealCF  (map negate x)
-
-\subsection{Multiplication}
-
-> realRealMultiply :: RealReal -> RealReal -> RealReal
-> realRealMultiply (RealInt x) (RealInt y) = RealInt (x*y)
-> realRealMultiply (RealInt x) (RealRat y) = wrapRat ((x%%1)*y)
-> realRealMultiply (RealInt x) (RealCF  y) = RealCF (algebraicAlgorithm h y)
->                                            where h  = ([x,0],[0,1])
-> realRealMultiply (RealRat x) (RealInt y) = wrapRat (x*(y%%1))
-> realRealMultiply (RealRat x) (RealRat y) = wrapRat (x*y)
-> realRealMultiply (RealRat x) (RealCF  y) = RealCF (algebraicAlgorithm h y)
->                                            where nx = qNumerator   x
->                                                  dx = qDenominator x
->                                                  h  = ([nx,0],[0,dx])
-> realRealMultiply (RealCF  x) (RealInt y) = RealCF (algebraicAlgorithm h x)
->                                            where h  = ([y,0],[0,1])
-> realRealMultiply (RealCF  x) (RealRat y) = RealCF (algebraicAlgorithm h x)
->                                            where ny = qNumerator   y
->                                                  dy = qDenominator y
->                                                  h  = ([ny,0],[0,dy])
-> realRealMultiply (RealCF  x) (RealCF  y) = RealCF (quadraticAlgorithm h x y)
->                                            where h  = ([1,0,0,0],[0,0,0,1])
-
-\subsection{Abs}
-
-We have a problem here. If @signum x == -1@ then we negate the
-continued fraction; otherwise we leave the continued fraction
-alone.
-
-> realRealAbs :: RealReal -> RealReal
-> realRealAbs x = if realRealSignum x < 0 then realRealNegate x else x
-
-\subsection{Signum}
-
-This isn't computable either, so again {\it caveat emptor}.
-
-> realRealSignum :: RealReal -> RealReal
-> realRealSignum x = wrapRat (signum (makeRational x))
-
-\section{@RealReal@ is an instance of @Real@}
-
-This is yet another non-computable operation that is fudged.
-
-> instance Real RealReal where
->   toRational rx = qNumerator x % qDenominator x
->                   where x = makeRational rx
-
-\section{@RealReal@ is an instance of @Fractional@}
-
-Fortunately these operations are computable.
-
-> instance Fractional RealReal where
->   x / y          = realRealDivide x y
->   fromRational x = wrapRat (numerator x %% denominator x)
-
-> realRealDivide :: RealReal -> RealReal -> RealReal
-> realRealDivide (RealInt x) (RealInt y) = wrapRat (x%%y)
-> realRealDivide (RealInt x) (RealRat y) = wrapRat ((x%%1)/y)
-> realRealDivide (RealInt x) (RealCF  y) = RealCF (algebraicAlgorithm h y)
->                                          where h  = ([0,x],[1,0])
-> realRealDivide (RealRat x) (RealInt y) = wrapRat (x/(y%%1))
-> realRealDivide (RealRat x) (RealRat y) = wrapRat (x/y)
-> realRealDivide (RealRat x) (RealCF  y) = RealCF (algebraicAlgorithm h y)
->                                          where nx = qNumerator   x
->                                                dx = qDenominator x
->                                                h  = ([0,nx],[dx,0])
-> realRealDivide (RealCF  x) (RealInt y) = RealCF (algebraicAlgorithm h x)
->                                          where h  = ([1,0],[0,y])
-> realRealDivide (RealCF  x) (RealRat y) = RealCF (algebraicAlgorithm h x)
->                                          where ny = qNumerator   y
->                                                dy = qDenominator y
->                                                h  = ([dy,0],[0,ny])
-> realRealDivide (RealCF  x) (RealCF  y) = RealCF (quadraticAlgorithm h x y)
->                                          where h  = ([0,0,1,0],[0,1,0,0])
-
-\section{@RealReal@ is an instance of @Floating@}
-
-> instance Floating RealReal where
->   pi    = realRealPi
->   exp   = realRealExp
->   log   = realRealLog
->   sqrt  = realRealSqrt
->   sin   = realRealSin
->   cos   = realRealCos
->   tan   = realRealTan
->   asin  = realRealAsin
->   acos  = realRealAcos
->   atan  = realRealAtan
->   sinh  = realRealSinh
->   cosh  = realRealCosh
->   tanh  = realRealTanh
->   asinh = realRealAsinh
->   acosh = realRealAcosh
->   atanh = realRealAtanh
-
-Here we notice that $\pi$ is a computable real number.
-
-> realRealPi :: RealReal
-> realRealPi = 4 * (RealCF atan1)
-
-\subsection{logarithms and exponentials}
-
-> realRealExp :: RealReal -> RealReal
-> realRealExp (RealInt x) = RealCF exp1 ^^ x
-> realRealExp (RealRat x) = RealCF exp1 ^^ i * RealCF (expQ (x-(i%%1)))
->                           where i = qRound x
-> realRealExp (RealCF  x) = RealCF exp1 ^^ i * RealCF (expR x')
->                           where (i,x') = integerFraction x
-
-> realRealLog :: RealReal -> RealReal
-> realRealLog (RealInt x)
->  = if x==1 then RealInt 0 else RealCF log10 * RealInt i + RealCF (logQ x')
->    where i  = genericLength (show x):: Integer
->          x' = x %% (10^i)
-> realRealLog (RealRat x)
->  = RealCF log10 * RealInt i + RealCF (logQ x')
->    where i   = if x > 1 then f x else negate (f (1/x))
->          f  :: QRational -> Integer
->          f x = genericLength (show (qRound x + 1))
->          x'  = x / ((10%%1)^^i)
-> realRealLog (RealCF  x)
->  = RealCF log10 * RealInt i + RealCF (logR x')
->    where i   = if head x /= 0 then f x else negate (f (tail x))
->          f  :: ContinuedFraction -> Integer
->          f x = genericLength (show (i + 1)) where (i,_) = integerFraction x
->          m   = (1%%10)^^i
->          x'  = algebraicAlgorithm ([qNumerator m,0],[0,qDenominator m]) x
-
-> realRealSqrt :: RealReal -> RealReal
-> realRealSqrt    (RealInt x) = RealCF (cfRat2CFSqrt (x%%1))
-> realRealSqrt    (RealRat x) = RealCF (cfRat2CFSqrt x)
-> realRealSqrt rx@(RealCF  _) = exp (log rx / 2)
-
-\subsection{Trigonometry}
-
-The basic idea for @sin@ and @cos@ is to use the following formul\ae .
-\[\sin(x) = \frac{2 \tan(x/2)}{1+\tan^2(x/2)} \mbox{ and }
-  \cos(x) = \frac{1-\tan^2(x/2)}{1+\tan^2(x/2)}\]
-
-Unfortunately, we observe that there is a problem at the points
-$(2k+1)\pi / 2$ since the @tan@ function switches from $+\infty$ to
-$-\infty$. To cure this we observe the following identities that are
-based on the value of the cotangent function:
-\[\sin(x) = \frac{2 \cot(x/2)}{1+\cot^2(x/2)} \mbox{ and }
-  \cos(x) = \frac{\cot^2(x/2)-1}{\cot^2(x/2)+1}\]
-
-> realRealSin :: RealReal -> RealReal
-> realRealSin rx@(RealCF _)
->  = 2*c/(1+c*c)
->    where (RealCF x2) = rx/2
->          cotx2       = cotR x2
->          tanx2       = tanR x2
->          cf = if head cotx2 == 0 then cotx2 else tanx2
->          c  = RealCF cf
-> realRealSin x = 2*t / (1+t*t)
->                 where t = tan (x/2)
-
-> realRealCos :: RealReal -> RealReal
-> realRealCos rx@(RealCF _)
->  = if head cotx2 == 0 then (c2-1)/(c2+1) else (1-c2)/(1+c2)
->    where (RealCF x2) = rx/2
->          cotx2       = cotR x2
->          tanx2       = tanR x2
->          cf = if head cotx2 == 0 then cotx2 else tanx2
->          c  = RealCF cf
->          c2 = c*c
-> realRealCos x = (1-t2) / (1+t2)
->                 where t  = tan (x/2)
->                       t2 = t*t
-
-> realRealTan :: RealReal -> RealReal
-> realRealTan (RealInt x) = if x==0 then RealInt 0 else RealCF (tanQ (x%%1))
-> realRealTan (RealRat x) = RealCF (tanQ x)
-> realRealTan rx@(RealCF _)
->  = RealCF (tanR x)
->    where (RealCF x') = rx / pi
->          (i,_)       = integerFraction x'
->          (RealCF x)  = rx - (RealInt i * pi)
-
-\subsection{Inverse Trig}
-
-> realRealAsin :: RealReal -> RealReal
-> realRealAsin x = atan (x / sqrt (1-x*x))
-
-> realRealAcos :: RealReal -> RealReal
-> realRealAcos x = atan (sqrt (1-x*x) / x)
-
-> realRealAtan :: RealReal -> RealReal
-> realRealAtan (RealInt x) = if x==0 then RealInt 0 else RealCF (atanQ (x%%1))
-> realRealAtan (RealRat x) = RealCF (atanQ x)
-> realRealAtan (RealCF  x) = RealCF (atanR x)
-
-\subsection{Hyperbolic Trig Functions}
-
-> realRealSinh :: RealReal -> RealReal
-> realRealSinh x = (ex - (1 / ex)) / 2
->                  where ex = exp x
-
-> realRealCosh :: RealReal -> RealReal
-> realRealCosh x = (ex + (1 / ex)) / 2
->                  where ex = exp x
-
-> realRealTanh :: RealReal -> RealReal
-> realRealTanh x = (ex - ex') / (ex + ex')
->                  where ex  = exp x
->                        ex' = 1/ex
-
-\subsection{Inverse Hyperbolic Trig Functions}
-
-> realRealAsinh :: RealReal -> RealReal
-> realRealAsinh x = log (x + sqrt (x*x+1))
-
-> realRealAcosh :: RealReal -> RealReal
-> realRealAcosh x = log (x + sqrt (x*x-1))
-
-This is not quite true as the root needs to be real.
-
-> realRealAtanh :: RealReal -> RealReal
-> realRealAtanh x = log ((1+x)/(1-x)) / 2
-
-\section{@RealReal@ is an instance of @Text@}
-
-> instance Read RealReal where
->   readsPrec p    = readSigned readRealReal
-> instance Show RealReal where
->   showsPrec p rx = showString (showRat p (makeRational rx))
-
-> makeRational :: RealReal -> QRational
-> makeRational (RealInt x) = (x%%1)
-> makeRational (RealRat x) = x
-> makeRational (RealCF  x) = cf2Rat x
-
-\subsection{reading @RealReal@'s}
-
-We can't guarantee that the invariant @read(show x) = x@ once we
-consider computable real numbers -- so we don't!
-
-> readRealReal :: ReadS RealReal
-> readRealReal r = [(RealInt x, s) | (x,s) <- reads r]
-
-\subsection{printing @RealReal@'s}
-
-There's an irritating feature in the Haskell standard prelude that
-requires us to deal with the minus sign here rather than leaving it to
-the standard mechanism.
-
-> showRat :: Int -> QRational -> String
-> showRat p x = if s == "-" && p > 6 then "(" ++ res ++ ")" else res
->               where z     = qRound (x/accuracy)
->                     s     = if z < 0 then "-" else ""
->                     zs    = show (abs z)
->                     l     = length zs
->                     zs'   = pad ++ zs
->                     pad   = if (l <= decimals)
->                             then take (decimals-l+1) zeros else ""
->                     zeros = repeat '0'
->                     (i,f) = splitAt (length zs' - decimals) zs'
->                     res   = s ++ i ++ "." ++ f
-
-
-
diff --git a/spectral/compreals/Transcendentals.lhs b/spectral/compreals/Transcendentals.lhs
deleted file mode 100644 (file)
index 7f6e860..0000000
+++ /dev/null
@@ -1,448 +0,0 @@
-\chapter{Transcendental Functions}
-
-\begin{verbatim}
-$Log: Transcendentals.lhs,v $
-Revision 1.2  1999/11/02 16:10:42  simonpj
-Haskell 98 changes
-
-Revision 1.1  1996/01/08 20:05:19  partain
-Initial revision
-
-\end{verbatim}
-
-> module Transcendentals
->  (module ContinuedFractions, 
->   exp1, expQ, expR, log10, logQ, logR, tanQ, tanR, cotR, atan1, atanQ, atanR)
-> where
-> import ContinuedFractions
-
-\section{Basic Transcendental Functions}
-
-In this section we show how to use Gauss' continued fraction
-representation of $F(\alpha , 1, \gamma ; x)$ to implement the basic
-transcendental functions: $\exp$, $\log$, $\tan$, and $\arctan$.
-
-As a reminder, $F(\alpha , \beta , \gamma ; x)$ is the following series:
-\[1 + \frac{\alpha \beta}{1! \gamma}x +
-\frac{\alpha(\alpha +1) \beta (\beta +1)}{2! \gamma(\gamma +1)}x^2 +
-\frac{\alpha(\alpha +1)(\alpha +2) \beta (\beta +1)(\beta +2)}%
-{3! \gamma(\gamma +1)(\gamma +2)}x^3 +
-\ldots \]
-
-Modifying Gauss' formulation of the continued fraction we get that
-$F(\alpha , 1, \gamma ; x)$ is:
-
-\[ [0,~1,~\frac{1}{\alpha}(\frac{\gamma}{-x}),
-~\frac{\alpha}{\gamma-\alpha}(\gamma+1),~
-p_0(\frac{\gamma +2}{-x}),~p'_0(\gamma +3)~\ldots~
-p_{n}(\frac{\gamma +2n+2}{-x}),~ p'_{n}(\gamma +2n+3),~ \ldots ] \]
-where
-\[p_n = \frac{1}{\alpha}\prod_{i=0}^n \frac{(\gamma - \alpha +i)(i+1)}%
-{(\alpha +i+1)(\gamma +i)}\]
-and
-\[p'_n = \frac{\alpha}{\gamma-\alpha}
-\prod_{i=0}^n \frac{(\alpha +i+1)(\gamma +i)}%
-{(\gamma - \alpha +i+1)(i+2)}\]
-
-This often simplifies.
-
-\subsection{Exponential}
-
-Gauss gives the following representation of $\exp(x)$:
-$\lim_{k\to\infty}F(k,1,1;\frac{x}{k})$. If we expand this we have
-the following continued fraction:
-
-\[ [0,~1,~1.(\frac{-1}{x}),-2,~3.(\frac{1}{x}),~\frac{4}{2},~
-5.(\frac{-1}{x}),~\frac{-6}{3},~\ldots~]\]
-We then observe that $\exp(x) = \frac{1}{\exp(-x)}$ getting:
-\[ [1,~\frac{1}{x},~-2,~\frac{-3}{x},~2,~\frac{5}{x},~-2,~\ldots ] \]
-
-Setting $x$ to $1$ we have the following continued fraction for $e$:
-
-> exp1CF :: [QRational]
-> exp1CF = 1: f 0 where f n = (4*n+1): (-2): (-(4*n+3)): 2: f (n+1)
-
-To use this continued fraction to generate the continued fraction for
-$\exp(x)$, where $x\in {\cal Q}$, we simply multiply the even terms of
-the continued fraction for $e$ by $\frac{1}{x}$ (treating the first
-term as an even one). The starting homography for the modified
-Algebraic Algorithm is generated by absorbing the first six terms of
-the continued fraction. We then multiply the resultant continued
-fraction by $x$. This is marginally more stable than Vuillemin's
-formulation.
-
-> expQ  :: QRational -> ContinuedFraction
-> expQ  x = aaQ hom ts2
->           where (ts1,ts2)   = splitAt 2 cf
->                 hom         = absorbQ startH ts1
->                 (startH,cf) = if abs x < 1 then
->                               (([1,0],[0,1]), oddQ x exp1CF) else
->                               (([qDenominator x,0],[0,qNumerator x]),
->                                 evenQ (1/x) exp1CF) 
-
-> exp1 :: ContinuedFraction
-> exp1 = expQ (1%%1)
-
-We use an exactly similar process for the case of continued fraction
-arguments.
-
-> expR  :: ContinuedFraction -> ContinuedFraction
-> expR  x = if mightBeZero x
->           then processQS x (splitAt 6 (oddR exp1CF))
->           else quadraticAlgorithm ([0,1,1,0],[0,1,0,0]) x cf2
->           where cf2 = processQS x (splitAt 6 (oddR (tail exp1CF)))
-
-\subsection{Logarithm}
-
-Vuillemin's continued fraction is too unstable for values of $x-1$
-which aren't near $0$. So instead we use Gauss' formula:
-\[\log(\frac{1+x}{1-x}) = 2xF(\frac{1}{2},1,\frac{3}{2}; x^2)\]
-It works because the transformation eliminates those values of $x$
-in the above expansion that would cause problems, {\it e.g.} $|x|
-\geq 0$.
-
-The continued fraction for $F(\frac{1}{2},1,\frac{3}{2}; x^2)$ is
-closely related to that for $\arctan(x)$, because:
-\[\tanh^{-1}(x) = \frac{1}{2} \log(\frac{1+x}{1-x})\]
-and
-\[\imath \tanh(x) = \tan(\imath x).\]
-
-> logQ  :: QRational -> ContinuedFraction
-> logQ  x = aaQ startH ts2
->           where (ts1,ts2) = splitAt 2 (oddQ (-t*t) atan1CF)
->                 t         = (x-1)/(x+1)
->                 tn        = qNumerator t
->                 td        = qDenominator t
->                 startH    = absorbQ ([-2*td,0],[0,tn]) ts1
-
-> log10 :: ContinuedFraction
-> log10 = logQ (10%%1)
-
-The nice part is that by requiring the multiplier @-t*t@ to be
-less than @1@, we place no constraints on @x@ provided that it is
-already in the range $(0,\infty )$.
-
-> logR  :: ContinuedFraction -> ContinuedFraction
-> logR  x = quadraticAlgorithm ([2,0,0,0],[0,0,0,1]) t cf1
->           where t   = algebraicAlgorithm ([1,(-1)],[1,1]) x
->                 t2  = quadraticAlgorithm ([-1,0,0,0],[0,0,0,1]) t t
->                 cf  = qsSetupAlt (drop 1 (oddR atan1CF)) t2
->                 cf1 = quadraticAlgorithm ([0,1,0,0],[0,1,1,0]) t2 cf
-
-\subsection{Tangent}
-
-For this function we use Lambert's Continued Fraction
-\[\tan(x) = [0,~\frac{1}{x},~\frac{-3}{x},~\frac{5}{x},~\frac{-7}{x},~\ldots]\]
-
-It is much more convenient if we can obtain a fraction with some
-constant terms in it so we observe the following transformation is
-possible:
-\[\tan(x) = \frac{1}{x}[0,~\frac{1}{x^2},~-3,~\frac{5}{x^2},~-7,~\ldots]\]
-provided that $x$ is non-zero.
-
-> tan1CF :: [QRational]
-> tan1CF = 0: f 0 where f n = (4*n+1) : (-(4*n+3)) : f (n+1)
-
-> tanQ  :: QRational -> ContinuedFraction
-> tanQ  x = if x /= 0 then aaQ startH ts2 else error "tanQ: 0 argument"
->           where n         = qRound (abs (x/2)) + 1
->                 (ts1,ts2) = splitAt (fromIntegral n) (bothQ x tan1CF)
->                 startH    = absorbQ ([1,0],[0,1]) ts1
-
-> tanR  :: ContinuedFraction -> ContinuedFraction
-> tanR  x = quadraticAlgorithm ([0,1,0,0],[0,0,1,0]) cfd cfn
->           where cf  = qsSetupAlt (drop 1 (oddR tan1CF)) x2
->                 x2  = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  x
->                 cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  cf
->                 cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf
-
-> cotR  :: ContinuedFraction -> ContinuedFraction
-> cotR  x = quadraticAlgorithm ([0,1,0,0],[0,0,1,0]) cfn cfd
->           where cf  = qsSetupAlt (drop 1 (oddR tan1CF)) x2
->                 x2  = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  x
->                 cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  cf
->                 cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf
-
-\subsection{Arctangent}
-
-Recall, from the discussion of the continued fraction for $\log$, that
-the $\arctan(x)$ is $xF(\frac{1}{2},1,\frac{3}{2};-x^2)$ 
-The terms in this continued fraction are tending towards $\pi$ and $1$.
-
-> atan1CF :: [QRational]
-> atan1CF = 0: 1: 3: 5/4: f (8/9) 0
->           where f :: QRational -> QRational -> [QRational]
->                 f p n = (7/2 + 2*n)*p: (9/2 + 2*n)*p': f p'' (n+1)
->                         where p'  = 1/(p*(n+2)*(n+2))
->                               p'' = 1/(p'*(n+5/2)*(n+5/2))
-
-The @atanQ@ function is valid for $x \in (-\infty,\,\infty)$.
-
-We note the following identities:
-\[\mbox{for $x>0$ } \arctan(x) = \frac{\pi}{2} -\arctan(\frac{1}{x})\]
-and
-\[\mbox{for $x<0$ } \arctan(x) = -\frac{\pi}{2} -\arctan(\frac{1}{x})\]
-
-
-> atanQ :: QRational -> ContinuedFraction
-> atanQ x
->  = if x > 1  then
->       quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1]) atan1 (cf (1/x)) else
->    if x < -1 then
->       quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1]) atan1 (cf (1/x)) else
->       cf x
->    where
->    cf x = aaQ startH ts2
->           where
->           (ts1,ts2) = splitAt 4 (oddQ (x*x) atan1CF)
->           startH    = absorbQ ([qDenominator x,0],[0,qNumerator x]) ts1
-
-Because we will be using it to help improve the stability of the
-$arctan$ function, we pre-define the constant continued fraction
-$\pi/4$.
-
-> atan1 :: ContinuedFraction
-> atan1 = atanQ (1%%1)
-
-> atanR :: ContinuedFraction -> ContinuedFraction
-> atanR []         = algebraicAlgorithm ([2,0],[0,1]) atan1
-> atanR xs@(x:xs')
->  = if x > 1 then
->       quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1]) atan1 ys else
->    if x < -1 then
->       quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1]) atan1 ys else
->       atanRdirect xs
->    where ys = atanRdirect (algebraicAlgorithm ([0,1],[1,0]) xs)
-
-> atanRdirect :: ContinuedFraction -> ContinuedFraction
-> atanRdirect x
->  = quadraticAlgorithm ([0,0,1,0],[0,1,0,0]) cfn cfd
->    where x2  = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x x
->          cf  = qsSetupAlt (drop 1 (oddR atan1CF)) x2
->          cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x  cf
->          cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf
-
-\subsection{Tests on CFs}
-
-We use @mightBeZero@ to test whether a continued fraction can be zero.
-
-> mightBeZero :: ContinuedFraction -> Bool
-> mightBeZero []     = False
-> mightBeZero (x:xs) = x == 0
-
-\section{Processing Rational Transcendentals}
-
-In this section we deal with the case of transcendental functions of
-rational arguments. We first define a variant of Gosper's Algenraic
-Algorithm that takes a homography and a continued fraction with
-rational terms, and returns an equivalent continued fraction with
-integer terms.
-
-> aaQ :: Homography -> [QRational] -> ContinuedFraction
-> aaQ homography@([n0,n1],[d0,d1]) xs
->  = if n0*d1 /= n1*d0 then aaQNoCheck homography xs
->    else if d1 /= 0   then rat2cf (n1%%d1)
->    else if d0 /= 0   then rat2cf (n0%%d0) 
->    else                   error "aaQ: undefined homography"
-
-As we did with the original formulation, we check the determinant of
-the homography only once.
-
-> aaQNoCheck homography@([n0,n1],[d0,d1]) []
->  = rat2cf (n0%%d0)
-> aaQNoCheck homography@([n0,n1],[d0,d1]) xs'@(x:xs)
->  = case algebraicOutput homography x of
->      Just o  -> o : aaQNoCheck ([d0,d1],          [n0-d0*o,n1-d1*o]) xs'
->      Nothing ->     aaQNoCheck ([n0*n+n1*d,n0*d], [d0*n+d1*d,d0*d])  xs
->    where (n,d) = (qNumerator x, qDenominator x)
-
-To absorb a finite sequence of terms from a rational continued fraction
-we use @absorbQ@.
-
-> absorbQ :: Homography -> [QRational] -> Homography
-> absorbQ homography        []     = homography
-> absorbQ ([n0,n1],[d0,d1]) (q:qs) = absorbQ ([n0*n+n1*d,n0*d],
->                                             [d0*n+d1*d,d0*d]) qs
->                                    where n = qNumerator   q
->                                          d = qDenominator q
-
-To set up the continued fraction with the argument we use @bothQ@,
-@oddQ@ and @evenQ@.
-
-> bothQ :: QRational -> [QRational] -> [QRational]
-> bothQ x (x1:xs) = x1/x:   bothQ x xs
-
-> evenQ :: QRational -> [QRational] -> [QRational]
-> evenQ x (x1:x2:xs) = x1/x: x2:   evenQ x xs
-
-> oddQ :: QRational  -> [QRational] -> [QRational]
-> oddQ  x (x1:x2:xs) = x1:   x2/x: oddQ  x xs
-
-\section{Processing Irrational Transcendentals}
-
-To produce continued fractions for the transcendental functions with
-continued fraction arguments we use an ``infinite state machine''.
-This infinite state is an infinite list of sub-states, each consisting
-of a rational, and a homography for the Quadratic Algorithm.
-
-> type State = (QRational, Homography)
-
-To set up the initial continued fraction for the ``infinite state
-machine'' from the initial continued fraction we use @oddR@.
-If the initial continued fraction is
-\[[\frac{a_0}{b_0},~\frac{a_1}{b_1},~\frac{a_2}{b_2},~\frac{a_3}{b_3},~\ldots]\]
-then the ``infinite state'' will be:
-\[[\left(\frac{a_0}{b_0},
-   \left(\begin{array}[c]{cccc}0&a_1&b_1&0\\b_1&0&0&0\end{array}\right),~
-\right)
-\left(\frac{a_2}{b_2},
-   \left(\begin{array}[c]{cccc}0&a_3&b_3&0\\b_3&0&0&0\end{array}\right)
-\right),~
-\ldots]\]
-
-> oddR :: [QRational] -> [State]
-> oddR (x1:x2:xs) = (x1,([0,xn,xd,0],[xd,0,0,0])): oddR xs
->                   where (xn,xd) = (qNumerator x2, qDenominator x2)
-
-The processing of the infinite state and the continued fraction is
-initiated by @processQS@.
-
-> processQS :: ContinuedFraction -> ([State],[State]) -> ContinuedFraction
-> processQS x (hs1,hs2) = preProcess x hs1 (qsSetupAlt hs2 x)
-
-This in turn makes use of @preProcess@ which is used to absorb some of
-the leading terms (@hs1@) and @qsSetupAlt@ which is used to process
-the remainder (@hs2@).
-
-> preProcess :: ContinuedFraction -> -- the argument CF
->               [State]           -> -- leading parts of infinite state
->               ContinuedFraction -> -- accumulating CF
->               ContinuedFraction    -- the resulting CF
-> preProcess _ []         c = c
-> preProcess x ((i,h):hs) c = algebraicAlgorithm ([ni,nd],[nd,0]) 
->                             (quadraticAlgorithm h x (preProcess x hs c))
->                             where ni = qNumerator i
->                                   nd = qDenominator i
-
-To process the remainder of the infinite state, we use @qsSetupAlt@.
-
-> qsSetupAlt :: [State] -> ContinuedFraction -> ContinuedFraction
-> qsSetupAlt hs x = aaQ ([1,0],[0,1]) (qsProcess (hs, x))
-
-\subsection{Quadratic Process}
-
-In essence the state for the algorithm holds:
-\begin{enumerate}
-\item the infinite list of @Rational@/@Homography@ pairs; and
-\item the infinite continued fraction of the argument.
-\end{enumerate}
-
-We advance by absorbing one term from the continued fraction or by
-considering one more term of the infinite state.
-
-> qsProcess :: ([State], ContinuedFraction) -> [QRational]
-> qsProcess    (ihs, [])    = asProcess ihs
-> qsProcess st@(ihs, (_:_)) = qsNext 1 st
-
-The decision about whether sufficient accuracy has been achieved is
-decided by @qsNext@. To do this a call to @qsTerm@ is made. This
-returns a list of output values and a modified state (@ihs'@).  If the
-output list is sufficiently long we can output them as part of the
-resulting continued fraction. Alternatively we can absorb another term
-of the argument continued fraction. In either case we consider one
-more term of the infinite state.
-
-> qsNext :: Int -> ([State], ContinuedFraction) -> [QRational]
-> qsNext n (ihs, []) = asProcess (qsEnd ihs)
-> qsNext n st@(((i,h):ihs), xs'@(x:xs))
->  = if length outs > 1
->       then init outs ++ qsNext (n+1) (ihs',xs')
->       else qsNext (n+1) (qsIn x ((i,h): ihs), xs)
->    where (outs, (ihs',_)) = qsTerm n st
-
-The @qsTerm@ function considers the initial @n@ terms of the infinite
-state, and the first term of the continued fraction for the argument.
-From this it calculates a modified infinite state, in which, what is
-effectively carry propogation has occurred.
-
-> qsTerm :: Int -> ([State], ContinuedFraction)
->          -> ([QRational],([State], ContinuedFraction))
-> qsTerm n (ihs,xs)
->  = (outs, (ihs',xs))
->    where (front,back@((i,_):_)) = splitAt n ihs
->          (outs,ihs') = foldr (qsTerm' (head xs)) ([i],back) front
-
-To deal with each individual term of the infinite state we use @qsTerm'@.
-
-> qsTerm' :: Integer -> State -> ([QRational], [State]) -> 
->                                ([QRational], [State])
-> qsTerm' x ih@(i,h) (is,ihs)
->  = (is',((last is',h''):ihs))
->    where h' = qsAbsorb (init is) h
->          i' = last is
->          (is',h'') = qsEmit x i' ([i],h')
-
-This in turn uses @qsAbsorb@ and @qsEmit@ to propogate the carry and
-emit output terms, respectively.
-
-> qsAbsorb :: [QRational] -> Homography -> Homography
-> qsAbsorb [] h = h
-> qsAbsorb (i:is) h@([n0,n1,n2,n3],[d0,d1,d2,d3])
->  = qsAbsorb is ([n0*n+n2*d, n1*n+n3*d, n0*d, n1*d],
->                 [d0*n+d2*d, d1*n+d3*d, d0*d, d1*d])
->    where n = qNumerator i
->          d = qDenominator i
-
-@qsEmit@ is essentially a variant of the Quadratic Algorithm, capable
-of dealing with the output of more than $1$ term at once.
-
-> qsEmit :: Integer -> QRational -> ([QRational], Homography) ->
->                                 ([QRational], Homography)
-> qsEmit x i ish@(is,h@([n0,n1,n2,n3],[d0,d1,d2,d3]))
->  = case quadraticOutput h [x%%1,i] of
->      Just o  -> qsEmit x i (is++[o%%1],
->                             ([d0,d1,d2,d3],
->                              [n0-d0*o, n1-d1*o, n2-d2*o, n3-d3*o]))
->      Nothing -> ish
-
-To absorb the leading term of the argument continued fraction into
-each Quadratic Homography we use @qsIn@.
-
-> qsIn :: Integer -> [State] -> [State]
-> qsIn = map . input
->        where input x (i,([n0,n1,n2,n3],[d0,d1,d2,d3]))
->                = (i, ([n0*x+n1, n0, n2*x+n3, n2],
->                       [d0*x+d1, d0, d2*x+d3, d2]))
-
-\subsection{Finite Continued Fractions}
-
-It is possible to generate finite continued fractions, so we need a
-way to deal with @[]@ in the @qsProcess@. This results in a
-degeneration to the algebraic process.
-
-@qsEnd@ turns each Quadratic Homography into the equivalent Algebraic
-Homography.
-
-> qsEnd :: [State] -> [State]
-> qsEnd = map input
->         where input (i,([n0,_,n2,_],[d0,_,d2,_]))
->                 = (i,([n0,n2],[d0,d2]))
-
-The function @asProcess@ munges these homographies to produce an answer.
-
-> asProcess :: [State] -> [QRational]
-> asProcess ihs = fst (head ihs) : asProcess (asNext ihs)
-
-> asNext :: [State] -> [State]
-> asNext ((i,h@([n0,n1],[d0,d1])):ihs'@((i',_):ihs))
->  = case algebraicOutput h i' of
->      Just o  -> (o%%1, ([d0,d1], [n0-d0*o,n1-d1*o])): ihs'
->      Nothing -> asNext ((i,hIn i' h): asNext ihs')
->    where hIn x ([n0,n1],[d0,d1]) = ([n0*xn+n1*xd,n0*xd],[d0*xn+d1*xd,d0*xd])
->                                    where xn = qNumerator x
->                                          xd = qDenominator x
-
-And that completes the implementation of the specified subset of
-transcendental functions.
-
-
diff --git a/spectral/compreals/paper.tex b/spectral/compreals/paper.tex
deleted file mode 100644 (file)
index 52f3106..0000000
+++ /dev/null
@@ -1,1296 +0,0 @@
-\documentstyle{jfp}
-%\documentstyle[11pt,fpart,fpbib]{article}
-% maths stuff
-\newtheorem{theorem}{Theorem}[section]
-\newtheorem{definition}[theorem]{Definition}
-\newtheorem{example}[theorem]{Example}
-\newtheorem{counterexample}[theorem]{Counter Example}
-\newtheorem{lemma}[theorem]{Lemma}
-\newtheorem{hypothesis}[theorem]{Hypothesis}
-\newtheorem{proposition}[theorem]{Proposition}
-\newtheorem{corollary}[theorem]{Corollary}
-%\def\proof[#1]{\begin{flushleft} \mbox{\bf Proof of #1} \end{flushleft}
-%               \begin{list}{}{} \item[]}
-%\def\endproof{\end{list}
-%              \begin{flushright} \mbox{$\Box$} \end{flushright}}
-
-%\def\author{David Lester}
-%\def\fromemail{dlester@cs.man.ac.uk}
-%\def\doctype{Report-0004}
-%\def\title{Computable Real Numbers for {\sc Haskell}}
-
-% Filled in on creation of document: \author, \fromemail, \doctype and \title 
-
-% RCS amended stuff
-%$Log: paper.tex,v $
-%Revision 1.1  1996/01/08 20:05:20  partain
-%Initial revision
-%
-%Revision 1.1  94/06/15  11:50:41  dlester
-%Initial
-%
-%$Author: partain $
-\def\fpstate{Draft}
-\def\fprevision{1.1}
-\def\fpdate{94/06/15} \def\fptime{11:50:41}
-%\def\acknowledge{}      % delete this if you have no affiliation acks
-%\def\fromaddress{Functional Programming Group,\\Department of Computer Science,
-%Manchester University,\\Oxford Road,
-%Manchester M13 9PL, UK.}
-
-%\rcsheadings{}
-\def\indexTT#1{}
-\def\indexDTT#1{}
-
-\begin{document}
-\title[Computable Reals]{Computable Real Numbers for {\sc Haskell}}
-\author[David Lester]{David Lester\\Functional Programming Group,
-Department of Computer Science,\\ Manchester University, Oxford Road,
-Manchester M13 9PL, UK}
-
-\bibliographystyle{fp}
-\maketitle
-%\fptitle{}
-\begin{abstract}
-
-This work arose from a desire to provide an implementation of
-computable real numbers for {\sc Haskell}. The starting point was
-Vuillemin's report on continued fraction arithmetic
-\shortcite{VUILLEMIN87a}. One advantage that the present
-implementation has over Vuillemin's original work is that -- being
-implemented in a lazy functional language -- we are able to take
-advantage of laziness to ensure that the stream processing takes place
-in a computable way.
-
-This paper goes into some of the details that were absent from
-Vuillemin's original work. Further, we show how to take advantage of
-Gauss' work on continued fractions, to generate the transcendental
-functions: \mbox{\tt exp}\indexTT{exp}, \mbox{\tt log}\indexTT{log}, \mbox{\tt 
-tan}\indexTT{tan}, and \mbox{\tt atan}\indexTT{atan}. This requires us to
-process continued fractions whose terms are themselves continued
-fractions; we show how to achieve this.
-
-\end{abstract}
-%\keywords{Computable Real Numbers, Computable Transcendental
-%          Functions, Continued Fractions, {\sc Haskell}}
-
-\section{Introduction}
-
-If we are content with using only rational numbers, accurate computer
-arithmetic is easy to implement. Unfortunately, as was discovered
-early in the history of mathematics, there are numbers, such as
-$\sqrt{2}$, which are {\em irrational}. The discovery of computable
-functions by Turing \shortcite{TURING36a} and Church
-\shortcite{CHURCH41a} led to work on computable numbers
-\cite{RICE54a,BISHOP85a}. This work used a radix representation along
-with intervals. This is fine from the point of correctness, but as the
-implementation is {\em memoryless}, it could be quite inefficient.  In
-their paper \shortcite{Boehm90a}, Boehm and Cartwright make use of a
-higher order implementation to ensure that a rational implementation
-carries out calculations to sufficient accuracy for the printing of
-the answer.
-
-This paper was written in an attempt to understand and interpret the
-work of Jean Vuillemin \shortcite{VUILLEMIN87a}, in which he describes a
-representation of computable real numbers using continued fractions.
-In this paper we will hopefully see the benefits of using a lazy
-functional language to implement algorithms that are intrinsically
-lazy.  We also explore some of the algorithms in more detail than
-\cite{VUILLEMIN87a}, paying particular care to ensure that laziness is
-preserved.
-
-The advantage with using continued fractions as a representation of
-numbers is that if, for example, we calculate $\pi$ to a high degree
-of accuracy at some stage, this accuracy is then available everywhere
-else that $\pi$ is used without further calculation being necessary.
-
-Because we will be using the terms frequently, and because the
-concepts may not be familiar, the introduction continues by informally
-introducing three key concepts used in this paper. They are {\em Real
-Numbers}, {\em Computable Real Numbers}, and {\em Continued
-Fractions}. Readers already familiar with this material should feel
-free to skip the remainder of the introduction.
-
-In Section~\ref{sc:reals} we discuss computable real numbers, and in
-Section~\ref{sc:computability} we investigate which
-computable/non-computable functions we are required to implement in
-{\sc Haskell}. In Section~\ref{sc:comp-cf} computable continued
-fractions are discussed.  A gentle introduction to Gosper's algorithm
-is provided in Section~\ref{sc:aa}; the main algorithm is given in
-Section~\ref{sc:qa}. The continued fraction representation of the
-exponential, logarithm, tangent and arctangent functions is given in
-Section~\ref{sc:tf}.
-
-\section{What are Real Numbers?}
-\label{sc:reals}
-
-The best way to think of the various stages in the development of
-number theory is that of a closure operation. If you start with the
-integer one, and the operation of $+$ you can generate all of the
-positive natural numbers ${\cal N}^{+}$. Augmenting the set of
-operations with $-$ it is possible to generate all of the integers,
-${\cal Z}$. We may therefore make the statement: the closure of
-operations $+$ and $-$ with the number one is ${\cal Z}$, the
-integers. Multiplication adds no new numbers, but when we consider the
-division operation we must include fractions. These constitute the
-rational numbers, which we represent by the symbol ${\cal Q}$. We can
-restate this as: the closure of the operations of $+$, $-$, $\times$,
-and $/$ on the number one is ${\cal Q}\,\cup\,\{\bot,\,\infty\}$.
-(We use $\bot$ to represent the number $0/0$ and $\infty$ to
-represent $n/0$, whenever $n\not=0$).
-
-If we think in terms of a pencil and paper implementation for the
-moment, then provided that we have a sufficiently large piece of paper
-we can perform operations on rational numbers to any desired degree of
-accuracy. In the context of the implementation of rational number
-arithmetic, this can be paraphrased as: provided that the memory does
-not overflow we can perform operations on rational numbers to any
-desired degree of accuracy.
-
-But of course -- as the Greek Mathematicians discovered\footnote{The
-theory of irrationals in \cite{EUCLID} is attributed to The\ae{}tus.}
--- this is not the whole story. There exist numbers (such as
-$\sqrt{2}$) that cannot be represented by any finite fraction. Such
-numbers are said to be {\em irrational}. There are a many beautiful
-generalisations that arise from this, but for the moment we simply
-observe that numbers such as square roots arise from finding the roots
-of quadratic polynomials. The rational and irrational numbers that are
-generated as roots to polynomial equations are called {\em algebraic
-numbers}.
-
-At this stage you may wonder where this process stops; {\it i.e.\/}
-how many operations must we use to generate all of the real numbers?
-The answer is that the process doesn't stop and no finite set of
-operations will generate all of the real numbers. The solution is to
-restrict attention to a subset of ${\cal R}$ which I shall call ${\cal
-A}$; which has only a limited set of operations.  The set ${\cal A}$
-will constitute a useful set of computable real numbers for
-engineering applications.
-
-\[\begin{array}{rcl}
-1 & \in & {\cal A}\\
-{\cal A}+{\cal A}& \subset & {\cal A}\\
--{\cal A}& \subset & {\cal A}\\
-{\cal A}/{\cal A}_{\not= 0}& \subset & {\cal A}\\
-e^{{\cal A}}& \subset & {\cal A}\\
-\log {\cal A}_{\, > 0}& \subset & {\cal A}
-\end{array}\]
-
-Notice that by removing the restrictions on the argument for the
-logarithm function we would obtain a subset of the complex numbers,
-${\cal C}$, which would include trigonometric functions. This should
-provide a sufficient range of basic arithmetic operations for most
-engineering applications. The next subsection high-lights the major
-problem involved in producing an arithmetic package for real numbers.
-
-\subsection{What are Computable Real Numbers?}
-
-Informally, a number is said to be computable if we can write a
-program which will print its value to some specified degree of
-accuracy. Following Pour-El and Richards \shortcite{POUR-EL89a}, we may
-formalise this definition.
-
-\begin{definition}\rm
-A real number $x$ is {\em computable} if there exists a computable
-sequence $\{r_k\}$ of rationals which converges effectively to $x$.
-\end{definition}
-
-\begin{definition}\rm
-A sequence $\{r_k\}$ of rational numbers is {\em computable} if there
-exist three recursive functions $a$, $b$, and $s$ from ${\cal N}$ to
-${\cal N}$ such that $b(k) \not= 0$ for all $k$ and
-\[r_k = (-1)^{s(k)}\frac{a(k)}{b(k)} \mbox{ for all } k.\]
-\end{definition}
-
-In other words: an infinite {\sc Haskell} sequence, not all of whose
-elements are $\bot$.
-
-\begin{definition}\rm
-A sequence $\{r_k\}$ of rational numbers {\em converges effectively}
-to a real number $x$ if there exists a recursive function $e:{\cal
-N}\to{\cal N}$ such that for all $N$:
-\[k \geq e(N) \mbox{ implies } |r_k - x| \leq 2^{-N}.\]
-\end{definition}
-
-This definition of computable real number is (inefficiently)
-implementable directly in {\sc Haskell}. What we require is an
-infinite list of \mbox{\tt Rational}\indexTT{Rational} numbers, not all of whic
-h are undefined;
-and a modulus of convergence function: \mbox{\tt e} from \mbox{\tt Integer}\indexTT{Integer} to
-\mbox{\tt Integer}\indexTT{Integer}.
-
-A cardinality argument suffices to show that non-computable numbers
-exist: the set of programs is countable, the set of real numbers is
-not. It seems reasonable that in dealing with numbers for computers,
-we should restrict our attention to just the subset of real numbers
-that is computable. In fact for the engineering applications that we
-have in mind this set may be further restricted to the closure of the
-arithmetic operations on the singleton set: $\{1\}$ -- as outlined
-above.
-
-Life is seldom simple. It is an unfortunate fact of life that it is
-not even possible to tell whether a computable real number is $0$.
-Let us suppose that we are given a monotonically decreasing sequence
-of rationals that converges to $0$. One such sequence is $1$, $0.1$,
-$0.01$, $0.001$, \ldots; proving the convergence establishes that this
-sequence represents the computable real number $0$.  Similarly we
-might wish to consider the monotonically increasing sequence: $-1$,
-$-0.1$, $-0.01$, $-0.001$, \ldots; again this is another
-representation of the real number $0$. But notice that the application
-of the \mbox{\tt signum}\indexTT{signum} function to the sequence gives.
-\[\mbox{\tt map}\indexTT{map}\,\mbox{\tt signum}\indexTT{signum}\,[-0.1,\, -0.0
-1,\, -0.001,\,\ldots] = 
-[-1,\, -1,\, -1,\,\ldots]\]
-But
-\[\mbox{\tt signum}\indexTT{signum}\, (\lim\{-0.1,\, -0.01,\, -0.001,\,\ldots\}
-) = 0\]
-
-There are more operations that are not computable; some of these are
-listed in Theorem~\ref{th:rice}, which is a combination of results
-from Turing \shortcite{TURING36a} and Rice \shortcite{RICE54a}. Proofs
-of these results may be found in \cite{VUILLEMIN87a}.
-
-\begin{theorem}[{Turing, Rice}]\label{th:rice}
-
-There is no algorithm, which for an arbitrary computable real number
-$r\in{\cal R}$, can compute:
-\begin{enumerate}
-\item whether the number $r$ is zero;
-\item whether $r > r'$ for any fixed $r'\in {\cal R}$;
-\item whether the number is rational;
-\item the first digit of the decimal expansion of $r$;
-\item the first term of the regular continued fraction expansion of
-$r$ (which is $\lfloor r\rfloor$);
-\item the value $f(r)\in {\cal R}$ of any function $f$ which is not
-continuous at $r$.
-\end{enumerate}
-\end{theorem}
-
-This is all very disturbing. However, facing facts is infinitely
-preferable to facing defeat, so we must cast around for ways to avoid
-the Turing/Rice Theorems.
-
-The first ``fix'' is reasonably obvious: we must give up any idea that
-exact comparisons are available. They are replaced by what I shall
-refer to as $\varepsilon$-comparisons. Two numbers $r$ and $r'$ are
-said to be ${=}_\varepsilon$ if $|r-r'| < \varepsilon$. Similar
-properties hold for the other comparison operations. This is not so
-very different to the properties of comparison operations available in
-floating point packages.
-
-The second ``fix'' is implied by the fourth observation in
-Theorem~\ref{th:rice}, where we notice that we can't even find the
-first digit of a computable number. To see this consider the sequence
-$\{0.9,\,0.99,\,0.999,\,\ldots\}$.  The first digit of each rational
-approximation is $0$, but of course the first digit of the number $1$
-is $1$. The solution: use a redundant representation of the computable
-numbers. This is implicit in the decimal representation of numbers: as
-we have seen $0.\dot{9}$ and $1$ are both representations of the same
-number.
-
-We now investigate precisely which operations we are required to
-implement for {\sc Haskell}.
-
-\section{{\sc Haskell} and Computability}
-\label{sc:computability}
-
-One of the more irritating features of {\sc Haskell}'s numeric type classes
-is that it makes no distinction between those methods that are
-computable and those which are not. It would be useful if we could at
-least ensure that each class includes methods that either all
-computable or all non-computable. As Table~\ref{tb:comp} shows, this
-is not likely to meet with general approval, because some inverses are
-non-computable {\it e.g.} \mbox{\tt read}\indexTT{read}/\mbox{\tt show}\indexTT
-{show}.
-
-\begin{table}[htbp]
-\begin{tabular}{l|p{2in}|p{1.5in}}
-\hline
-{\tt Class} & {\tt Computable} & {\tt Non-computable} \\
-\hline
-{\tt Eq} & & \mbox{\tt ==}, \mbox{\tt /=} \\
-{\tt Ord} & \mbox{\tt max}\indexTT{max}, \mbox{\tt min}\indexTT{min} & \mbox{\t
-t <=}, \mbox{\tt <}, \mbox{\tt >}, \mbox{\tt >=}  \\
-{\tt Text} & \mbox{\tt readSigned}\indexTT{readSigned} & \mbox{\tt showSigned}\
-indexTT{showSigned} \\
-{\tt Num} & \mbox{\tt +}, \mbox{\tt -}, \mbox{\tt *}, \mbox{\tt negate}\indexTT
-{negate}, \mbox{\tt abs}\indexTT{abs}, \mbox{\tt fromInteger}\indexTT{fromInteg
-er} & \mbox{\tt signum}\indexTT{signum} \\
-{\tt Real} & & \mbox{\tt toRational}\indexTT{toRational} \\
-{\tt Fractional} & \mbox{\tt /}, \mbox{\tt fromRational}\indexTT{fromRational} 
-& \\
-{\tt Floating} & \mbox{\tt pi}, \mbox{\tt exp}\indexTT{exp}, \mbox{\tt log}\indexTT{log}, \mbox{\tt sqrt}\indexTT{sqrt}, \mbox{\tt **}, \mbox{\tt logBase}\indexTT{logBase}, \mbox{\tt sin}\indexTT{sin}, \mbox{\tt cos}\indexTT{cos}, \mbox{
-\tt tan}\indexTT{tan},
-                 \mbox{\tt asin}\indexTT{asin}, \mbox{\tt acos}\indexTT{acos}, 
-\mbox{\tt atan}\indexTT{atan}, \mbox{\tt sinh}\indexTT{sinh}, \mbox{\tt cosh}\indexTT{cosh}, \mbox{\tt tanh}\indexTT{tanh}
-                 \mbox{\tt asinh}\indexTT{asinh}, \mbox{\tt acosh}\indexTT{acosh}, \mbox{\tt atanh}\indexTT{atanh} & \\
-\hline
-\end{tabular}
-\caption{Computable and non-computable Classes is Standard {\sc Haskell}}
-\label{tb:comp}
-\end{table}
-
-\subsection{Class \mbox{\tt Ord}\indexTT{Ord}}
-
-It is a shame that the standard prelude defines \mbox{\tt max}\indexTT{max} and
- \mbox{\tt min}\indexTT{min} in
-terms of the \mbox{\tt <=} function; from the point of view of computability it
-would be more natural for the default methods to be:
-\begin{verbatim}
-> x <= y  = x == min x y
-> x <  y  = x <= y && x /= y
-> x >  y  = y < x
-> x >= y  = y <= x
-\end{verbatim}
-This properly places the blame for the non-computability onto the test
-for equality and not on the ordering itself.
-
-\subsection{Class \mbox{\tt Text}\indexTT{Text}}
-
-The first, and most obvious, problem is that we must forego the
-condition that \mbox{\tt read(show\ x)\ =\ x}. This is a natural consequence of
-the fact that we have printed out only a finite approximation to \mbox{\tt x}
-using \mbox{\tt show}\indexTT{show}.
-
-There is a problem with \mbox{\tt showSigned}\indexTT{showSigned}. It is define
-d as:
-\begin{verbatim}
-> showSigned f p x = showParen (x < 0 && p > 6) (f x)
-\end{verbatim}
-\indexDTT{showSigned}%
-What happens when \mbox{\tt x} is \mbox{\tt 0}? We first observe that the compa
-rison
-\mbox{\tt x\ <\ 0} is performed using a rational approximation to \mbox{\tt x},
- which might
-not be \mbox{\tt 0}, but some small negative number. The \mbox{\tt f} function,
- however
-merely outputs \mbox{\tt "0.000000"}, and we have needlessly added
-parenthesis. The cure is to rewrite the function from scratch.
-
-\subsection{Class \mbox{\tt Num}\indexTT{Num}}
-
-The problem with this class revolve aroung the defintion of \mbox{\tt abs}\indexTT{abs} and
-\mbox{\tt signum}\indexTT{signum}. One choice is to define \mbox{\tt abs}\indexTT{abs} as:
-\begin{verbatim}
-> abs x = max x (-x)
-\end{verbatim}
-\indexDTT{abs}%
-which is then computable. The problem with \mbox{\tt signum}\indexTT{signum} is
- (again) that we
-must test \mbox{\tt x} for equality with \mbox{\tt 0}. A consequence of this de
-finition
-is that it is no longer true that: \mbox{\tt abs\ x\ *\ signum\ x\ =\ x}, which
- is
-unfortunate.
-
-The alternative which I have chosen to implement preserves the
-condition that \mbox{\tt abs\ x\ *\ signum\ x\ =\ x}, by calculating the
-approximation to \mbox{\tt signum\ x} first and then using this to decide
-whether to negate \mbox{\tt x}.
-
-\subsection{Classes \mbox{\tt Fractional}\indexTT{Fractional} and \mbox{\tt Flo
-ating}\indexTT{Floating}}
-
-Provided that the arguments remain within the domain of definition for
-each operator, all the operations in these classes are computable.
-
-There are however a number of extra default methods that could be
-added for the \mbox{\tt Real}\indexTT{Real} case.
-\begin{verbatim}
-> pi = 4 * atan 1
-\end{verbatim}
-The basic idea for \mbox{\tt sin}\indexTT{sin} and \mbox{\tt cos}\indexTT{cos} 
-is to use the following formul\ae .
-\[\sin(x) = \frac{2 \tan(x/2)}{1+\tan^2(x/2)} \mbox{ and }
-  \cos(x) = \frac{1-\tan^2(x/2)}{1+\tan^2(x/2)}\]
-
-Unfortunately, we observe that there is a problem at the points
-$(2k+1)\pi / 2$ since the \mbox{\tt tan}\indexTT{tan} function switches from $+
-\infty$ to
-$-\infty$. To cure this we observe the following identities that are
-based on the value of the cotangent function:
-\[\sin(x) = \frac{2 \cot(x/2)}{1+\cot^2(x/2)} \mbox{ and }
-  \cos(x) = \frac{\cot^2(x/2)-1}{\cot^2(x/2)+1}\]
-
-We use these alternative definitions when we expect that $x$ might be
-close to $(2k+1)\pi /2$.
-\begin{verbatim}
-> asin x = atan (x / sqrt (1-x*x))
-> acos x = atan (sqrt (1-x*x) / x)
-> sinh x = (ex - (1 / ex)) / 2 where ex = exp x
-> cosh x = (ex + (1 / ex)) / 2 where ex = exp x
-\end{verbatim}
-\indexDTT{asin}%
-\indexDTT{acos}%
-\indexDTT{sinh}%
-\indexDTT{cosh}%
-It is also marginally more efficient to calculate \mbox{\tt tanh}\indexTT{tanh}
- using the
-following identity rather than the one in the standard prelude, which
-involves the recomputation of \mbox{\tt exp\ x}.
-\begin{verbatim}
-> tanh x = (ex - ex') / (ex + ex') where ex  = exp x
->                                        ex' = 1/ex
-\end{verbatim}
-\indexDTT{tanh}%
-\begin{verbatim}
-> asinh x = log (x + sqrt (x*x+1))
-> acosh x = log (x + sqrt (x*x-1))
-> atanh x = log ((1+x)/(1-x)) / 2
-\end{verbatim}
-\indexDTT{asinh}%
-\indexDTT{acosh}%
-\indexDTT{atanh}%
->From the foregoing we can see that the transcendental functions
-required of a {\sc Haskell} implementation are reduced to: \mbox{\tt exp}\index
-TT{exp},
-\mbox{\tt log}\indexTT{log}, \mbox{\tt tan}\indexTT{tan} (and \mbox{\tt cot}\indexTT{cot}), and \mbox{\tt atan}\indexTT{atan}.
-
-We now move on to consider continued fractions.
-
-\section{What are Continued Fractions?}
-In this section we will see how to represent numbers by continued fractions.
-Continued fractions are an attractive representation of real numbers
-because: we are able to take advantage of what functional programmers
-would call {\em laziness}, and there exist algorithms to perform the
-arithmetic operations we require.
-
-\begin{definition}\rm A finite continued fraction representation of
-the number $r$ is a sequence of numbers $[x_0,\,x_1,\,\ldots \, x_n]$,
-satisfying:
-\[r=x_0+\frac{1}{x_1 +\frac{1}{\ddots + \frac{1}{x_n}}}\]
-\end{definition}
-
-There is no reason why we should only consider finite continued
-fractions, so infinite fractions are also permitted.
-
-For the time being, in this paper, we will make the numbers in the
-sequence $x_0\ldots x_n$ be integers, where it is to be understood
-that these integers may be arbitrarily large. A finite continued
-fraction with integer terms can be used as a compact method to represent
-a rational number.  It is possible to calculate the values of the
-integers $x_0,\,x_1,\,\ldots \, x_n$ from $r$ in the following way.
-Let $r = r_0$; we compute
-\[r_0 = x_0 + \frac{1}{r_1},\,
-  r_1 = x_1 + \frac{1}{r_2},\,\ldots,\,r_n = x_n \]
-where at each step we choose an integer $x_n$ that is {\em close} to
-$r_n$. Let us give $x_n$ the value $x_n = f(r_n)$. There are a number
-of choices for the function $f$. Obvious candidates are the rounding
-operations.
-\[\begin{array}{llllll}
-\mbox{{\bf floor}} &\lfloor r\rfloor & = & z&\mbox{such that} &
-0 \leq r-z < 1;\\
-\mbox{{\bf round}} &\lfloor r\rceil & = &z&\mbox{such that} &
--\frac{1}{2} < r-z \leq \frac{1}{2};\\
-\mbox{{\bf ceiling}} &\lceil r\rceil & = &z&  \mbox{such that} &
--1 < r-z \leq 0.
-\end{array}\]
-The choice of \mbox{{\bf floor}} leads to Gosper's {\bf N}-fractions.
-As an alternative, we can give the {\bf Z}-fractions for the following
-numbers, which are obtained by using the \mbox{\tt round}\indexTT{round} operat
-ion to generate
-the continued fraction.
-\[\begin{array}{lll}
-\sqrt{2} & = & [1,\,2,\,\ldots,\,2,\,\ldots]\\
-\sqrt{3} & = & [2,\,-4,\,4,\,\ldots,\,-4,\,4,\,\ldots]\\
-\phi     & = & [2,\,-3,\,3,\,\ldots,\,-3,\,3,\,\ldots]\\
-e & = & [3,\,-4,\,2,\,5,\,-2,\,-7,\,\ldots,\,2,\,4n+5,\,-2,\,-4n-7,\,\ldots]
-\end{array}\]
-The interested reader can use the definitions given earlier, in
-conjunction with a series expansion for each irrational number, to
-prove that these continued fractions are correct.
-
-\section{Computable Continued Fractions}
-\label{sc:comp-cf}
-
-Although continued fractions are a good representation because they
-can be used lazily, there remains a problem with selecting either {\bf
-N}-fractions or {\bf Z}-fractions as representations of real numbers:
-neither {\bf floor} nor {\bf round} are computable operations. This
-means that for some numbers it is not possible for a computer to
-construct the continued fraction, even though it exists. What is
-needed is a computable analogue to {\bf floor} or {\bf round}; this is
-implicitly defined by the function \mbox{\tt cfBound}\indexTT{cfBound}, which i
-ndicates the
-range in which a continued fraction must lie given its leading term
-\mbox{\tt x}.
-\begin{verbatim}
-> cfBound :: QRational -> Interval
-> cfBound x
->  = if qInfinite x then (1%%0,   1%%0)    else
->    if abs x < 2   then (-2%%1,  2%%1)    else
->    if x < 0       then (1%%2-x, x+1%%2)  else
->                        (x-1%%2, -x-1%%2)
-\end{verbatim}
-\indexDTT{cfBound}%
-\indexDTT{cfBound}%
-Notice that we have had to extend {\sc Haskell}'s \mbox{\tt Rational}\indexTT{R
-ational} numeric
-type with infinite values, so that we can effectively deal with open
-intervals of the form: $(x,\,\infty )$. It is also worth noting that
--- following Vuillemin's lead -- open intervals of the form $(x,\,y)$
-in which $x>y$ are to be interpreted as the interval: $(-\infty
-,\,y)\,\cup\,(x,\,\infty )$
-
-The important point to note about {\bf E}-fractions is that although
-we have weakened the conditions upon the terms in the continued
-fraction so that there is redundancy, the conditions are still strong
-enough to ensure effective convergence of the rational approximations.
-
-We now consider ways to implement arithmetic on continued fractions.
-
-\section{The Algebraic Algorithm}
-\label{sc:aa}
-
-Having seen various sorts of continued fractions, and how to turn them
-into rational approximations, we move on to consider the
-implementation of ordinary arithmetic operations on them.  The
-simplest is the algebraic algorithm \mbox{\tt algebraicAlgorithm}\indexTT{algeb
-raicAlgorithm}, which we
-investigate in this section.
-
-The algebraic algorithm satisfies the following condition:
-\[\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}n_0&n_1\\d_0&d_1\end{array}\right)\, r =
-\frac{n_0r+n_1}{d_0r+d_1}\]
-
-That is: \mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm} is a functio
-n that takes a state or
-homography (consisting of four integers) and a real number
-(represented by a continued fraction) and returns the real number
-$(n_0r+n_1)/(d_0r+d_1)$ (again represented as a continued fraction).
-
-Let us briefly look at the domain of definition for this
-algorithm. The first constraint is that if the determinant of the
-homography be non-zero; if the determinant {\em is} zero then the
-result is independent of the value of $r$. This check need only be
-applied once, as the state transformations do not affect the
-determinant of the homography.
-
-The main use of the \mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm} i
-s to perform operations on
-real numbers, such as adding a rational to a real or multiplying a
-rational by a real. Adding the rational $n/d$ is achieved by using the
-following initial configuration:
-\[\frac{n}{d}+r = \mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}d&n\\0&d\end{array}\right)\, r,\]
-and multiplication by $n/d$ is given by:
-\[\frac{n}{d}\times r = \mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorith
-m}\,
-\left(\begin{array}{cc}n&0\\0&d\end{array}\right)\, r.\]
-We observe that the following properties arise with the continued
-fraction representation of the number $r$ and the result. Firstly if
-the continued fraction is finite, then at some stage we will need to
-consider the case where $r$ is infinite; we define the algebraic
-algorithm in this case in the following manner\footnote{The fraction
-$n_0/d_0$ will need to be converted to a continued fraction -- but this is
-easy.}:
-\[\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}n_0&n_1\\d_0&d_1\end{array}\right)\,\infty =
-\frac{n_0}{d_0}.\]
-
-Suppose that we are applying the function
-\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm} to a number of the fo
-rm $x_0+\frac{1}{r}$. Then it is
-easy to verify the following identity:
-\begin{eqnarray*}
-\lefteqn{\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}n_0&n_1\\d_0&d_1\end{array}\right)\,(x_0+\frac{1}{r})} \
-\
-& = & \mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}n_0x_0+n_1&n_0\\d_0x_0+d_1&d_0\end{array}\right)\,r.
-\end{eqnarray*}
-
-Another identity that holds for the function \mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm} is that for any
-$y_0$:
-
-\begin{eqnarray*}
-\lefteqn{
-\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,\left(\begin{array}{
-cc}n_0&n_1\\d_0&d_1\end{array}\right)\,r}\\
-&= & y_0 + \frac{1}{\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}d_0&d_1\\n_0-d_0y_0&n_1-d_1y_0\end{array}\right)\,r }.
-\end{eqnarray*}
-
-These two identities define the {\em consumption} and {\em emission}
-of terms in the continued fraction representations of the initial real
-number $r$ or the resultant real number. The critical question is
-therefore: which of these operations should we choose to perform,
-given a choice? The answer is straightforward if we wish to be as lazy
-as possible: we emit a term provided it is sufficiently close to the
-value we require for an {\bf E}-fraction, otherwise we consume a term
-from the input.
-
-Let's be a bit more explicit: to evaluate the following state
-\[r =\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}n_0&n_1\\d_0&d_1\end{array}\right)\,(x+\frac{1}{x'}),\]
-the Algebraic Algorithm goes through the following steps.
-
-\begin{enumerate}
-\item After checking that the input continued fraction hasn't terminated,
-we use the first term $x$ to determine the bounds for $x+1/x'$, using
-the function \mbox{\tt cfBound}\indexTT{cfBound}.
-\item From this we can determine the bounds on the value of $r$, by
-applying the homography.
-\item If the bounds $(i,\,s)$ on $r$ are such that \mbox{\tt cfBound\ o}
-contains the interval $(i,\,s)$, then we may output \mbox{\tt o} as the next
-term in the continued fraction.  Alternatively, we absorb the term
-$x$.
-\end{enumerate}
-
-\subsection{Generating Output}
-
-One alternative use of the algebraic algorithm is to generate output;
-here we replace the test that determines whether we have sufficient
-accuracy to emit a term, by one that determines whether we have
-calculated the answer to sufficient accuracy yet. The alternative is
-to continue absorbing terms.
-\begin{verbatim}
-> outputAlgorithm homography@([n0,n1],[d0,d1]) [] = n0 %% d0
-> outputAlgorithm homography@([n0,n1],[d0,d1]) xs'@(x:xs)
->  = if i <= s && s-i < accuracy
->    then algebraicEval homography [x%%1]
->    else outputAlgorithm ([n0*x+n1,n0],[d0*x+d1,d0]) xs
->    where (i,s) = algebraicBound homography [x%%1]
-\end{verbatim}
-\indexDTT{outputAlgorithm}%
-\indexDTT{outputAlgorithm}%
-The main use to which we will put the Algebraic Algorithm is as a base
-case for the more useful Quadratic Algorithm, which we now consider.
-
-\section{Quadratic Algorithm}
-\label{sc:qa}
-
-The problem at the moment is that we are limited to rational
-operations on any particular real number we care to consider. What we
-would really like to do is to add and multiply two real numbers
-represented as continued fractions.
-
-This is what Gosper's Quadratic Algorithm does. His original papers
-\shortcite{GOSPER80a,GOSPER81a} are unfortunately a little difficult to
-obtain but a pr\'{e}cis is available in \cite{PEYTON-JONES84b}. This
-paper has as an added bonus, an implementation of continued fractions
-in SASL. The Quadratic Algorithm calculates the following function
-(which, as we shall see, can be used to implement the basic arithmetic
-operations):
-
-\[\mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}\,
-\left(\begin{array}{cccc}n_0&n_1&n_2&n_3\\d_0&d_1&d_2&d_3\end{array}\right)\, x
-\,y = 
-\frac{(n_0x+n_1)y+(n_2x+n_3)}{(d_0x+d_1)y+(d_2x+d_3)}.\]
-Notice that this time the state consists of eight integers and that
-the function \mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm} takes tw
-o real numbers (represented as continued
-fractions) returning another real number.
-
-We first observe that the following base cases can be given:
-\begin{eqnarray*}
-\lefteqn{\mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}\,
-\left(\begin{array}{cccc}n_0&n_1&n_2&n_3\\d_0&d_1&d_2&d_3\end{array}\right)\,\i
-nfty\,y}\\
- &= &
-\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}n_0&n_2\\d_0&d_2\end{array}\right)\,y,
-\end{eqnarray*}
-and
-\begin{eqnarray*}
-\lefteqn{\mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}\,
-\left(\begin{array}{cccc}n_0&n_1&n_2&n_3\\d_0&d_1&d_2&d_3\end{array}\right)\,x\
-,\infty}\\
-& =&
-\mbox{\tt algebraicAlgorithm}\indexTT{algebraicAlgorithm}\,
-\left(\begin{array}{cc}n_0&n_1\\d_0&d_1\end{array}\right)\,x.
-\end{eqnarray*}
-
-As we would expect, we can derive similar identities to those deduced
-for the algebraic algorithm. Firstly, absorbing a term from each
-continued fraction\footnote{To be as lazy as possible, we should
-follow Gosper\shortcite{GOSPER81a} and only absorb a term from one of the
-continued fractions at a time. The idea is to select the one which
-will cause the most rapid convergence.}, we observe that the following
-holds:
-\begin{eqnarray*}
-\lefteqn{\mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}\,
-\left(\begin{array}{cccc}n_0&n_1&n_2&n_3\\d_0&d_1&d_2&d_3\end{array}\right)
-\,(x_0+\frac{1}{x})
-\,(y_0+\frac{1}{y})=}\\
-\lefteqn{\mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}}\\
-&&\left(\begin{array}{cccc}
-n_0x_0y_0+n_1y_0+n_2x_0+n_3&n_0y_0+n_2&n_0x_0+n_1&n_0\\
-d_0x_0y_0+d_1y_0+d_2x_0+d_3&d_0y_0+d_2&d_0x_0+d_1&d_0
-\end{array}\right)\,x\,y.\end{eqnarray*}
-
-Emitting a term, $y_0$, remains essentially as described for the Algebraic
-Algorithm:
-\begin{eqnarray*}
-\lefteqn{\mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}\,
-\left(\begin{array}{cccc}n_0&n_1&n_2&n_3\\d_0&d_1&d_2&d_3\end{array}\right)
-\,x\,y =}\\
-&&z_0 + \frac{1}{\mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}\,\mbox{\tt hom}\indexTT{hom}\,x\,y}\\
-\lefteqn{\mbox{where}}\\
-\mbox{\tt hom}\indexTT{hom} & = &
- \left(\begin{array}{cccc}d_0&d_1&d_2&d_3\\n_0-d_0z_0&n_1-d_1z_0&n_2-d_2z_0&n_3
--d_3z_0
-             \end{array}\right)
-\end{eqnarray*}
-
-To use the Quadratic Algorithm to add two numbers $x$ and $y$ we
-observe that
-\[x+y = \mbox{\tt quadraticAlgorithm}\indexTT{quadraticAlgorithm}\,
-\left(\begin{array}{cccc}0&1&1&0\\0&0&0&1\end{array}\right)
-\,x\,y\]
-Similar initial states may be computed for the other arithmetic
-operations. We next consider the problem of transcendental functions.
-\section{Transcendental Functions}
-\label{sc:tf}
-
-If we are not able to compute logarithms, exponentials, and
-trigonometric functions then our real number implementation is
-severely limited.  Fortunately these functions can be represented as
-continued fractions, as has long been known
-\cite{GAUSS1812a,GAUSS1876a}. In this section we show how to use
-Gauss' continued fraction representation of $F(\alpha , 1, \gamma ;
-x)$ to implement the basic transcendental functions: $\exp$, $\log$,
-$\tan$, and $\arctan$.
-
-As a reminder, $F(\alpha , \beta , \gamma ; x)$ is the following series:
-\[
-[1 + \frac{\alpha \beta}{1! \gamma}x +
-\frac{\alpha(\alpha +1) \beta (\beta +1)}{2! \gamma(\gamma +1)}x^2 +
-\frac{\alpha(\alpha +1)(\alpha +2) \beta (\beta +1)(\beta +2)}%
-{3! \gamma(\gamma +1)(\gamma +2)}x^3 +
-\ldots \]
-
-Modifying Gauss' formulation of the continued fraction we get that
-$F(\alpha , 1, \gamma ; x)$ is:
-
-\begin{eqnarray*}
- [0,~1,~\frac{1}{\alpha}(\frac{\gamma}{-x}),
-~\frac{\alpha}{\gamma-\alpha}(\gamma+1),~
-p_0(\frac{\gamma +2}{-x}),~p'_0(\gamma +3)~\ldots~ &&\\
-\qquad p_{n}(\frac{\gamma +2n+2}{-x}),~ p'_{n}(\gamma +2n+3),~ \ldots ]&&
-\end{eqnarray*}
-where
-\[p_n = \frac{1}{\alpha}\prod_{i=0}^n \frac{(\gamma - \alpha +i)(i+1)}%
-{(\alpha +i+1)(\gamma +i)}\]
-and
-\[p'_n = \frac{\alpha}{\gamma-\alpha}
-\prod_{i=0}^n \frac{(\alpha +i+1)(\gamma +i)}%
-{(\gamma - \alpha +i+1)(i+2)}\]
-
-Fortunately, the $p_n$ and $p_n'$ coefficients often simplify.
-
-\subsection{Exponential}
-
-Gauss gives the following representation of $\exp(x)$:
-$\lim_{k\to\infty}F(k,1,1;\frac{x}{k})$. If we expand this we have
-the following continued fraction:
-
-\[ [0,~1,~1(\frac{-1}{x}),-2,~3(\frac{1}{x}),~\frac{4}{2},~
-5(\frac{-1}{x}),~\frac{-6}{3},~\ldots~]\]
-We then observe that $\exp(x) = \frac{1}{\exp(-x)}$ getting:
-\[ [1,~\frac{1}{x},~-2,~\frac{-3}{x},~2,~\frac{5}{x},~-2,~\ldots ] \]
-
-Setting $x$ to $1$ we have the following continued fraction for $e$:
-\begin{verbatim}
-> exp1CF :: [QRational]
-> exp1CF = 1: f 0
->          where f n = (4*n+1): (-2): (-(4*n+3)): 2: f (n+1)
-\end{verbatim}
-\indexDTT{exp1CF}%
-\indexDTT{exp1CF}%
-To use this continued fraction to generate the continued fraction for
-$\exp(x)$, where $x\in {\cal Q}$, we first test to see whether
-$\mbox{\tt abs}\indexTT{abs}\,\mbox{\tt x}< 1$.  If it is then we multiply each
- odd term of \mbox{\tt exp1CF}\indexTT{exp1CF}
-by $1/x$ which results in each term of the continued fraction (apart
-from the second one) will have an absolute magnitude greater than $2$.
-Alternatively we can multiply the even terms by $x$ and divide the
-result by $x$. This is marginally more stable than Vuillemin's
-formulation.
-\begin{verbatim}
-> expQ  :: QRational -> ContinuedFraction
-> expQ  x
->  = aaQ hom ts2
->    where (ts1,ts2)   = splitAt 2 cf
->          hom         = absorbQ startH ts1
->          (startH,cf) = if abs x < 1 then
->                         (([1,0],[0,1]), oddQ x exp1CF) else
->                         (([qDenominator x,0],
->                           [0,qNumerator x]),
->                                    evenQ (1/x) exp1CF) 
-\end{verbatim}
-\indexDTT{expQ}%
-\indexDTT{expQ}%
-The \mbox{\tt aaQ}\indexTT{aaQ} function is a modified version of the \mbox{\tt
- algebraicAlgorithm}\indexTT{algebraicAlgorithm}
-that accepts continued fractions with rational terms; it still outputs
-integer-termed continued fractions. The \mbox{\tt absorbQ}\indexTT{absorbQ} fun
-ction absorbs
-terms of the finite continued fraction, to create a starting
-homography.
-
-We use an exactly similar process for the case of continued fraction
-arguments. This time, however, we are simply concerned with whether
-$x$ might be $0$. If it is, then calculating $1/x$ will likely lead to
-disaster, and is therefore avoided.
-\begin{verbatim}
-> expR :: ContinuedFraction -> ContinuedFraction
-> expR x = if mightBeZero x
->          then processQS x (splitAt 6 (oddR exp1CF))
->          else quadraticAlgorithm ([0,1,1,0],[0,1,0,0]) x cf2
->          where cf2 = processQS x (splitAt 6
->                                       (oddR (tail exp1CF)))
-\end{verbatim}
-\indexDTT{expR}%
-\indexDTT{expR}%
-\subsection{Logarithm}
-
-Vuillemin's continued fraction is too unstable for values of $x-1$
-which aren't near $0$. So instead we use Gauss' formula:
-\[\log(\frac{1+x}{1-x}) = 2xF(\frac{1}{2},1,\frac{3}{2}; x^2)\]
-It works because the transformation eliminates those values of $x$
-in the above expansion that would cause problems, {\it e.g.} $|x|
-\geq 1$.
-
-The continued fraction for $F(\frac{1}{2},1,\frac{3}{2}; x^2)$ is
-closely related to that for $\arctan(x)$, because:
-\[\tanh^{-1}(x) = \frac{1}{2} \log(\frac{1+x}{1-x})\]
-and
-\[\imath \tanh(x) = \tan(\imath x).\]
-\begin{verbatim}
-> logQ :: QRational -> ContinuedFraction
-> logQ x = aaQ startH ts2
->          where (ts1,ts2) = splitAt 2 (oddQ (-t*t) atan1CF)
->                t         = (x-1)/(x+1)
->                tn        = qNumerator t
->                td        = qDenominator t
->                startH    = absorbQ ([-2*td,0],[0,tn]) ts1
-\end{verbatim}
-\indexDTT{logQ}%
-\indexDTT{logQ}%
-The nice part is that by requiring the multiplier \mbox{\tt -t*t} to be less
-than \mbox{\tt 1}, we place no constraints on \mbox{\tt x} provided that it is 
-already
-in the range $(0,\infty )$. Because there are no range decisions to be
-taken, the continued fraction argument function \mbox{\tt logR}\indexTT{logR} i
-s analogous to
-that of \mbox{\tt logQ}\indexTT{logQ}.
-
-\subsection{Tangent}
-
-For this function we use Lambert's Continued Fraction:
-\[\tan(x) = [0,~\frac{1}{x},~\frac{-3}{x},~\frac{5}{x},~\frac{-7}{x},~\ldots]\]
-
-It is much more convenient if we can obtain a fraction with some
-constant terms in it so we observe the following transformation is
-possible:
-\[\tan(x) = \frac{1}{x}[0,~\frac{1}{x^2},~-3,~\frac{5}{x^2},~-7,~\ldots]\]
-provided that $x$ is non-zero.
-\begin{verbatim}
-> tan1CF :: [QRational]
-> tan1CF = 0: f 0 where f n = (4*n+1) : (-(4*n+3)) : f (n+1)
-\end{verbatim}
-\indexDTT{tan1CF}%
-\indexDTT{tan1CF}%
-\begin{verbatim}
-> tanQ  :: QRational -> ContinuedFraction
-> tanQ  x = if x /= 0 then aaQ startH ts2
->                     else error "tanQ: 0 argument"
->           where n         = qRound (abs (x/2)) + 1
->                 (ts1,ts2) = splitAt n (bothQ x tan1CF)
->                 startH    = absorbQ ([1,0],[0,1]) ts1
-\end{verbatim}
-\indexDTT{tanQ}%
-\indexDTT{tanQ}%
-The case for continued fraction arguments is \mbox{\tt tanR}\indexTT{tanR}, and
- works
-similarly; the cotangent function \mbox{\tt cotR}\indexTT{cotR} merely inverts 
-the top-level
-division.
-
-\subsection{Arctangent}
-
-Recall, from the discussion of the continued fraction for $\log$, that
-the $\arctan(x)$ is $xF(\frac{1}{2},1,\frac{3}{2};-x^2)$ 
-The terms in this continued fraction are tending towards $\pi$ and $1$.
-\begin{verbatim}
-> atan1CF :: [QRational]
-> atan1CF
->   = 0: 1: 3: 5/4: f (8/9) 0
->     where f :: QRational -> QRational -> [QRational]
->           f p n = (7/2 + 2*n)*p: (9/2 + 2*n)*p': f p'' (n+1)
->                   where p'  = 1/(p*(n+2)*(n+2))
->                         p'' = 1/(p'*(n+5/2)*(n+5/2))
-\end{verbatim}
-\indexDTT{atan1CF}%
-The \mbox{\tt atanQ}\indexTT{atanQ} function is valid for $x \in (-\infty,\,\infty)$.
-
-We note the following identities:
-\[\mbox{for $x>0$ } \arctan(x) = \frac{\pi}{2} -\arctan(\frac{1}{x})\]
-and
-\[\mbox{for $x<0$ } \arctan(x) = -\frac{\pi}{2} -\arctan(\frac{1}{x})\]
-
-\begin{verbatim}
-> atanQ :: QRational -> ContinuedFraction
-> atanQ x
->  = if x > 1  then
->       quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1])
->                          atan1 (cf (1/x))
->    else if x < -1 then
->       quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1])
->                          atan1 (cf (1/x))
->    else  cf x
->    where
->    cf x = aaQ startH ts2
->           where
->           (ts1,ts2) = splitAt 4 (oddQ (x*x) atan1CF)
->           startH    = absorbQ ([qDenominator x,0],
->                                [0,qNumerator x]) ts1
-\end{verbatim}
-\indexDTT{atanQ}%
-\indexDTT{atanQ}%
-Because we will be using it to help improve the stability of the
-$arctan$ function \mbox{\tt atanR}\indexTT{atanR}, we pre-define the constant c
-ontinued
-fraction $\mbox{\tt atan1}\indexTT{atan1} = \pi/4$.
-\begin{verbatim}
-> atan1 :: ContinuedFraction
-> atan1 = atanQ (1%%1)
-\end{verbatim}
-\indexDTT{atan1}%
-\indexDTT{atan1}%
-As in the case of \mbox{\tt expR}\indexTT{expR} -- because we are now distingui
-shing between
-different domains of definition -- we inspect the leading term of the
-argument continued fraction.
-\begin{verbatim}
-> atanR :: ContinuedFraction -> ContinuedFraction
-> atanR []         = algebraicAlgorithm ([2,0],[0,1]) atan1
-> atanR xs@(x:xs')
->  = if x > 1 then
->       quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1])
->                          atan1 ys else
->    if x < -1 then
->       quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1])
->                          atan1 ys else
->       atanRdirect xs
->    where ys = atanRdirect (algebraicAlgorithm ([0,1],[1,0]) xs)
-\end{verbatim}
-\indexDTT{atanR}%
-\indexDTT{atanR}%
-\indexDTT{atanR}%
-The function \mbox{\tt atanRdirect}\indexTT{atanRdirect} deals with the straigh
-tforward case where
-the argument is within the domain of definition.
-
-\subsection{Processing Irrational Transcendentals}
-
-To produce continued fractions for the transcendental functions with
-continued fraction arguments we use an ``infinite state machine''.
-This infinite state is an infinite list of sub-states, each consisting
-of a rational, and a homography for the Quadratic Algorithm.
-\begin{verbatim}
-> type State = (QRational, Homography)
-\end{verbatim}
-\indexDTT{type}%
-To set up the initial continued fraction for the ``infinite state
-machine'' from the initial continued fraction we use \mbox{\tt oddR}\indexTT{od
-dR}.
-If the initial continued fraction is
-\[[\frac{a_0}{b_0},~\frac{a_1}{b_1},~\frac{a_2}{b_2},~\frac{a_3}{b_3},~\ldots]\]
-then the ``infinite state'' will be:
-\[[\left(\frac{a_0}{b_0},
-   \left(\begin{array}[c]{cccc}0&a_1&b_1&0\\b_1&0&0&0\end{array}\right),~
-\right)
-\left(\frac{a_2}{b_2},
-   \left(\begin{array}[c]{cccc}0&a_3&b_3&0\\b_3&0&0&0\end{array}\right)
-\right),~
-\ldots]\]
-
-The processing of the infinite state and the continued fraction is
-initiated by \mbox{\tt processQS}\indexTT{processQS}. What it does is to absorb
- leading terms of the state
-
-\begin{verbatim}
-> processQS :: ContinuedFraction -> ([State],[State])
->                                -> ContinuedFraction
-> processQS x (hs1,hs2) = preProcess x hs1 (qsSetupAlt hs2 x)
-\end{verbatim}
-\indexDTT{processQS}%
-\indexDTT{processQS}%
-This in turn makes use of \mbox{\tt preProcess}\indexTT{preProcess} which is us
-ed to absorb some of
-the leading terms (\mbox{\tt hs1}\indexTT{hs1}) and \mbox{\tt qsSetupAlt}\indexTT{qsSetupAlt} which is used to process
-the remainder (\mbox{\tt hs2}\indexTT{hs2}).
-\begin{verbatim}
-> preProcess :: ContinuedFraction -> -- the argument CF
->               [State]           -> -- leading parts of inf state
->               ContinuedFraction -> -- accumulating CF
->               ContinuedFraction    -- the resulting CF
-> preProcess _ []         c = c
-> preProcess x ((i,h):hs) c
->  = algebraicAlgorithm ([ni,nd],[nd,0]) 
->                       (quadraticAlgorithm h x (preProcess x hs c))
->    where ni = qNumerator i
->          nd = qDenominator i
-\end{verbatim}
-\indexDTT{preProcess}%
-\indexDTT{preProcess}%
-\indexDTT{preProcess}%
-To process the remainder of the infinite state, we use \mbox{\tt qsSetupAlt}\indexTT{qsSetupAlt}.
-\begin{verbatim}
-> qsSetupAlt :: [State] -> ContinuedFraction
->                       -> ContinuedFraction
-> qsSetupAlt hs x = aaQ ([1,0],[0,1]) (qsProcess (hs, x))
-\end{verbatim}
-\indexDTT{qsSetupAlt}%
-\indexDTT{qsSetupAlt}%
-\subsection{Quadratic Process}
-
-In essence the state for the algorithm holds:
-\begin{enumerate}
-\item the infinite list of \mbox{\tt Rational}\indexTT{Rational}/\mbox{\tt Homo
-graphy}\indexTT{Homography} pairs; and
-\item the infinite continued fraction of the argument.
-\end{enumerate}
-
-We advance by absorbing one term from the continued fraction or by
-considering one more term of the infinite state.
-\begin{verbatim}
-> qsProcess :: ([State], ContinuedFraction) -> [QRational]
-> qsProcess    (ihs, [])    = asProcess ihs
-> qsProcess st@(ihs, (_:_)) = qsNext 1 st
-\end{verbatim}
-\indexDTT{qsProcess}%
-\indexDTT{qsProcess}%
-\indexDTT{qsProcess}%
-The decision about whether sufficient accuracy has been achieved is
-decided by \mbox{\tt qsNext}\indexTT{qsNext}. To do this a call to \mbox{\tt qs
-Term}\indexTT{qsTerm} is made. This
-returns a list of output values and a modified state (\mbox{\tt ihs'}).  If the
-output list is sufficiently long we can output them as part of the
-resulting continued fraction. Alternatively we can absorb another term
-of the argument continued fraction. In either case we consider one
-more term of the infinite state.
-\begin{verbatim}
-> qsNext :: Int -> ([State], ContinuedFraction) -> [QRational]
-> qsNext n (ihs, []) = asProcess (qsEnd ihs)
-> qsNext n st@(((i,h):ihs), xs'@(x:xs))
->  = if length outs > 1
->       then init outs ++ qsNext (n+1) (ihs',xs')
->       else qsNext (n+1) (qsIn x ((i,h): ihs), xs)
->    where (outs, (ihs',_)) = qsTerm n st
-\end{verbatim}
-\indexDTT{qsNext}%
-\indexDTT{qsNext}%
-\indexDTT{qsNext}%
-The \mbox{\tt qsTerm}\indexTT{qsTerm} function considers the initial \mbox{\tt 
-n} terms of the infinite
-state, and the first term of the continued fraction for the argument.
->From this it calculates a modified infinite state, in which, what is
-effectively carry propogation has occurred.
-\begin{verbatim}
-> qsTerm :: Int -> ([State], ContinuedFraction)
->          -> ([QRational],([State], ContinuedFraction))
-> qsTerm n (ihs,xs)
->  = (outs, (ihs',xs))
->    where (front,back@((i,_):_)) = splitAt n ihs
->          (outs,ihs') = foldr (qsTerm' (head xs))
->                                       ([i],back) front
-\end{verbatim}
-\indexDTT{qsTerm}%
-\indexDTT{qsTerm}%
-To deal with each individual term of the infinite state we use \mbox{\tt qsTerm
-'}.
-\begin{verbatim}
-> qsTerm' :: Integer -> State -> ([QRational], [State]) -> 
->                                ([QRational], [State])
-> qsTerm' x ih@(i,h) (is,ihs)
->  = (is',((last is',h''):ihs))
->    where h' = qsAbsorb (init is) h
->          i' = last is
->          (is',h'') = qsEmit x i' ([i],h')
-\end{verbatim}
-This in turn uses \mbox{\tt qsAbsorb}\indexTT{qsAbsorb} and \mbox{\tt qsEmit}\i
-ndexTT{qsEmit} to propogate the carry and
-emit output terms, respectively.
-\begin{verbatim}
-> qsAbsorb :: [QRational] -> Homography -> Homography
-> qsAbsorb [] h = h
-> qsAbsorb (i:is) h@([n0,n1,n2,n3],[d0,d1,d2,d3])
->  = qsAbsorb is ([n0*n+n2*d, n1*n+n3*d, n0*d, n1*d],
->                 [d0*n+d2*d, d1*n+d3*d, d0*d, d1*d])
->    where n = qNumerator i
->          d = qDenominator i
-\end{verbatim}
-\indexDTT{qsAbsorb}%
-\indexDTT{qsAbsorb}%
-\indexDTT{qsAbsorb}%
-\mbox{\tt qsEmit}\indexTT{qsEmit} is essentially a variant of the Quadratic Alg
-orithm, capable
-of dealing with the output of more than $1$ term at once.
-\begin{verbatim}
-> qsEmit :: Integer -> QRational -> ([QRational], Homography)
->                                -> ([QRational], Homography)
-> qsEmit x i ish@(is,h@([n0,n1,n2,n3],[d0,d1,d2,d3]))
->  = case quadraticOutput h [x%%1,i] of
->      Just o  -> qsEmit x i (is++[o%%1],
->                             ([d0,d1,d2,d3],
->                              [n0-d0*o, n1-d1*o,
->                                   n2-d2*o, n3-d3*o]))
->      Nothing -> ish
-\end{verbatim}
-\indexDTT{qsEmit}%
-\indexDTT{qsEmit}%
-To absorb the leading term of the argument continued fraction into
-each Quadratic Homography we use \mbox{\tt qsIn}\indexTT{qsIn}.
-\begin{verbatim}
-> qsIn :: Integer -> [State] -> [State]
-> qsIn = map . input
->        where input x (i,([n0,n1,n2,n3],[d0,d1,d2,d3]))
->                = (i, ([n0*x+n1, n0, n2*x+n3, n2],
->                       [d0*x+d1, d0, d2*x+d3, d2]))
-\end{verbatim}
-\indexDTT{qsIn}%
-\indexDTT{qsIn}%
-\subsection{Finite Continued Fractions}
-
-It is possible to generate finite continued fractions, so we need a
-way to deal with \mbox{\tt []} in the \mbox{\tt qsProcess}\indexTT{qsProcess}. 
-This results in a
-degeneration to the algebraic process.
-
-\mbox{\tt qsEnd}\indexTT{qsEnd} turns each Quadratic Homography into the equiva
-lent Algebraic
-Homography.
-\begin{verbatim}
-> qsEnd :: [State] -> [State]
-> qsEnd = map input
->         where input (i,([n0,_,n2,_],[d0,_,d2,_]))
->                 = (i,([n0,n2],[d0,d2]))
-\end{verbatim}
-\indexDTT{qsEnd}%
-\indexDTT{qsEnd}%
-The function \mbox{\tt asProcess}\indexTT{asProcess} munges these homographies 
-to produce an answer.
-\begin{verbatim}
-> asProcess :: [State] -> [QRational]
-> asProcess ihs = fst (head ihs) : asProcess (asNext ihs)
-\end{verbatim}
-\indexDTT{asProcess}%
-\indexDTT{asProcess}%
-\begin{verbatim}
-> asNext :: [State] -> [State]
-> asNext ((i,h@([n0,n1],[d0,d1])):ihs'@((i',_):ihs))
->  = case algebraicOutput h i' of
->      Just o  -> (o%%1, ([d0,d1], [n0-d0*o,n1-d1*o])): ihs'
->      Nothing -> asNext ((i,hIn i' h): asNext ihs')
->    where hIn x ([n0,n1],[d0,d1]) = ([n0*xn+n1*xd,n0*xd],
->                                     [d0*xn+d1*xd,d0*xd])
->                                    where xn = qNumerator x
->                                          xd = qDenominator x
-\end{verbatim}
-\indexDTT{asNext}%
-\indexDTT{asNext}%
-And that completes the implementation of the specified subset of
-transcendental functions.
-
-\section{Conclusion}
-
-The system only computes the answer to some arithmetic expression to
-whatever accuracy you specify. It may run out of memory (in which case
-the program will inform you of this fact); it may fail to compute any
-answer whatsoever (in which case you will wait for ever to observe
-the output of the program). It does have major advantages over
-traditional floating point numbers.  Firstly, you may be sure that any
-answer it provides is accurate. Secondly, there is no need to re-code
-an algorithm when more accuracy is required.
-
-\section{Acknowledgements}
-
-Jean Vuillemin's work \shortcite{VUILLEMIN87a} is clearly the major
-source for the implementation work that I have undertaken. I would
-also like to thank Geoffrey Burn, Simon Peyton Jones, John Robson,
-Peter Sestoft, and Geraint Jones who have all contributed comments
-upon preliminary versions of this paper. Of course, all errors and
-indiscretions remain my own work.
-
-\begin{thebibliography}{}
-
-\bibitem[\protect\citename{Bishop and Bridges, }1985]{BISHOP85a}
-Bishop,~E. and Bridges,~D. 1985.
-\newblock {\em Constructive Analysis}.
-\newblock Springer Verlag.
-
-\bibitem[\protect\citename{Boehm and Cartwright, }1990]{Boehm90a}
-Boehm,~H. and Cartwright,~R. 1990.
-\newblock Exact Real Arithmetic: Formulating Real Numbers as Functions
-\newblock In {\em Research Topics in Functional Programming}, pages 43--64.
-  Addison-Wesley.
-
-\bibitem[\protect\citename{Church, }1941]{CHURCH41a}
-Church,~A. 1941.
-\newblock {\em The Calculi of Lambda-Conversion}.
-\newblock Princeton University Press, Princeton, N.J.
-
-\bibitem[\protect\citename{Euclid}]{EUCLID}
-Euclid.
-\newblock Elements.
-
-\bibitem[\protect\citename{Gauss, }1812]{GAUSS1812a}
-Gauss,~C.F. 1812.
-\newblock Disquisitiones generales circa serium infinitam
-  $1+\frac{\alpha\beta}{1.\gamma}x + $ etc. pars prior.
-\newblock In {\em Werke}, volume~3, pages 123--162. K\"{o}niglichen
-  Gesellschaft der Wissenschaften, G\"{o}ttingen.
-
-\bibitem[\protect\citename{Gauss }1876]{GAUSS1876a}
-Gauss,~C.F. 1876.
-\newblock {\em Werke}, volume~3.
-\newblock K\"{o}niglichen Gesellschaft der Wissenschaften, G\"{o}ttingen.
-
-\bibitem[\protect\citename{Gosper }1980]{GOSPER80a}
-Gosper,~R.W. 1980.
-\newblock Continued fraction arithmetic.
-\newblock HAKMEM 101b, MIT.
-
-\bibitem[\protect\citename{Gosper }1981]{GOSPER81a}
-Gosper,~R.W. 1981.
-\newblock Continued fraction arithmetic.
-\newblock Unpublished Draft Paper.
-
-\bibitem[\protect\citename{Peyton Jones }1984]{PEYTON-JONES84b}
-{Peyton Jones},~S.L. 1984.
-\newblock Arbitrary precision arithmetic using continued fractions.
-\newblock Internal Note 1530, Deptartment of Computer Science, University
-  College London.
-
-\bibitem[\protect\citename{Pour El and Richards }1989]{POUR-EL89a}
-{Pour-El},~M.B. and Richards,~J.I. 1989.
-\newblock {\em Computability in Analysis and Physics}.
-\newblock Perspectives in Mathematical Logic. Springer-Verlag, Heidelberg.
-
-\bibitem[\protect\citename{Rice }1954]{RICE54a}
-Rice,~H.G. 1954.
-\newblock Recursive real numbers.
-\newblock {\em Proceedings of the American Mathematical Society},
-  5(5):784--791.
-
-\bibitem[\protect\citename{Turing }1936]{TURING36a}
-Turing,~A.M. 1936.
-\newblock On computable numbers, with an application to the
-  {E}ntscheidungsproblem.
-\newblock {\em Proceedings of the London Mathematical Society}, 42(3):230--265.
-
-\bibitem[\protect\citename{Vuillemin }1987]{VUILLEMIN87a}
-Vuillemin,~J. 1987.
-\newblock Arithm\'{e}tic r\'{e}elle exacte par les fractions continues.
-\newblock Technical Report 760, Instituit National de Recherche en Informatique
-  et en Automatique, Domaine de Voluceau, Roquencourt, BP105, 78153 Le Chesnay
-  Cedex, France.
-
-\end{thebibliography}
-\end{document}