1 module RealM (RealT, evalReal, int2Real, addReal, subReal, mulReal, 2 divReal, sqrtReal) where 3 4 data Tree a = TreeC a (Tree a) (Tree a) 5 6 data RealT = RealC (Tree Integer) deriving () 7 8 9 treeFrom :: Integer -> Tree Integer 10 treeFrom n = TreeC n (treeFrom (2*n-1)) (treeFrom (2*n-2)) 11 ------------------------------------------------------------------------------- 12 13 14 numberTree :: Tree Integer 15 numberTree = treeFrom 0 16 ------------------------------------------------------------------------------- 17 18 19 -- multiply an integer, x, by 10^(-n). 20 scale :: Integer -> Integer -> Integer 21 scale x n | n <= 0 = (x * 10^(1-n) + 5) `quot` 10 22 scale x n | n > 0 = (x `quot` 10^(n-1) + 5) `quot` 10 23 ------------------------------------------------------------------------------- 24 25 26 -- return the position of the most significant digit of a real. The second 27 -- argument is an accumulator in which the value is collected. This should be 0 28 -- in the initial call. 29 msd :: Tree Integer -> Integer -> Integer 30 msd x k = if (xk >= 1) && (xk < 10) 31 then k 32 else if (xk == 10) then (k+1) 33 else if xk >= 1 34 then (msd x (k+1)) 35 else (msd x (k-1)) 36 37 where xk = abs (evalReal (RealC x) k) 38 ------------------------------------------------------------------------------- 39 40 41 -- evaluate a real with a tolerance of tol. 42 evalReal :: RealT -> Integer -> Integer 43 evalReal (RealC tree) tol | tol <= 0 = lookupTree tree tol 44 evalReal (RealC tree) tol | tol > 0 = scale (lookupTree tree 0) tol 45 ------------------------------------------------------------------------------- 46 47 48 lookupTree :: Tree Integer -> Integer -> Integer 49 lookupTree tree n = followPath (getPath (1-n) []) tree 50 where 51 getPath 1 p = p 52 getPath n p | even n = getPath (n `quot` 2) (0:p) 53 getPath n p | odd n = getPath (n `quot` 2) (1:p) 54 55 followPath [] (TreeC x _ _) = x 56 followPath (0:ps) (TreeC x lc _) = followPath ps lc 57 followPath (1:ps) (TreeC x _ rc) = followPath ps rc 58 ------------------------------------------------------------------------------- 59 60 61 mapTree :: (a -> b) -> Tree a -> Tree b 62 mapTree f (TreeC x lc rc) = TreeC (f x) (mapTree f lc) (mapTree f rc) 63 ------------------------------------------------------------------------------- 64 65 66 -- convert an integer to a real. 67 int2Real :: Integer -> RealT 68 int2Real x = RealC (mapTree (\n -> scale x n) numberTree) 69 ------------------------------------------------------------------------------- 70 71 72 -- add two real numbers. 73 addReal :: RealT -> RealT -> RealT 74 addReal x y = 75 RealC (mapTree (\n -> ((evalReal x (n-1)) + 76 (evalReal y (n-1)) + 5) `quot` 10) 77 numberTree) 78 ------------------------------------------------------------------------------- 79 80 81 -- subtract two real numbers. 82 subReal :: RealT -> RealT -> RealT 83 subReal x y = 84 RealC (mapTree (\n -> ((evalReal x (n-1)) - 85 (evalReal y (n-1)) + 5) `quot` 10) 86 numberTree) 87 ------------------------------------------------------------------------------- 88 89 90 -- multiply two real numbers. 91 mulReal :: RealT -> RealT -> RealT 92 mulReal (xr@(RealC x)) (yr@(RealC y)) = 93 RealC (mapTree (\n -> let 94 m = if n >= 0 then (n+1) `quot` 2 95 else (n-1) `quot` 2 96 x0 = evalReal xr m 97 y0 = evalReal yr m 98 99 mulAux x y = if y1 == 0 100 then 0 101 else scale (x1 * y1) 102 (4 + (msd x 0) + 103 (msd y 0) - n) 104 105 where 106 y1 = (evalReal yr 107 (n - (msd x 0) - 2)) 108 x1 = (evalReal xr 109 (n - (msd y 0) - 2)) 110 in if (x0 == 0) && (y0 == 0) 111 then 0 112 else if x0 == 0 113 then mulAux y x 114 else mulAux x y) 115 numberTree) 116 ------------------------------------------------------------------------------- 117 118 119 -- divide two real numbers. 120 divReal :: RealT -> RealT -> RealT 121 divReal x y = 122 mulReal x (reciprocalReal y) 123 ------------------------------------------------------------------------------- 124 125 126 -- find the reciprocal of a real number. 127 reciprocalReal :: RealT -> RealT 128 reciprocalReal xr@(RealC x) = 129 RealC (mapTree 130 (\n -> ((f (dx + n)) + 5) `quot` 10) 131 numberTree) 132 where 133 dx = msd x 0 134 135 f m = if m >= 0 136 then 0 137 else if m == -1 138 then 10000 `quot` (evalReal xr (dx - 2)) 139 else (scale ((fm) * ((scale 2 (m + hm - 2)) - 140 (evalReal xr (dx + m - 1)) * (fm))) 141 (2 - 2*hm)) 142 143 where hm = (m-1) `quot` 2 144 fm = f hm 145 ------------------------------------------------------------------------------- 146 147 148 -- find the square of a real number. 149 sqrReal :: RealT -> RealT 150 sqrReal (xr@(RealC x)) = 151 RealC (mapTree (\n -> let 152 m = if n >= 0 then (n+1) `quot` 2 153 else (n-1) `quot` 2 154 x0 = evalReal xr m 155 x1 = evalReal xr (n - (msd x 0) - 2) 156 in if (x0 == 0) 157 then 0 158 else scale (x1 * x1) (4 + 2*(msd x 0) - n)) 159 numberTree) 160 ------------------------------------------------------------------------------- 161 162 163 -- find the square root of a real number. 164 sqrtReal :: RealT -> RealT 165 sqrtReal xr@(RealC x) = 166 RealC (mapTree 167 (\n -> if (evalReal xr (2*n)) == 0 168 then 0 169 else f n (2*hdx+2) xr (-2*hdx-2+n)) 170 numberTree) 171 172 where 173 dx = msd x 0 174 hdx = if dx >= 0 then (dx+1) `quot` 2 else dx `quot` 2 175 176 f n dx xr m = if m >= 0 177 then 0 178 else if m == -1 179 then scale ((round . sqrt . fromInteger) 180 181 (evalReal xr (dx - 10))) 182 (5-hx+n) 183 else (fm + xv `quot` fm) `quot` 2 184 where hm = (m-1) `quot` 2 185 fm = (f n dx xr hm) 186 xv = evalReal xr (n+n) 187 hx = dx `quot` 2 188 189 ------------------------------------------------------------------------------- 190 191 192 {- 193 -- representation for pi. This representation is grossly inefficient. 194 pi :: RealT 195 pi = RealC (mapTree 196 (\n -> if n > 0 then 0 else evalReal f n 197 where 198 f = iter (int2Real 1) 199 (sqrtReal (RealC (mapTree (\n -> scale 5 (n+1)) 200 numberTree))) 201 (RealC (mapTree (\n -> scale 25 (n+2)) numberTree)) 202 (int2Real 1) 203 iter a b t x = if evalReal (subReal a b) (n-1) <= 1 204 then divReal (sqrReal a) t 205 else iter y 206 (sqrtReal (mulReal a b)) 207 (subReal t (mulReal x (sqrReal d))) 208 (mulReal x (int2Real 2)) 209 where y = divReal (addReal a b) (int2Real 2) 210 d = subReal y a) 211 numberTree) 212 ------------------------------------------------------------------------------- 213 -}