```module Data.Number.Real (
-- | show x will output as much decimalas as
-- a standard IEEE 754 double if possible.

-- | (==) and (/=) should not be used as x == y will diverge if
-- two reals should be equal.

CReal(), Nat, Chain,
PBool (..),
min, max,
lim, limRec, limRat, infSum, infSumRec,
approx,
pCompare, (<.), (>.), sqrt, exp, log,
fromDyadic, fromInt, fromWord, fromString, toString, toStringDec
) where

import qualified Data.Number.Ball as B

import Data.Order

import Data.Word(Word)
import Prelude hiding (min, max, sqrt, log, exp)

import System.IO.Unsafe(unsafePerformIO)

import Data.Maybe(isNothing, fromMaybe)

import Data.Ratio(numerator, denominator)

type Nat = Word

type Chain = Nat -> DI.Interval

-- | Real number is represented as a chain of dyadic intervals which
-- are neither necessarily nested nor bounded away from 0.
--
-- On n-th stage computations are performed with precision of n bits.
data CReal = CReal { state :: IORef (Nat, DI.Interval),
eval :: Nat -> CReal -> DI.Interval }

{-# INLINE make #-}
make   :: Chain -> CReal
make c = CReal { state = unsafePerformIO \$ newIORef (D.minPrec, c D.minPrec) ,
eval = \n (CReal s _) -> unsafePerformIO \$
do (n', i) <- readIORef s
if n' == n then return i
else do let i' = c n
writeIORef s (n, i')
return i'
}

{-# INLINE represent #-}
represent     :: (D.Precision -> DI.Interval -> DI.Interval) -> CReal -> CReal
represent f r = make (\n -> f n (eval r n r))

{-# INLINE represent2 #-}
represent2        :: (D.Precision -> DI.Interval -> DI.Interval -> DI.Interval)
-> CReal -> CReal -> CReal
represent2 f r r' =  make (\n -> f n (eval r n r) (eval r' n r'))

max :: CReal -> CReal -> CReal
max = represent2 DI.maxI

min :: CReal -> CReal -> CReal
min = represent2 DI.minI

instance Eq CReal where
r /= r' = or (map (isNothing . (\n -> DI.intersect (eval r n r) (eval r' n r'))) [1..])

instance Show CReal where
show = toStringDec 16

readsPrec _ s = [(fromString s, "")]

instance Num CReal where
(-) = represent2 DI.sub
(*) = represent2 DI.mul
negate = represent DI.neg
abs r = max r (negate r)
signum r = make \$ \n -> case DI.compareI (eval r n r) (eval 0 n 0) of
Less    -> DI.fromInt D.minPrec (negate 1)
Greater -> DI.fromInt D.minPrec 1
_       -> DI.fromBall (B.Ball 0 1)

instance Fractional CReal where
(/) = represent2 DI.div
recip r = 1 / r
fromRational r = fromIntegral (numerator r) / fromIntegral (denominator r)

sqrt :: CReal -> CReal
sqrt = represent DI.sqrt

exp :: CReal -> CReal
exp = represent DI.exp

log :: CReal -> CReal
log = represent DI.log

-- | A basic general limit which takes as arguments a sequence of reals and a sequence of
-- error bounds.
lim       :: (Nat -> CReal) -- ^ Sequence
-> (Nat -> CReal) -- ^ Error bounds
-> CReal
lim am rm = make limStage
where limStage n =  foldl1 DI.intersect lst
where lst = take (fromIntegral n) .
map (\k -> let n' = if k < div n 2 then k else n
(a, r) = (am k, rm k) -- get k-th element of the sequence
(an, rn) = (eval a n' a, eval r n' r) -- get the n-th approximation of the k-th element
i = case (an, rn) of
(Just b, Just b') -> DI.fromBall (B.Ball (B.center b) (B.radius b + B.upper_ b'))
_                 -> Nothing
in i) \$ [1..]

-- | Similar to lim, but can sometimes be more convenient for some sequences
limRec      :: CReal -- ^ initial value
-> (CReal -> Nat -> (CReal, CReal)) -- ^ a function which produces a pair, (next element, error estimate)
-- from previous one and location
-> CReal
limRec st f = make limStage
where limStage n = limStage' 1 st (eval st n st)
where limStage' k st' acc =
let (an, rn) = f st' k -- n-th element of the sequence
(ak, rk) = (eval an n an, eval rn n rn) -- k-th approximation
i = case (ak, rk) of
(Just b, Just b') -> DI.fromBall (B.Ball (B.center b) (B.radius b + B.upper_ b'))
_                 -> Nothing
in if k == n then DI.intersect acc i
else limStage' (succ k) an (DI.intersect acc i)

-- | Limit of a sequence of rationals.
-> CReal
limRat an rn = make (\n -> DI.fromBall (B.Ball (an n) (rn n)))

-- | Computes an infinite sum of a series
infSum      :: (Nat -> CReal) -- ^ Sequence of reals
-> (Nat -> CReal) -- ^ Sequence of series remainders
-> CReal
infSum am rm = make partialsum
where partialsum k = psum 1 (eval a0 k a0) Nothing
where psum n acc res =
let (an,rn) = (am n, rm n)
err = eval rn k rn
acc' = DI.add k acc (eval an k an)
(res', p) = case (acc', err) of
(Just  b, Just b') ->
let (cac,rac) = (B.center b, B.radius b)
(ler, uer) = (B.lower_ b', B.upper_ b')
in (DI.intersect res (DI.fromBall (B.Ball cac (rac + uer))), rac <= ler)
(Nothing, _) -> (Nothing, False)
(_, Nothing) -> (res, True)
in if p then psum (succ n) acc' res'
else res'
a0 = am 0

-- | Similar to infSum but can sometimes be more convenient
-- Second argument is a_0
infSumRec      :: CReal
-> (CReal -> Nat -> (CReal, CReal))
-> CReal
infSumRec st f = make partialsum
where partialsum k = psum 1 (eval a0 k a0) Nothing a0
where psum n acc res t =
let (an, rn) = f t n
err = eval rn k rn
acc' = DI.add k acc (eval an k an)
(res', p) = case (acc', err) of
(Just  b, Just b') ->
let (cac,rac) = (B.center b, B.radius b)
(ler, uer) = (B.lower_ b', B.upper_ b')
in (DI.intersect res (DI.fromBall (B.Ball cac (rac + uer))), rac <= ler)
(Nothing, _) -> (Nothing, False)
(_, Nothing) -> (res, True)
in if p then psum (succ n) acc' res' an
else res'
a0 = st

-- comparison functions

-- | @ pCompare x y @ returns a function @ Nat -> POrdering @ which
-- when applied to some @ n @ computes approximates with precision @ n @
-- and then compares the resulting intervals
pCompare      :: CReal -> CReal -> Nat -> POrdering
pCompare r r' = \n -> DI.compareI (eval r n r) (eval r' n r')

-- | @ x \<. y @ is a function @ Nat -> PBool @ which, when
-- applied to some @ n @, computes the approximation with precision @ n @
-- and then compares the intervals. If intervals are disjoint then result is
-- either PTrue or PFalse, otherwise result is Indeterminate.
infix 4 <.
(<.) :: CReal -> CReal -> Nat -> PBool
(<.) r r' = \n -> case pCompare r r' n of
Less    -> PTrue
Greater -> PFalse
_       -> Indeterminate

-- | Similar to (<.)
infix 4 >.
(>.) :: CReal -> CReal -> Nat -> PBool
(>.) r r' = \n -> case pCompare r r' n of
Less    -> PFalse
Greater -> PTrue
_       -> Indeterminate

-- | @ approx x n @ tries to compute a dyadic approximation to x so than @ |x - d| <= 10^(-n) @
-- If it succeeds it returns @ Right d @ where d is a dyadic rational, otherwise it returns
-- Left (d, n) where d is a dyadic rational and n is the number of accurate decimal places
--
-- Approx succeeds if result can be computed with precision less than the square of the number
-- of required bits of precision.
approx r k = approx' n
approx' n' | cp >= fromIntegral n = Right c
| threshold = Left (c, floor (logBase 10 2 * fromIntegral cp :: Double))
| otherwise = approx' \$ 2 * n'
where cp = if r' == 0 then fromIntegral n
else let t = negate . D.getExp \$ r'
in if t >= 0 then t else 0
B.Ball c r' = fromMaybe (B.Ball 0 (D.pow2 31)) (eval r n' r)
threshold = n * n < n'
n = ceiling ((logBase 2 10 :: Double) * fromIntegral k) + 1

fromDyadic d = make \$ \_ -> DI.fromBall (B.Ball d \$ 0)

-- | fromInt should be preferred over fromIntegral where applicable
fromInt   :: Int -> CReal
fromInt i = make \$ \_ -> DI.fromBall \$ B.Ball (D.fromInt D.Near 32 i) \$ 0

-- | fromWord should be preferred over fromIntegral where applicable
fromWord   :: Word -> CReal
fromWord i = make \$ \_ -> DI.fromBall \$ B.Ball (D.fromWord D.Near 32 i) \$ 0

fromString   :: String -> CReal
fromString s = make (\_ -> let l = length s
n = ceiling (logBase 2 10 * fromIntegral (if elem '.' s then pred l else l) :: Double)
cen = D.fromString s n 10 in
DI.fromBall (B.Ball cen 0))

-- | toStringDec tries to compute the result to the number of specified significand digits
toStringDec     :: Nat -> CReal ->  String
toStringDec n r = inf ++ s
where (inf, s) = case approx r n of
Right d    -> ("",D.toString n d)
Left (d,k) -> ("Could not compute to desired accuracy, only to " ++ show k ++ " significand digits : ",
D.toString k d)

-- | toString computes the result with specified precision.
toString     ::  Nat -> CReal -> String
toString n r = DI.toString (eval r n r)
```