module Data.Number.CReal(CReal, showCReal) where
import Data.Ratio
import Numeric(readFloat, readSigned)
data CReal = CR (Int -> Integer)
instance Eq CReal where
x == y = s' (digitsToBits digits) == 0 where (CR s') = xy
instance Ord CReal where
x <= y = s' (digitsToBits digits) <= 0 where (CR s') = xy
x < y = s' (digitsToBits digits) < 0 where (CR s') = xy
x >= y = s' (digitsToBits digits) >= 0 where (CR s') = xy
x > y = s' (digitsToBits digits) > 0 where (CR s') = xy
max (CR x') (CR y') = CR (\p -> max (x' p) (y' p))
min (CR x') (CR y') = CR (\p -> min (x' p) (y' p))
instance Num CReal where
(CR x') + (CR y') = CR (\p -> round_uk ((x' (p+2) + y' (p+2)) % 4))
(CR x') * (CR y') = CR (\p -> round_uk ((x' (p+sy)*y' (p+sx)) % 2^(p+sx+sy)))
where x0 = abs (x' 0)+2; y0 = abs (y' 0)+2
sx = sizeinbase x0 2+3; sy = sizeinbase y0 2+3
negate (CR x') = CR (\p -> negate (x' p))
abs x = max x (negate x)
signum (CR x') = fromInteger (signum (x' (digitsToBits digits)))
fromInteger n = CR (\p -> n*2^p)
instance Fractional CReal where
recip (CR x') = CR (\p -> let s = head [n | n <- [0..], 3 <= abs (x' n)]
in round_uk (2^(2*p+2*s+2) % (x' (p+2*s+2))))
fromRational x = fromInteger (numerator x) / fromInteger (denominator x)
div2n :: CReal -> Int -> CReal
div2n (CR x') n = CR (\p -> if p >= n then x' (pn) else round_uk (x' p % 2^n))
mul2n :: CReal -> Int -> CReal
mul2n (CR x') n = CR (\p -> x' (p+n))
instance Floating CReal where
pi = 16 * atan (fromRational (1 % 5))
4 * atan (fromRational (1 % 239))
sqrt x = CR (\p -> floorsqrt (x' (2*p))) where (CR x') = x
log x = if t < 0 then error "log of negative number\n" else
if t < 4 then log (recip x) else
if t < 8 then log_dr x else
log_dr (div2n x n) + fromIntegral n * log2
where (CR x') = x; t = x' 2; n = sizeinbase t 2 3
exp x = if n < 0 then div2n (exp_dr s) (fromInteger (n)) else
if n > 0 then mul2n (exp_dr s) (fromInteger n) else exp_dr s
where (CR u') = x/log2; n = u' 0; s = xfromInteger n*log2
sin x = if n == 0 then sin_dr y else
if n == 1 then sqrt1By2 * (cos_dr y + sin_dr y) else
if n == 2 then cos_dr y else
if n == 3 then sqrt1By2 * (cos_dr y sin_dr y) else
if n == 4 then sin_dr y else
if n == 5 then sqrt1By2 * (cos_dr y + sin_dr y) else
if n == 6 then cos_dr y else
sqrt1By2 * (cos_dr y sin_dr y)
where (CR z') = x/piBy4; s = round_uk (z' 2 % 4); n = s `mod` 8
y = x piBy4 * fromInteger s
cos x = if n == 0 then cos_dr y else
if n == 1 then sqrt1By2 * (cos_dr y sin_dr y) else
if n == 2 then sin_dr y else
if n == 3 then sqrt1By2 * (cos_dr y + sin_dr y) else
if n == 4 then cos_dr y else
if n == 5 then sqrt1By2 * (cos_dr y sin_dr y) else
if n == 6 then sin_dr y else
sqrt1By2 * (cos_dr y + sin_dr y)
where (CR z') = x/piBy4; s = round_uk (z' 2 % 4); n = s `mod` 8
y = x piBy4 * fromInteger s
atan x = if t < 5 then atan_dr (negate (recip x)) piBy2 else
if t == 4 then piBy4 atan_dr (xp1/xm1) else
if t < 4 then atan_dr x else
if t == 4 then piBy4 + atan_dr (xm1/xp1) else
piBy2 atan_dr (recip x)
where (CR x') = x; t = x' 2
xp1 = x+1; xm1 = x1
asin x = if x0 > 0 then pi / 2 atan (s/x) else
if x0 == 0 then atan (x/s) else
atan (s/x) pi / 2
where (CR x') = x; x0 = x' 0; s = sqrt (1 x*x)
acos x = pi / 2 asin x
sinh x = (y recip y) / 2 where y = exp x
cosh x = (y + recip y) / 2 where y = exp x
tanh x = (y y') / (y + y') where y = exp x; y' = recip y
asinh x = log (x + sqrt (x*x + 1))
acosh x = log (x + sqrt (x*x 1))
atanh x = log ((1 + x) / (1 x)) / 2
acc_seq :: (Rational -> Integer -> Rational) -> [Rational]
acc_seq f = scanl f (1 % 1) [1..]
exp_dr :: CReal -> CReal
exp_dr = power_series (acc_seq (\a n -> a*(1 % n))) id
log_dr :: CReal -> CReal
log_dr x = y * log_drx y where y = (x 1) / x
log_drx :: CReal -> CReal
log_drx = power_series [1 % n | n <- [1..]] (+1)
sin_dr :: CReal -> CReal
sin_dr x = x*power_series (acc_seq (\a n -> a*(1 % (2*n*(2*n+1))))) id (x*x)
cos_dr :: CReal -> CReal
cos_dr x = power_series (acc_seq (\a n -> a*(1 % (2*n*(2*n1))))) id (x*x)
atan_dr :: CReal -> CReal
atan_dr x = (x/y) * atan_drx ((x*x)/y) where y = x*x+1
atan_drx :: CReal -> CReal
atan_drx = power_series (acc_seq (\a n -> a*((2*n) % (2*n+1)))) (+1)
power_series :: [Rational] -> (Int -> Int) -> CReal -> CReal
power_series ps terms (CR x')
= CR (\p -> let t = terms p; l2t = 2*sizeinbase (toInteger t+1) 2+6; p' = p + l2t
xr = x' p'; xn = 2^p'; g yn = round_uk ((yn*xr) % (2^p'))
in round_uk (accumulate (iterate g xn) (take t ps) % (2^l2t)))
where accumulate _ [] = 0
accumulate [] _ = error "CReal.power_series.accumulate"
accumulate (x:xs) (c:cs) = let t = round_uk (c*(x % 1)) in
if t == 0 then 0 else t + accumulate xs cs
piBy2 :: CReal
piBy2 = div2n pi 1
piBy4 :: CReal
piBy4 = div2n pi 2
log2 :: CReal
log2 = div2n (log_drx (recip 2)) 1
sqrt1By2 :: CReal
sqrt1By2 = sqrt (recip 2)
instance Enum CReal where
toEnum i = fromIntegral i
fromEnum _ = error "Cannot fromEnum CReal"
enumFrom = iterate (+ 1)
enumFromTo n e = takeWhile (<= e) $ iterate (+ 1)n
enumFromThen n m = iterate (+(mn)) n
enumFromThenTo n m e = if m >= n then takeWhile (<= e) $ iterate (+(mn)) n
else takeWhile (>= e) $ iterate (+(mn)) n
instance Real CReal where
toRational _ = error "CReal.toRational"
instance RealFrac CReal where
properFraction x@(CR x') = (fromInteger n, x fromInteger n) where n = x' 0
instance RealFloat CReal where
floatRadix _ = error "CCeal.floatRadix"
floatDigits _ = error "CReal.floatDigits"
floatRange _ = error "CReal.floatRange"
decodeFloat _ = error "CReal.decodeFloat"
encodeFloat _ _ = error "CReal.encodeFloat"
exponent _ = 0
scaleFloat 0 x = x
significand x = x
isNaN _ = False
isInfinite _ = False
isDenormalized _ = False
isNegativeZero _ = False
isIEEE _ = False
showCReal :: Int
-> CReal
-> String
showCReal d (CR x')
= (if s then "-" else "") ++ zs ++ (if d /= 0 then '.':fs' else "")
where b = digitsToBits d
n = x' b
ds = show (round_uk ((n*10^d) % 2^b))
(s,ds') = let sgn = head ds == '-' in (sgn, if sgn then tail ds else ds)
ds'' = take (max (d+1length ds') 0) (repeat '0') ++ ds'
(zs,fs) = splitAt (length ds'' d) ds''
fs' = case reverse $ dropWhile (== '0') $ reverse fs of
"" -> "0"
xs -> xs
digitsToBits :: Int -> Int
digitsToBits d = ceiling (fromIntegral d * (logBase 2.0 10.0 :: Double)) + 4
digits :: Int
digits = 40
instance Read CReal where
readsPrec _p = readSigned readFloat
instance Show CReal where
showsPrec p x = let xs = showCReal digits x in
if head xs == '-' then showParen (p > 6) (showString xs)
else showString xs
sizeinbase :: Integer -> Int -> Int
sizeinbase i b = f (abs i)
where f n = if n <= 1 then 1 else 1 + f (n `div` toInteger b)
floorsqrt :: Integer -> Integer
floorsqrt x = until satisfy improve x
where improve y = floor ((y*y+x) % (2*y))
satisfy y = y*y <= x && x <= (y+1)*(y+1)
round_uk :: Rational -> Integer
round_uk x = floor (x+1 % 2)