1 2 3 4 5 6 7 8 9 10 11 12 13 module ContinuedFractions 14 (module Maybe, module QRationals, 15 ContinuedFraction(..), Homography(..), Interval(..), rat2cf, 16 algebraicAlgorithm, algebraicOutput, 17 quadraticAlgorithm, quadraticOutput, 18 cfRat2CFSqrt, decimals, accuracy, cf2Rat, integerFraction) 19 where 20 import Maybe 21 import List( sort ) 22 import QRationals 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 type ContinuedFraction = [Integer] 41 type Homography = ([Integer],[Integer]) 42 type Interval = (QRational,QRational) 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 rat2cf :: QRational -> ContinuedFraction 59 rat2cf q 60 = if qUndefined q then error "rat2cf: undefined" else 61 if qInfinite q then [] else 62 rat2cf' q 63 where rat2cf' q = x : if q' == 0 then [] else rat2cf' (1 / q') 64 where x = qRound q 65 q' = q - (x%%1) 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 cfBound :: QRational -> Interval 81 cfBound x 82 = if qInfinite x then (1%%0, 1%%0) else 83 if abs x < 2 then (-2%%1, 2%%1) else 84 if x < 0 then (1%%2-x, x+1%%2) else 85 (x-1%%2, -x-1%%2) 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 rootThreeMinus = 1732050807 %% 1000000000 :: QRational 107 rootThreePlus = 1732050808 %% 1000000000 :: QRational 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 generateOutput 123 :: (Homography -> [QRational] -> Interval) -> -- Bound 124 (Homography -> [QRational] -> QRational) -> -- Evaluator 125 Homography -> -- Homography 126 [QRational] -> -- Leading Input Terms 127 Maybe Integer -- Possible Output 128 129 generateOutput bound eval homography qs 130 = if not (qInfinite q0 && qInfinite q1) -- && signum (o%%1) == signum o' 131 && acceptable int (o%%1) then Just o else Nothing 132 where int@(q0,q1) = bound homography qs 133 o' = eval homography qs 134 o = qRound (if abs o' < 2 then o' else 135 if abs q0 < abs q1 then q0 else q1) 136 137 acceptable :: Interval -> QRational -> Bool 138 acceptable int o = cfBound o `encloses` int 139 140 141 142 143 144 145 encloses :: Interval -> Interval -> Bool 146 (i0,s0) `encloses` (i1,s1) 147 = if i0 > s0 then 148 if i1 < s1 then i0 <= i1 || s1 <= s0 149 else i0 <= i1 && s1 <= s0 150 else i1 < s1 && i0 <= i1 && s1 <= s0 151 152 153 154 intersect :: Interval -> Interval -> Interval 155 intersect (i,s) (i',s') = (rs,ri) where (ri,rs) = union [(s,i),(s',i')] 156 157 union :: [Interval] -> Interval 158 union = result . foldl join [] . sort . finitize 159 where result [int] = int 160 result [(i,_),(_,s)] = (i,s) 161 result ints = (0,0) 162 163 164 165 166 167 join :: [Interval] -> Interval -> [Interval] 168 join [] y = [y] 169 join xs@((i,s):iss) y@(i',s') = if s >= s' then xs else 170 if s >= i' then (i,s'):iss else y:xs 171 172 finitize 173 = concat . map split 174 where split x@(i,s) = if i < s then [x] else [(-1%%0,s),(i,1%%0)] 175 176 finiteInterval (i,s) = i < s 177 infiniteInterval (i,s) = i > s 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 algebraicAlgorithm :: Homography -> -- Homography 196 ContinuedFraction -> -- Input continued fraction 197 ContinuedFraction -- Output continued fraction 198 199 algebraicAlgorithm homography@([n0,n1],[d0,d1]) xs 200 = if n0*d1 /= n1*d0 then algebraicAlgorithmNoCheck homography xs 201 else if d1 /= 0 then rat2cf (n1%%d1) 202 else if d0 /= 0 then rat2cf (n0%%d0) 203 else error "algebraicAlgorithm: undefined homography" 204 205 algebraicAlgorithmNoCheck :: Homography -> -- Homography 206 ContinuedFraction -> -- Input CF 207 ContinuedFraction -- Output CF 208 209 algebraicAlgorithmNoCheck homography@([n0,n1],[d0,d1]) [] 210 = rat2cf (n0%%d0) 211 algebraicAlgorithmNoCheck homography@([n0,n1],[d0,d1]) xs'@(x:xs) 212 = case algebraicOutput homography (x%%1) of 213 Just o -> o : algebraicAlgorithmNoCheck ([d0,d1],[n0-d0*o,n1-d1*o]) xs' 214 Nothing -> algebraicAlgorithmNoCheck ([n0*x+n1,n0],[d0*x+d1,d0]) xs 215 216 algebraicOutput :: Homography -> -- Homography 217 QRational -> -- Leading Term of Input CF 218 Maybe Integer -- Leading Term of Output CF 219 220 algebraicOutput homography q 221 = if q == 0 then Nothing else a 222 where a = generateOutput algebraicBound algebraicEval homography [q] 223 224 225 226 227 228 229 230 algebraicBound :: Homography -> -- Homography 231 [QRational] -> -- Leading term of input CF 232 Interval -- Interval 233 234 algebraicBound homography@([n0,n1],[d0,d1]) [q] 235 = if n0*d1 == n1*d0 236 then (indp, indp) -- bound independent of x 237 else if xi == dcross || xs == dcross -- infinity at endpoint 238 then if xi > xs && (ncross <= xs || xi <= ncross) || 239 xi < xs && ncross <= xs && xi <= ncross 240 then zint -- zero in interval 241 else nzint -- zero not in interval 242 else if xi > xs && (d0 == 0 || xi < dcross || dcross < xs) || 243 xi < xs && d0 /= 0 && xi < dcross && dcross < xs 244 then (ys,yi) -- infinity in interval 245 else (yi,ys) -- infinity not in interval 246 where dcross = (-d1) %% d0 247 ncross = (-n1) %% n0 248 (xi,xs) = cfBound q 249 yi' = algebraicEval homography [xi] 250 ys' = algebraicEval homography [xs] 251 (yi,ys) = (min yi' ys', max yi' ys') 252 indp = if d0 /= 0 then n0%%d0 else 253 if d1 /= 0 then n1%%d1 else 254 error "algebraicBound: denominator 0" 255 zint = if qInfinite yi then if ys < 0%%1 then (ys, 1%%0) 256 else (-1%%0, ys) 257 else if yi < 0%%1 then (yi, 1%%0) 258 else (-1%%0, yi) 259 nzint = if qInfinite yi then if ys < 0%%1 then (-1%%0, ys) 260 else (ys, 1%%0) 261 else if yi < 0%%1 then (-1%%0, yi) 262 else (yi, 1%%0) 263 264 265 266 267 268 269 algebraicEval :: Homography -> -- Homography 270 [QRational] -> -- Input QRational 271 QRational -- Output QRational 272 273 algebraicEval ([n0,n1],[d0,d1]) [x] = (n0*xn + n1*xd) %% (d0*xn + d1*xd) 274 where xn = qNumerator x 275 xd = qDenominator x 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 quadraticAlgorithm :: Homography -> -- Homography 294 ContinuedFraction -> -- First Input CF 295 ContinuedFraction -> -- Second Input CF 296 ContinuedFraction -- Output CF 297 298 quadraticAlgorithm homography@([n0,n1,n2,n3],[d0,d1,d2,d3]) [] [] 299 = rat2cf (n0 %% d0) 300 quadraticAlgorithm homography@([n0,n1,n2,n3],[d0,d1,d2,d3]) [] ys'@(y:ys) 301 = if aaDet homography' == 0 302 then quadraticAlgorithm (qaInRight homography y) [] ys 303 else algebraicAlgorithmNoCheck homography' ys' 304 where homography' = ([n0,n2],[d0,d2]) 305 quadraticAlgorithm homography@([n0,n1,n2,n3],[d0,d1,d2,d3]) xs'@(x:xs) [] 306 = if aaDet homography' == 0 307 then quadraticAlgorithm (qaInLeft homography x) xs [] 308 else algebraicAlgorithmNoCheck homography' xs' 309 where homography' = ([n0,n1],[d0,d1]) 310 quadraticAlgorithm homography xs'@(x:xs) ys'@(y:ys) 311 = case quadraticOutput homography [x%%1,y%%1] of 312 Just o -> o : quadraticAlgorithm (qaOut homography o) xs' ys' 313 Nothing -> quadraticAlgorithm (qaIn homography x y) xs ys 314 315 quadraticOutput :: Homography -> -- Homography 316 [QRational] -> -- Leading Term of Input CF 317 Maybe Integer -- Leading Term of Output CF 318 quadraticOutput = generateOutput quadraticBound quadraticEval 319 320 qaIn :: Homography -> Integer -> Integer -> Homography 321 qaIn ([n0,n1,n2,n3],[d0,d1,d2,d3]) x y 322 = ([(n0*x+n1)*y+n2*x+n3, n0*y+n2, n0*x+n1, n0], 323 [(d0*x+d1)*y+d2*x+d3, d0*y+d2, d0*x+d1, d0]) 324 325 qaOut :: Homography -> Integer -> Homography 326 qaOut ([n0,n1,n2,n3],[d0,d1,d2,d3]) o 327 = ([d0,d1,d2,d3], [n0-d0*o, n1-d1*o, n2-d2*o, n3-d3*o]) 328 329 qaInLeft :: Homography -> Integer -> Homography 330 qaInLeft ([n0,n1,n2,n3],[d0,d1,d2,d3]) x 331 = ([n0*x+n1, n0, n2*x+n3, n2], [d0*x+d1, d0, d2*x+d3, d2]) 332 333 qaInRight :: Homography -> Integer -> Homography 334 qaInRight ([n0,n1,n2,n3],[d0,d1,d2,d3]) y 335 = ([n0*y+n2, n1*y+n3, n0, n1], [d0*y+d2, d1*y+d3, d0, d1]) 336 337 aaDet :: Homography -> Integer 338 aaDet ([n0,n1],[d0,d1]) = n0*d1 - n1*d0 339 340 341 342 quadraticEval :: Homography -> -- homography 343 [QRational] -> -- QRational 344 QRational -- result 345 346 quadraticEval (ns,ds) [x, y] 347 = ((n0*x+n1)*y + (n2*x+n3)) / ((d0*x+d1)*y + (d2*x+d3)) 348 where [n0,n1,n2,n3] = map (%% 1) ns 349 [d0,d1,d2,d3] = map (%% 1) ds 350 351 quadraticBound :: Homography -> -- homography 352 [QRational] -> -- terms of CFs 353 Interval -- interval 354 355 quadraticBound ([n0,n1,n2,n3],[d0,d1,d2,d3]) [x, y] 356 = union ints 357 where (xi,xs) = cfBound x 358 (yi,ys) = cfBound y 359 ints = [algebraicBound (aax x) [y] | x <- [xi,xs]] ++ 360 [algebraicBound (aay y) [x] | y <- [yi,ys]] 361 aax x = ([n0*xn+n1*xd,n2*xn+n3*xd],[d0*xn+d1*xd,d2*xn+d3*xd]) 362 where xn = qNumerator x 363 xd = qDenominator x 364 aay y = ([n0*yn+n2*yd,n1*yn+n3*yd],[d0*yn+d2*yd,d1*yn+d3*yd]) 365 where yn = qNumerator y 366 yd = qDenominator y 367 368 369 370 371 372 373 374 375 cfRat2CFSqrt :: QRational -> ContinuedFraction 376 cfRat2CFSqrt q 377 = map snd (next ([2*d*z,d,n-d*z*z,0], z)) 378 where n = qNumerator q 379 d = qDenominator q 380 s = intSquareRoot (n*d) 381 z = qRound (s%%d) 382 t = s - z*d 383 next :: ([Integer],Integer) -> [([Integer],Integer)] 384 next (s,z) = if d0 == 0 385 then [(s,z)] 386 else (s,z): next ([n0',n1',d0',d1'],z') 387 where [n0,n1,d0,d1] = s 388 n0' = d0*z'+d1 389 n1' = d0 390 d0' = n0*z'+n1-z'*n0' 391 d1' = n0-z'*d0 392 z' = qRound ((n0+t) %% d0) 393 394 395 396 397 398 399 intSquareRoot :: Integer -> Integer 400 intSquareRoot x = until satisfy (improve x) x 401 where satisfy y = y2+y >= x && y2-y < x 402 where y2 = y*y 403 improve x y = (y*y+x) `ddiv` (2*y) 404 ddiv x y = if (r <= y `div` 2) then q else q+1 405 where (q,r) = divMod x y 406 407 root2 = cfRat2CFSqrt (2%%1) 408 root3 = cfRat2CFSqrt (3%%1) 409 test0 :: ContinuedFraction 410 test0 = f 1 where f n = 0: 10^n: f (n+1) 411 412 413 414 415 416 417 418 419 decimals = 10 :: Int 420 421 accuracy = 1 %% (10^(toInteger decimals)) :: QRational 422 423 424 425 426 cf2Rat :: ContinuedFraction -> QRational 427 cf2Rat = outputAlgorithm ([1,0],[0,1]) 428 429 outputAlgorithm :: Homography -> -- Homography 430 ContinuedFraction -> -- Input CF 431 QRational -- Final Approximation 432 433 outputAlgorithm homography@([n0,n1],[d0,d1]) [] = n0 %% d0 434 outputAlgorithm homography@([n0,n1],[d0,d1]) xs'@(x:xs) 435 = if i <= s && s-i < accuracy 436 then algebraicEval homography [x%%1] 437 else outputAlgorithm ([n0*x+n1,n0],[d0*x+d1,d0]) xs 438 where (i,s) = algebraicBound homography [x%%1] 439 440 441 integerFraction :: ContinuedFraction -> -- Input CF 442 (Integer,ContinuedFraction) -- Answer 443 integerFraction x = (o, algebraicAlgorithm h' f) 444 where (h,f) = intFrac (([1,0],[0,1]),x) 445 o = qRound (algebraicEval h [head f%%1]) 446 h' = ([n0-d0*o,n1-d1*o],[d0,d1]) 447 where ([n0,n1],[d0,d1]) = h 448 449 intFrac :: (Homography,ContinuedFraction) -> (Homography,ContinuedFraction) 450 intFrac st@(_,[_]) = st 451 intFrac st@(homography@([n0,n1],[d0,d1]),(x:xs)) 452 = if i <= s && s-i < 1 %% 10 453 then st 454 else intFrac (([n0*x+n1,n0],[d0*x+d1,d0]), xs) 455 where (i,s) = algebraicBound homography [x%%1] 456 457 458