Real.hs

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.DyadicInterval as DI
import qualified Data.Number.Ball as B
import qualified Data.Number.Dyadic as D

import Data.Order

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

import Data.IORef(IORef, newIORef, writeIORef, readIORef)
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 

instance Read CReal where
    readsPrec _ s = [(fromString s, "")]

instance Num CReal where
    (+) = represent2 DI.add
    (-) = 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)
    fromInteger = fromDyadic . fromInteger

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.
limRat :: (Nat -> D.Dyadic) -- ^ Sequence of dyadics
          -> (Nat -> D.Dyadic) -- ^ Sequence of (dyadic) error bounds
          -> 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     :: CReal -> Nat -> Either (D.Dyadic, Word) D.Dyadic
approx r k = approx' n
             where approx' :: Nat -> Either (D.Dyadic, Word) D.Dyadic
                   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.Dyadic -> CReal
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)