module RealM (RealT, evalReal, int2Real, addReal, subReal, mulReal,
divReal, sqrtReal) where
data Tree a = TreeC a (Tree a) (Tree a)
data RealT = RealC (Tree Integer) -- deriving ()
treeFrom :: Integer -> Tree Integer
treeFrom n = TreeC n (treeFrom (2*n-1)) (treeFrom (2*n-2))
-------------------------------------------------------------------------------
numberTree :: Tree Integer
numberTree = treeFrom 0
-------------------------------------------------------------------------------
-- multiply an integer, x, by 10^(-n).
scale :: Integer -> Integer -> Integer
scale x n | n <= 0 = (x * 10^(1-n) + 5) `quot` 10
scale x n | n > 0 = (x `quot` 10^(n-1) + 5) `quot` 10
-------------------------------------------------------------------------------
-- return the position of the most significant digit of a real. The second
-- argument is an accumulator in which the value is collected. This should be 0
-- in the initial call.
msd :: Tree Integer -> Integer -> Integer
msd x k = if (xk >= 1) && (xk < 10)
then k
else if (xk == 10) then (k+1)
else if xk >= 1
then (msd x (k+1))
else (msd x (k-1))
where xk = abs (evalReal (RealC x) k)
-------------------------------------------------------------------------------
-- evaluate a real with a tolerance of tol.
evalReal :: RealT -> Integer -> Integer
evalReal (RealC tree) tol | tol <= 0 = lookupTree tree tol
evalReal (RealC tree) tol | tol > 0 = scale (lookupTree tree 0) tol
-------------------------------------------------------------------------------
lookupTree :: Tree Integer -> Integer -> Integer
lookupTree tree n = followPath (getPath (1-n) []) tree
where
getPath 1 p = p
getPath n p | even n = getPath (n `quot` 2) (0:p)
getPath n p | odd n = getPath (n `quot` 2) (1:p)
followPath [] (TreeC x _ _) = x
followPath (0:ps) (TreeC x lc _) = followPath ps lc
followPath (1:ps) (TreeC x _ rc) = followPath ps rc
-------------------------------------------------------------------------------
mapTree :: (a -> b) -> Tree a -> Tree b
mapTree f (TreeC x lc rc) = TreeC (f x) (mapTree f lc) (mapTree f rc)
-------------------------------------------------------------------------------
-- convert an integer to a real.
int2Real :: Integer -> RealT
int2Real x = RealC (mapTree (\n -> scale x n) numberTree)
-------------------------------------------------------------------------------
-- add two real numbers.
addReal :: RealT -> RealT -> RealT
addReal x y =
RealC (mapTree (\n -> ((evalReal x (n-1)) +
(evalReal y (n-1)) + 5) `quot` 10)
numberTree)
-------------------------------------------------------------------------------
-- subtract two real numbers.
subReal :: RealT -> RealT -> RealT
subReal x y =
RealC (mapTree (\n -> ((evalReal x (n-1)) -
(evalReal y (n-1)) + 5) `quot` 10)
numberTree)
-------------------------------------------------------------------------------
-- multiply two real numbers.
mulReal :: RealT -> RealT -> RealT
mulReal (xr@(RealC x)) (yr@(RealC y)) =
RealC (mapTree (\n -> let
m = if n >= 0 then (n+1) `quot` 2
else (n-1) `quot` 2
x0 = evalReal xr m
y0 = evalReal yr m
mulAux x y = if y1 == 0
then 0
else scale (x1 * y1)
(4 + (msd x 0) +
(msd y 0) - n)
where
y1 = (evalReal yr
(n - (msd x 0) - 2))
x1 = (evalReal xr
(n - (msd y 0) - 2))
in if (x0 == 0) && (y0 == 0)
then 0
else if x0 == 0
then mulAux y x
else mulAux x y)
numberTree)
-------------------------------------------------------------------------------
-- divide two real numbers.
divReal :: RealT -> RealT -> RealT
divReal x y =
mulReal x (reciprocalReal y)
-------------------------------------------------------------------------------
-- find the reciprocal of a real number.
reciprocalReal :: RealT -> RealT
reciprocalReal xr@(RealC x) =
RealC (mapTree
(\n -> ((f (dx + n)) + 5) `quot` 10)
numberTree)
where
dx = msd x 0
f m = if m >= 0
then 0
else if m == -1
then 10000 `quot` (evalReal xr (dx - 2))
else (scale ((fm) * ((scale 2 (m + hm - 2)) -
(evalReal xr (dx + m - 1)) * (fm)))
(2 - 2*hm))
where hm = (m-1) `quot` 2
fm = f hm
-------------------------------------------------------------------------------
-- find the square of a real number.
sqrReal :: RealT -> RealT
sqrReal (xr@(RealC x)) =
RealC (mapTree (\n -> let
m = if n >= 0 then (n+1) `quot` 2
else (n-1) `quot` 2
x0 = evalReal xr m
x1 = evalReal xr (n - (msd x 0) - 2)
in if (x0 == 0)
then 0
else scale (x1 * x1) (4 + 2*(msd x 0) - n))
numberTree)
-------------------------------------------------------------------------------
-- find the square root of a real number.
sqrtReal :: RealT -> RealT
sqrtReal xr@(RealC x) =
RealC (mapTree
(\n -> if (evalReal xr (2*n)) == 0
then 0
else f n (2*hdx+2) xr (-2*hdx-2+n))
numberTree)
where
dx = msd x 0
hdx = if dx >= 0 then (dx+1) `quot` 2 else dx `quot` 2
f n dx xr m = if m >= 0
then 0
else if m == -1
then scale ((round . sqrt . fromInteger)
(evalReal xr (dx - 10)))
(5-hx+n)
else (fm + xv `quot` fm) `quot` 2
where hm = (m-1) `quot` 2
fm = (f n dx xr hm)
xv = evalReal xr (n+n)
hx = dx `quot` 2
-------------------------------------------------------------------------------
{-
-- representation for pi. This representation is grossly inefficient.
pi :: RealT
pi = RealC (mapTree
(\n -> if n > 0 then 0 else evalReal f n
where
f = iter (int2Real 1)
(sqrtReal (RealC (mapTree (\n -> scale 5 (n+1))
numberTree)))
(RealC (mapTree (\n -> scale 25 (n+2)) numberTree))
(int2Real 1)
iter a b t x = if evalReal (subReal a b) (n-1) <= 1
then divReal (sqrReal a) t
else iter y
(sqrtReal (mulReal a b))
(subReal t (mulReal x (sqrReal d)))
(mulReal x (int2Real 2))
where y = divReal (addReal a b) (int2Real 2)
d = subReal y a)
numberTree)
-------------------------------------------------------------------------------
-}
|
