module Data.Number.Real (
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
data CReal = CReal { state :: IORef (Nat, DI.Interval),
eval :: Nat -> CReal -> DI.Interval }
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'
}
represent :: (D.Precision -> DI.Interval -> DI.Interval) -> CReal -> CReal
represent f r = make (\n -> f n (eval r n r))
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
lim :: (Nat -> CReal)
-> (Nat -> CReal)
-> 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)
(an, rn) = (eval a n' a, eval r n' r)
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..]
limRec :: CReal
-> (CReal -> Nat -> (CReal, CReal))
-> 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
(ak, rk) = (eval an n an, eval rn n rn)
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)
limRat :: (Nat -> D.Dyadic)
-> (Nat -> D.Dyadic)
-> CReal
limRat an rn = make (\n -> DI.fromBall (B.Ball (an n) (rn n)))
infSum :: (Nat -> CReal)
-> (Nat -> CReal)
-> 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
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
pCompare :: CReal -> CReal -> Nat -> POrdering
pCompare r r' = \n -> DI.compareI (eval r n r) (eval r' n r')
infix 4 <.
(<.) :: CReal -> CReal -> Nat -> PBool
(<.) r r' = \n -> case pCompare r r' n of
Less -> PTrue
Greater -> PFalse
_ -> Indeterminate
infix 4 >.
(>.) :: CReal -> CReal -> Nat -> PBool
(>.) r r' = \n -> case pCompare r r' n of
Less -> PFalse
Greater -> PTrue
_ -> Indeterminate
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 :: Int -> CReal
fromInt i = make $ \_ -> DI.fromBall $ B.Ball (D.fromInt D.Near 32 i) $ 0
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 :: 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 :: Nat -> CReal -> String
toString n r = DI.toString (eval r n r)