Commit 6b234b1c authored by Michal Terepeta's avatar Michal Terepeta Committed by Ben Gamari

spectral: remove compreals

It's been disabled for a long time and there's no input to actually
run it.
Signed-off-by: Michal Terepeta's avatarMichal Terepeta <michal.terepeta@gmail.com>

Test Plan: build & run

Reviewers: bgamari

Differential Revision: https://phabricator.haskell.org/D3330
parent 435f60aa
......@@ -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
\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]
\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}
\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
\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