1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 module RealReals (RealReal) where 17 import Transcendentals 18 import Ratio 19 import List( genericLength ) 20 import Numeric( readSigned ) 21 22 data RealReal = RealInt Integer | 23 RealRat QRational | 24 RealCF ContinuedFraction 25 26 27 28 29 30 31 instance Eq RealReal where 32 (==) = realRealComp (==) 33 (/=) = realRealComp (/=) 34 35 36 37 38 realRealComp :: (QRational -> QRational -> Bool) -> -- Operator 39 RealReal -> -- First Argument 40 RealReal -> -- Second Argument 41 Bool -- Result 42 realRealComp op x y = op (makeRational (x-y)) 0 43 44 45 46 47 48 instance Ord RealReal where 49 (<=) = realRealComp (<=) 50 (<) = realRealComp (<) 51 (>=) = realRealComp (>=) 52 (>) = realRealComp (>) 53 54 55 56 57 58 instance Enum RealReal where 59 enumFrom = iterate (+1) 60 enumFromThen n m = iterate (+(m-n)) n 61 62 63 64 65 66 instance Num RealReal where 67 x + y = realRealAdd x y 68 negate x = realRealNegate x 69 x * y = realRealMultiply x y 70 abs x = realRealAbs x 71 signum x = realRealSignum x 72 fromInteger x = RealInt x 73 74 75 76 realRealAdd :: RealReal -> RealReal -> RealReal 77 realRealAdd (RealInt x) (RealInt y) = RealInt (x+y) 78 realRealAdd (RealInt x) (RealRat y) = wrapRat ((x%%1)+y) 79 realRealAdd (RealInt x) (RealCF y) = RealCF (algebraicAlgorithm h y) 80 where h = ([1,x],[0,1]) 81 realRealAdd (RealRat x) (RealInt y) = wrapRat (x+(y%%1)) 82 realRealAdd (RealRat x) (RealRat y) = wrapRat (x+y) 83 realRealAdd (RealRat x) (RealCF y) = RealCF (algebraicAlgorithm h y) 84 where nx = qNumerator x 85 dx = qDenominator x 86 h = ([dx,nx],[0,dx]) 87 realRealAdd (RealCF x) (RealInt y) = RealCF (algebraicAlgorithm h x) 88 where h = ([1,y],[0,1]) 89 realRealAdd (RealCF x) (RealRat y) = RealCF (algebraicAlgorithm h x) 90 where ny = qNumerator y 91 dy = qDenominator y 92 h = ([dy,ny],[0,dy]) 93 realRealAdd (RealCF x) (RealCF y) = RealCF (quadraticAlgorithm h x y) 94 where h = ([0,1,1,0],[0,0,0,1]) 95 96 97 98 99 wrapRat :: QRational -> RealReal 100 wrapRat x = if qDenominator x == 1 then RealInt (qNumerator x) else 101 RealRat x 102 103 104 105 106 107 108 109 realRealNegate :: RealReal -> RealReal 110 realRealNegate (RealInt x) = RealInt (negate x) 111 realRealNegate (RealRat x) = RealRat (negate x) 112 realRealNegate (RealCF x) = RealCF (map negate x) 113 114 115 116 realRealMultiply :: RealReal -> RealReal -> RealReal 117 realRealMultiply (RealInt x) (RealInt y) = RealInt (x*y) 118 realRealMultiply (RealInt x) (RealRat y) = wrapRat ((x%%1)*y) 119 realRealMultiply (RealInt x) (RealCF y) = RealCF (algebraicAlgorithm h y) 120 where h = ([x,0],[0,1]) 121 realRealMultiply (RealRat x) (RealInt y) = wrapRat (x*(y%%1)) 122 realRealMultiply (RealRat x) (RealRat y) = wrapRat (x*y) 123 realRealMultiply (RealRat x) (RealCF y) = RealCF (algebraicAlgorithm h y) 124 where nx = qNumerator x 125 dx = qDenominator x 126 h = ([nx,0],[0,dx]) 127 realRealMultiply (RealCF x) (RealInt y) = RealCF (algebraicAlgorithm h x) 128 where h = ([y,0],[0,1]) 129 realRealMultiply (RealCF x) (RealRat y) = RealCF (algebraicAlgorithm h x) 130 where ny = qNumerator y 131 dy = qDenominator y 132 h = ([ny,0],[0,dy]) 133 realRealMultiply (RealCF x) (RealCF y) = RealCF (quadraticAlgorithm h x y) 134 where h = ([1,0,0,0],[0,0,0,1]) 135 136 137 138 139 140 141 142 realRealAbs :: RealReal -> RealReal 143 realRealAbs x = if realRealSignum x < 0 then realRealNegate x else x 144 145 146 147 148 149 realRealSignum :: RealReal -> RealReal 150 realRealSignum x = wrapRat (signum (makeRational x)) 151 152 153 154 155 156 instance Real RealReal where 157 toRational rx = qNumerator x % qDenominator x 158 where x = makeRational rx 159 160 161 162 163 164 instance Fractional RealReal where 165 x / y = realRealDivide x y 166 fromRational x = wrapRat (numerator x %% denominator x) 167 168 realRealDivide :: RealReal -> RealReal -> RealReal 169 realRealDivide (RealInt x) (RealInt y) = wrapRat (x%%y) 170 realRealDivide (RealInt x) (RealRat y) = wrapRat ((x%%1)/y) 171 realRealDivide (RealInt x) (RealCF y) = RealCF (algebraicAlgorithm h y) 172 where h = ([0,x],[1,0]) 173 realRealDivide (RealRat x) (RealInt y) = wrapRat (x/(y%%1)) 174 realRealDivide (RealRat x) (RealRat y) = wrapRat (x/y) 175 realRealDivide (RealRat x) (RealCF y) = RealCF (algebraicAlgorithm h y) 176 where nx = qNumerator x 177 dx = qDenominator x 178 h = ([0,nx],[dx,0]) 179 realRealDivide (RealCF x) (RealInt y) = RealCF (algebraicAlgorithm h x) 180 where h = ([1,0],[0,y]) 181 realRealDivide (RealCF x) (RealRat y) = RealCF (algebraicAlgorithm h x) 182 where ny = qNumerator y 183 dy = qDenominator y 184 h = ([dy,0],[0,ny]) 185 realRealDivide (RealCF x) (RealCF y) = RealCF (quadraticAlgorithm h x y) 186 where h = ([0,0,1,0],[0,1,0,0]) 187 188 189 190 instance Floating RealReal where 191 pi = realRealPi 192 exp = realRealExp 193 log = realRealLog 194 sqrt = realRealSqrt 195 sin = realRealSin 196 cos = realRealCos 197 tan = realRealTan 198 asin = realRealAsin 199 acos = realRealAcos 200 atan = realRealAtan 201 sinh = realRealSinh 202 cosh = realRealCosh 203 tanh = realRealTanh 204 asinh = realRealAsinh 205 acosh = realRealAcosh 206 atanh = realRealAtanh 207 208 209 210 realRealPi :: RealReal 211 realRealPi = 4 * (RealCF atan1) 212 213 214 215 realRealExp :: RealReal -> RealReal 216 realRealExp (RealInt x) = RealCF exp1 ^^ x 217 realRealExp (RealRat x) = RealCF exp1 ^^ i * RealCF (expQ (x-(i%%1))) 218 where i = qRound x 219 realRealExp (RealCF x) = RealCF exp1 ^^ i * RealCF (expR x') 220 where (i,x') = integerFraction x 221 222 realRealLog :: RealReal -> RealReal 223 realRealLog (RealInt x) 224 = if x==1 then RealInt 0 else RealCF log10 * RealInt i + RealCF (logQ x') 225 where i = genericLength (show x):: Integer 226 x' = x %% (10^i) 227 realRealLog (RealRat x) 228 = RealCF log10 * RealInt i + RealCF (logQ x') 229 where i = if x > 1 then f x else negate (f (1/x)) 230 f :: QRational -> Integer 231 f x = genericLength (show (qRound x + 1)) 232 x' = x / ((10%%1)^^i) 233 realRealLog (RealCF x) 234 = RealCF log10 * RealInt i + RealCF (logR x') 235 where i = if head x /= 0 then f x else negate (f (tail x)) 236 f :: ContinuedFraction -> Integer 237 f x = genericLength (show (i + 1)) where (i,_) = integerFraction x 238 m = (1%%10)^^i 239 x' = algebraicAlgorithm ([qNumerator m,0],[0,qDenominator m]) x 240 241 realRealSqrt :: RealReal -> RealReal 242 realRealSqrt (RealInt x) = RealCF (cfRat2CFSqrt (x%%1)) 243 realRealSqrt (RealRat x) = RealCF (cfRat2CFSqrt x) 244 realRealSqrt rx@(RealCF _) = exp (log rx / 2) 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 realRealSin :: RealReal -> RealReal 260 realRealSin rx@(RealCF _) 261 = 2*c/(1+c*c) 262 where (RealCF x2) = rx/2 263 cotx2 = cotR x2 264 tanx2 = tanR x2 265 cf = if head cotx2 == 0 then cotx2 else tanx2 266 c = RealCF cf 267 realRealSin x = 2*t / (1+t*t) 268 where t = tan (x/2) 269 270 realRealCos :: RealReal -> RealReal 271 realRealCos rx@(RealCF _) 272 = if head cotx2 == 0 then (c2-1)/(c2+1) else (1-c2)/(1+c2) 273 where (RealCF x2) = rx/2 274 cotx2 = cotR x2 275 tanx2 = tanR x2 276 cf = if head cotx2 == 0 then cotx2 else tanx2 277 c = RealCF cf 278 c2 = c*c 279 realRealCos x = (1-t2) / (1+t2) 280 where t = tan (x/2) 281 t2 = t*t 282 283 realRealTan :: RealReal -> RealReal 284 realRealTan (RealInt x) = if x==0 then RealInt 0 else RealCF (tanQ (x%%1)) 285 realRealTan (RealRat x) = RealCF (tanQ x) 286 realRealTan rx@(RealCF _) 287 = RealCF (tanR x) 288 where (RealCF x') = rx / pi 289 (i,_) = integerFraction x' 290 (RealCF x) = rx - (RealInt i * pi) 291 292 293 294 realRealAsin :: RealReal -> RealReal 295 realRealAsin x = atan (x / sqrt (1-x*x)) 296 297 realRealAcos :: RealReal -> RealReal 298 realRealAcos x = atan (sqrt (1-x*x) / x) 299 300 realRealAtan :: RealReal -> RealReal 301 realRealAtan (RealInt x) = if x==0 then RealInt 0 else RealCF (atanQ (x%%1)) 302 realRealAtan (RealRat x) = RealCF (atanQ x) 303 realRealAtan (RealCF x) = RealCF (atanR x) 304 305 306 307 realRealSinh :: RealReal -> RealReal 308 realRealSinh x = (ex - (1 / ex)) / 2 309 where ex = exp x 310 311 realRealCosh :: RealReal -> RealReal 312 realRealCosh x = (ex + (1 / ex)) / 2 313 where ex = exp x 314 315 realRealTanh :: RealReal -> RealReal 316 realRealTanh x = (ex - ex') / (ex + ex') 317 where ex = exp x 318 ex' = 1/ex 319 320 321 322 realRealAsinh :: RealReal -> RealReal 323 realRealAsinh x = log (x + sqrt (x*x+1)) 324 325 realRealAcosh :: RealReal -> RealReal 326 realRealAcosh x = log (x + sqrt (x*x-1)) 327 328 329 330 realRealAtanh :: RealReal -> RealReal 331 realRealAtanh x = log ((1+x)/(1-x)) / 2 332 333 334 335 instance Read RealReal where 336 readsPrec p = readSigned readRealReal 337 instance Show RealReal where 338 showsPrec p rx = showString (showRat p (makeRational rx)) 339 340 makeRational :: RealReal -> QRational 341 makeRational (RealInt x) = (x%%1) 342 makeRational (RealRat x) = x 343 makeRational (RealCF x) = cf2Rat x 344 345 346 347 348 349 350 readRealReal :: ReadS RealReal 351 readRealReal r = [(RealInt x, s) | (x,s) <- reads r] 352 353 354 355 356 357 358 359 showRat :: Int -> QRational -> String 360 showRat p x = if s == "-" && p > 6 then "(" ++ res ++ ")" else res 361 where z = qRound (x/accuracy) 362 s = if z < 0 then "-" else "" 363 zs = show (abs z) 364 l = length zs 365 zs' = pad ++ zs 366 pad = if (l <= decimals) 367 then take (decimals-l+1) zeros else "" 368 zeros = repeat '0' 369 (i,f) = splitAt (length zs' - decimals) zs' 370 res = s ++ i ++ "." ++ f 371 372 373