-- ERA: Exact Real Arithmetic (version 1.0)

--

-- A tolerably efficient and possibly correct implementation of the computable

-- reals using Haskell 1.2.

--

-- David Lester, Department of Computer Science, Manchester University, M13 9PL.

--           (2000-2001)



module Data.Number.CReal(CReal, showCReal) where

import Data.Ratio

import Numeric(readFloat, readSigned)



-- |The 'CReal' type implements (constructive) real numbers.

--

-- Note that the comparison operations on 'CReal' may diverge

-- since it is (by necessity) impossible to implementent them

-- correctly and always terminating.

--

-- This implementation is really David Lester's ERA package.

data CReal = CR (Int -> Integer)



instance Eq  CReal where

  x == y = s' (digitsToBits digits) == 0 where (CR s') = x-y



instance Ord CReal where

  x <= y = s' (digitsToBits digits) <= 0 where (CR s') = x-y

  x <  y = s' (digitsToBits digits) <  0 where (CR s') = x-y

  x >= y = s' (digitsToBits digits) >= 0 where (CR s') = x-y

  x >  y = s' (digitsToBits digits) >  0 where (CR s') = x-y

  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)



-- two useful scaling functions:



div2n :: CReal -> Int -> CReal

div2n (CR x') n = CR (\p -> if p >= n then x' (p-n) else round_uk (x' p % 2^n))



mul2n :: CReal -> Int -> CReal

mul2n (CR x') n = CR (\p -> x' (p+n))



-- transcendental functions (mostly range reductions):



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

            {- 7 < t -}   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 = x-fromInteger 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

            {- n == 7 -}   - 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

            {- n == 7 -}   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

            {- t >   4 -}   piBy2 - atan_dr (recip x)

            where (CR x') = x; t = x' 2

                  xp1 = x+1; xm1 = x-1

  asin x  = if x0 >  0 then pi / 2 - atan (s/x) else

            if x0 == 0 then atan (x/s)                      else

            {- x0 <  0 -}   - 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*n-1))))) 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 takes as arguments:

--   a (rational) list of the coefficients of the power series

--   a function from the desired accuracy to the number of terms needed

--   the argument x



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



-- Some useful constants:



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 (+(m-n)) n

  enumFromThenTo n m e = if m >= n then takeWhile (<= e) $ iterate (+(m-n)) n

                          else takeWhile (>= e) $ iterate (+(m-n)) n



instance Real CReal where

 -- toRational x@(CR x') = x' n % 2^n where n = digitsToBits digits

  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 "CReal.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



-- printing and reading the reals:



-- |The 'showCReal' function connverts a 'CReal' to a 'String'.

showCReal :: Int                -- ^ The number of decimals

          -> CReal              -- ^ The real number

          -> String             -- ^ The resulting 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+1-length 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



-- GMP functions not provided by Haskell



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)