1 2 3 4 5 6 7 8 9 10 11 12 13 module Transcendentals 14 (module ContinuedFractions, 15 exp1, expQ, expR, log10, logQ, logR, tanQ, tanR, cotR, atan1, atanQ, atanR) 16 where 17 import ContinuedFractions 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 exp1CF :: [QRational] 63 exp1CF = 1: f 0 where f n = (4*n+1): (-2): (-(4*n+3)): 2: f (n+1) 64 65 66 67 68 69 70 71 72 73 74 expQ :: QRational -> ContinuedFraction 75 expQ x = aaQ hom ts2 76 where (ts1,ts2) = splitAt 2 cf 77 hom = absorbQ startH ts1 78 (startH,cf) = if abs x < 1 then 79 (([1,0],[0,1]), oddQ x exp1CF) else 80 (([qDenominator x,0],[0,qNumerator x]), 81 evenQ (1/x) exp1CF) 82 83 exp1 :: ContinuedFraction 84 exp1 = expQ (1%%1) 85 86 87 88 89 expR :: ContinuedFraction -> ContinuedFraction 90 expR x = if mightBeZero x 91 then processQS x (splitAt 6 (oddR exp1CF)) 92 else quadraticAlgorithm ([0,1,1,0],[0,1,0,0]) x cf2 93 where cf2 = processQS x (splitAt 6 (oddR (tail exp1CF))) 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 logQ :: QRational -> ContinuedFraction 111 logQ x = aaQ startH ts2 112 where (ts1,ts2) = splitAt 2 (oddQ (-t*t) atan1CF) 113 t = (x-1)/(x+1) 114 tn = qNumerator t 115 td = qDenominator t 116 startH = absorbQ ([-2*td,0],[0,tn]) ts1 117 118 log10 :: ContinuedFraction 119 log10 = logQ (10%%1) 120 121 122 123 124 125 logR :: ContinuedFraction -> ContinuedFraction 126 logR x = quadraticAlgorithm ([2,0,0,0],[0,0,0,1]) t cf1 127 where t = algebraicAlgorithm ([1,(-1)],[1,1]) x 128 t2 = quadraticAlgorithm ([-1,0,0,0],[0,0,0,1]) t t 129 cf = qsSetupAlt (drop 1 (oddR atan1CF)) t2 130 cf1 = quadraticAlgorithm ([0,1,0,0],[0,1,1,0]) t2 cf 131 132 133 134 135 136 137 138 139 140 141 142 143 tan1CF :: [QRational] 144 tan1CF = 0: f 0 where f n = (4*n+1) : (-(4*n+3)) : f (n+1) 145 146 tanQ :: QRational -> ContinuedFraction 147 tanQ x = if x /= 0 then aaQ startH ts2 else error "tanQ: 0 argument" 148 where n = qRound (abs (x/2)) + 1 149 (ts1,ts2) = splitAt (fromIntegral n) (bothQ x tan1CF) 150 startH = absorbQ ([1,0],[0,1]) ts1 151 152 tanR :: ContinuedFraction -> ContinuedFraction 153 tanR x = quadraticAlgorithm ([0,1,0,0],[0,0,1,0]) cfd cfn 154 where cf = qsSetupAlt (drop 1 (oddR tan1CF)) x2 155 x2 = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x x 156 cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x cf 157 cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf 158 159 cotR :: ContinuedFraction -> ContinuedFraction 160 cotR x = quadraticAlgorithm ([0,1,0,0],[0,0,1,0]) cfn cfd 161 where cf = qsSetupAlt (drop 1 (oddR tan1CF)) x2 162 x2 = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x x 163 cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x cf 164 cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf 165 166 167 168 169 170 171 172 atan1CF :: [QRational] 173 atan1CF = 0: 1: 3: 5/4: f (8/9) 0 174 where f :: QRational -> QRational -> [QRational] 175 f p n = (7/2 + 2*n)*p: (9/2 + 2*n)*p': f p'' (n+1) 176 where p' = 1/(p*(n+2)*(n+2)) 177 p'' = 1/(p'*(n+5/2)*(n+5/2)) 178 179 180 181 182 183 184 185 186 187 atanQ :: QRational -> ContinuedFraction 188 atanQ x 189 = if x > 1 then 190 quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1]) atan1 (cf (1/x)) else 191 if x < -1 then 192 quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1]) atan1 (cf (1/x)) else 193 cf x 194 where 195 cf x = aaQ startH ts2 196 where 197 (ts1,ts2) = splitAt 4 (oddQ (x*x) atan1CF) 198 startH = absorbQ ([qDenominator x,0],[0,qNumerator x]) ts1 199 200 201 202 203 204 atan1 :: ContinuedFraction 205 atan1 = atanQ (1%%1) 206 207 atanR :: ContinuedFraction -> ContinuedFraction 208 atanR [] = algebraicAlgorithm ([2,0],[0,1]) atan1 209 atanR xs@(x:xs') 210 = if x > 1 then 211 quadraticAlgorithm ([0,-1,2,0] ,[0,0,0,1]) atan1 ys else 212 if x < -1 then 213 quadraticAlgorithm ([0,-1,-2,0],[0,0,0,1]) atan1 ys else 214 atanRdirect xs 215 where ys = atanRdirect (algebraicAlgorithm ([0,1],[1,0]) xs) 216 217 atanRdirect :: ContinuedFraction -> ContinuedFraction 218 atanRdirect x 219 = quadraticAlgorithm ([0,0,1,0],[0,1,0,0]) cfn cfd 220 where x2 = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x x 221 cf = qsSetupAlt (drop 1 (oddR atan1CF)) x2 222 cfn = quadraticAlgorithm ([1,0,0,0],[0,0,0,1]) x cf 223 cfd = quadraticAlgorithm ([0,1,1,0],[0,0,0,1]) x2 cf 224 225 226 227 228 229 mightBeZero :: ContinuedFraction -> Bool 230 mightBeZero [] = False 231 mightBeZero (x:xs) = x == 0 232 233 234 235 236 237 238 239 240 241 aaQ :: Homography -> [QRational] -> ContinuedFraction 242 aaQ homography@([n0,n1],[d0,d1]) xs 243 = if n0*d1 /= n1*d0 then aaQNoCheck homography xs 244 else if d1 /= 0 then rat2cf (n1%%d1) 245 else if d0 /= 0 then rat2cf (n0%%d0) 246 else error "aaQ: undefined homography" 247 248 249 250 251 aaQNoCheck homography@([n0,n1],[d0,d1]) [] 252 = rat2cf (n0%%d0) 253 aaQNoCheck homography@([n0,n1],[d0,d1]) xs'@(x:xs) 254 = case algebraicOutput homography x of 255 Just o -> o : aaQNoCheck ([d0,d1], [n0-d0*o,n1-d1*o]) xs' 256 Nothing -> aaQNoCheck ([n0*n+n1*d,n0*d], [d0*n+d1*d,d0*d]) xs 257 where (n,d) = (qNumerator x, qDenominator x) 258 259 260 261 262 absorbQ :: Homography -> [QRational] -> Homography 263 absorbQ homography [] = homography 264 absorbQ ([n0,n1],[d0,d1]) (q:qs) = absorbQ ([n0*n+n1*d,n0*d], 265 [d0*n+d1*d,d0*d]) qs 266 where n = qNumerator q 267 d = qDenominator q 268 269 270 271 272 bothQ :: QRational -> [QRational] -> [QRational] 273 bothQ x (x1:xs) = x1/x: bothQ x xs 274 275 evenQ :: QRational -> [QRational] -> [QRational] 276 evenQ x (x1:x2:xs) = x1/x: x2: evenQ x xs 277 278 oddQ :: QRational -> [QRational] -> [QRational] 279 oddQ x (x1:x2:xs) = x1: x2/x: oddQ x xs 280 281 282 283 284 285 286 287 288 type State = (QRational, Homography) 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 oddR :: [QRational] -> [State] 304 oddR (x1:x2:xs) = (x1,([0,xn,xd,0],[xd,0,0,0])): oddR xs 305 where (xn,xd) = (qNumerator x2, qDenominator x2) 306 307 308 309 310 processQS :: ContinuedFraction -> ([State],[State]) -> ContinuedFraction 311 processQS x (hs1,hs2) = preProcess x hs1 (qsSetupAlt hs2 x) 312 313 314 315 316 317 preProcess :: ContinuedFraction -> -- the argument CF 318 [State] -> -- leading parts of infinite state 319 ContinuedFraction -> -- accumulating CF 320 ContinuedFraction -- the resulting CF 321 preProcess _ [] c = c 322 preProcess x ((i,h):hs) c = algebraicAlgorithm ([ni,nd],[nd,0]) 323 (quadraticAlgorithm h x (preProcess x hs c)) 324 where ni = qNumerator i 325 nd = qDenominator i 326 327 328 329 qsSetupAlt :: [State] -> ContinuedFraction -> ContinuedFraction 330 qsSetupAlt hs x = aaQ ([1,0],[0,1]) (qsProcess (hs, x)) 331 332 333 334 335 336 337 338 339 340 341 342 343 qsProcess :: ([State], ContinuedFraction) -> [QRational] 344 qsProcess (ihs, []) = asProcess ihs 345 qsProcess st@(ihs, (_:_)) = qsNext 1 st 346 347 348 349 350 351 352 353 354 355 qsNext :: Int -> ([State], ContinuedFraction) -> [QRational] 356 qsNext n (ihs, []) = asProcess (qsEnd ihs) 357 qsNext n st@(((i,h):ihs), xs'@(x:xs)) 358 = if length outs > 1 359 then init outs ++ qsNext (n+1) (ihs',xs') 360 else qsNext (n+1) (qsIn x ((i,h): ihs), xs) 361 where (outs, (ihs',_)) = qsTerm n st 362 363 364 365 366 367 368 qsTerm :: Int -> ([State], ContinuedFraction) 369 -> ([QRational],([State], ContinuedFraction)) 370 qsTerm n (ihs,xs) 371 = (outs, (ihs',xs)) 372 where (front,back@((i,_):_)) = splitAt n ihs 373 (outs,ihs') = foldr (qsTerm' (head xs)) ([i],back) front 374 375 376 377 qsTerm' :: Integer -> State -> ([QRational], [State]) -> 378 ([QRational], [State]) 379 qsTerm' x ih@(i,h) (is,ihs) 380 = (is',((last is',h''):ihs)) 381 where h' = qsAbsorb (init is) h 382 i' = last is 383 (is',h'') = qsEmit x i' ([i],h') 384 385 386 387 388 qsAbsorb :: [QRational] -> Homography -> Homography 389 qsAbsorb [] h = h 390 qsAbsorb (i:is) h@([n0,n1,n2,n3],[d0,d1,d2,d3]) 391 = qsAbsorb is ([n0*n+n2*d, n1*n+n3*d, n0*d, n1*d], 392 [d0*n+d2*d, d1*n+d3*d, d0*d, d1*d]) 393 where n = qNumerator i 394 d = qDenominator i 395 396 397 398 399 qsEmit :: Integer -> QRational -> ([QRational], Homography) -> 400 ([QRational], Homography) 401 qsEmit x i ish@(is,h@([n0,n1,n2,n3],[d0,d1,d2,d3])) 402 = case quadraticOutput h [x%%1,i] of 403 Just o -> qsEmit x i (is++[o%%1], 404 ([d0,d1,d2,d3], 405 [n0-d0*o, n1-d1*o, n2-d2*o, n3-d3*o])) 406 Nothing -> ish 407 408 409 410 411 qsIn :: Integer -> [State] -> [State] 412 qsIn = map . input 413 where input x (i,([n0,n1,n2,n3],[d0,d1,d2,d3])) 414 = (i, ([n0*x+n1, n0, n2*x+n3, n2], 415 [d0*x+d1, d0, d2*x+d3, d2])) 416 417 418 419 420 421 422 423 424 425 426 qsEnd :: [State] -> [State] 427 qsEnd = map input 428 where input (i,([n0,_,n2,_],[d0,_,d2,_])) 429 = (i,([n0,n2],[d0,d2])) 430 431 432 433 asProcess :: [State] -> [QRational] 434 asProcess ihs = fst (head ihs) : asProcess (asNext ihs) 435 436 asNext :: [State] -> [State] 437 asNext ((i,h@([n0,n1],[d0,d1])):ihs'@((i',_):ihs)) 438 = case algebraicOutput h i' of 439 Just o -> (o%%1, ([d0,d1], [n0-d0*o,n1-d1*o])): ihs' 440 Nothing -> asNext ((i,hIn i' h): asNext ihs') 441 where hIn x ([n0,n1],[d0,d1]) = ([n0*xn+n1*xd,n0*xd],[d0*xn+d1*xd,d0*xd]) 442 where xn = qNumerator x 443 xd = qDenominator x 444 445 446 447 448